diff --git a/DCHdigi/include/AlgData.h b/DCHdigi/include/AlgData.h index 5289459..d43f2b3 100644 --- a/DCHdigi/include/AlgData.h +++ b/DCHdigi/include/AlgData.h @@ -3,591 +3,535 @@ #include -#include "TTree.h" +#include +#include +#include +#include "Math/InterpolationTypes.h" +#include "Math/Interpolator.h" +#include "TCanvas.h" #include "TChain.h" #include "TF1.h" -#include "TCanvas.h" -#include "TGraph.h" #include "TFile.h" -#include -#include "Math/Interpolator.h" -#include "Math/InterpolationTypes.h" -#include -#include - - +#include "TGraph.h" +#include "TTree.h" -class AlgData{ +class AlgData { public: - inline AlgData(); - inline ~AlgData(); + inline AlgData(); + inline ~AlgData(); public: - inline void read_file(TString dataAlg); - inline TF1* read_graph(TString dataAlg, TString cvName, TString fitName); - - - inline TF1* get_fit(); - inline TFormula* get_formula(); - inline double get_Fitvalue(double betagamma) const; - inline ROOT::Math::Interpolator* Interpvalue(std::vector bg, std::vector Ydata,int Intertype); - inline void Load_file(TString dataAlg); - inline void Load_interp(); - inline double get_MPVExtra (double betagamma) const; - inline double get_SgmExtra (double betagamma) const; - inline double get_MeanExtra1 (double betagamma) const; - inline double get_SgmExtra1 (double betagamma) const; - inline double get_SlopeExtra1 (double betagamma)const; - inline double get_FracExtra1(double betagamma) const; - inline double get_FfracExtra(double betagamma) ; - inline double get_maxEx0(double betagamma); - inline double get_maxExSlp(); - inline double get_ExSgmlep(); - inline double get_ExSgmhad(); - inline double get_Ffrac(double betagamma); - inline double get_Fmpv1(double betagamma); - inline double get_Fsgm1(double betagamma); - inline double get_Fmpv2(double betagamma); - inline double get_Fsgm2(double betagamma); - inline double get_ClSzCorrInt(double betagamma); - inline double get_ClSzCorrSlp(double betagamma); - inline std::vector get_ClSzCorrpmean(double betagamma); - inline std::vector get_ClSzCorrpsgm(double betagamma); - - inline std::vector get_ClSzCorrdgmean(double betagamma); - inline std::vector get_ClSzCorrdgsgm(double betagamma); - - inline std::vector get_ClSzCorrdglfrac(double betagamma); - inline std::vector get_ClSzCorrdglmpvl(double betagamma); - inline std::vector get_ClSzCorrdglsgml(double betagamma); - inline std::vector get_ClSzCorrdglmeang(double betagamma); - inline std::vector get_ClSzCorrdglsgmg(double betagamma); - + inline void read_file(TString dataAlg); + inline TF1* read_graph(TString dataAlg, TString cvName, TString fitName); + + inline TF1* get_fit(); + inline TFormula* get_formula(); + inline double get_Fitvalue(double betagamma) const; + inline ROOT::Math::Interpolator* Interpvalue(std::vector bg, std::vector Ydata, int Intertype); + inline void Load_file(TString dataAlg); + inline void Load_interp(); + inline double get_MPVExtra(double betagamma) const; + inline double get_SgmExtra(double betagamma) const; + inline double get_MeanExtra1(double betagamma) const; + inline double get_SgmExtra1(double betagamma) const; + inline double get_SlopeExtra1(double betagamma) const; + inline double get_FracExtra1(double betagamma) const; + inline double get_FfracExtra(double betagamma); + inline double get_maxEx0(double betagamma); + inline double get_maxExSlp(); + inline double get_ExSgmlep(); + inline double get_ExSgmhad(); + inline double get_Ffrac(double betagamma); + inline double get_Fmpv1(double betagamma); + inline double get_Fsgm1(double betagamma); + inline double get_Fmpv2(double betagamma); + inline double get_Fsgm2(double betagamma); + inline double get_ClSzCorrInt(double betagamma); + inline double get_ClSzCorrSlp(double betagamma); + inline std::vector get_ClSzCorrpmean(double betagamma); + inline std::vector get_ClSzCorrpsgm(double betagamma); + + inline std::vector get_ClSzCorrdgmean(double betagamma); + inline std::vector get_ClSzCorrdgsgm(double betagamma); + + inline std::vector get_ClSzCorrdglfrac(double betagamma); + inline std::vector get_ClSzCorrdglmpvl(double betagamma); + inline std::vector get_ClSzCorrdglsgml(double betagamma); + inline std::vector get_ClSzCorrdglmeang(double betagamma); + inline std::vector get_ClSzCorrdglsgmg(double betagamma); private: - TString DataAlg; - TChain *data1; - TChain *datalep; - TChain *datahad; - TFile *file; - - double bgT,bgTlep,bgThad,tmpmaxEx0Tot,tmpErrmaxEx0Tot,tmpmaxExSlpTot,tmpErrmaxExSlpTot,tmpExSgmTotlep,tmpExSgmTothad; - double tmptFfracTot,tmptFerrfracTot,tmptFmpv1Tot,tmptFerrmpv1Tot,tmptFsgm1Tot,tmptFerrsgm1Tot,tmptFmpv2Tot,tmptFerrmpv2Tot,tmptFsgm2Tot,tmptFerrsgm2Tot; - double tmpCorrIntTot,tmpErrCorrIntTot,tmpCorrSlpTot,tmpErrCorrSlpTot; - - std::vector *tmpCorrmeanSliceTot=nullptr; - std::vector *tmpErrCorrmeanSliceTot=nullptr; - std::vector *tmpCorrsgmSliceTot=nullptr; - std::vector *tmpErrCorrsgmSliceTot=nullptr; - - std::vector *tmpCorrdglfracSliceTot=nullptr; - std::vector *tmpErrCorrdglfracSliceTot=nullptr; - std::vector *tmpCorrdglmeangSliceTot=nullptr; - std::vector *tmpErrCorrdglmeangSliceTot=nullptr; - std::vector *tmpCorrdglsgmgSliceTot=nullptr; - std::vector *tmpErrCorrdglsgmgSliceTot=nullptr; - std::vector *tmpCorrdglmpvSliceTot=nullptr; - std::vector *tmpErrCorrdglmpvSliceTot=nullptr; - std::vector *tmpCorrdglsgmlSliceTot=nullptr; - std::vector *tmpErrCorrdglsgmlSliceTot=nullptr; - - std::vector *tmpCorrdgmeanSliceTot=nullptr; - std::vector *tmpErrCorrdgmeanSliceTot=nullptr; - std::vector *tmpCorrdgsgmSliceTot=nullptr; - std::vector *tmpErrCorrdgsgmSliceTot=nullptr; - - //____________________________________________________// - std::vector bgv,bglep,bghad; - std::vector maxEx0Tot,ErrmaxEx0Tot,maxExSlpTot,ErrmaxExSlpTot,ExSgmTotlep,ExSgmTothad; - std::vector tFfracTot, tFerrfracTot, tFmpv1Tot, tFerrmpv1Tot, tFsgm1Tot, tFerrsgm1Tot, tFmpv2Tot, tFerrmpv2Tot, tFsgm2Tot, tFerrsgm2Tot; - std::vector CorrmeanSliceTot[20],ErrCorrmeanSliceTot[20],CorrsgmSliceTot[20],ErrCorrsgmSliceTot[20]; - std::vector CorrIntTot,ErrCorrIntTot,CorrSlpTot,ErrCorrSlpTot; - std::vector CorrdglfracSliceTot[20],ErrCorrdglfracSliceTot[20],CorrdglmeangSliceTot[20],ErrCorrdglmeangSliceTot[20],CorrdglsgmgSliceTot[20],ErrCorrdglsgmgSliceTot[20],CorrdglmpvSliceTot[20],ErrCorrdglmpvSliceTot[20],CorrdglsgmlSliceTot[20],ErrCorrdglsgmlSliceTot[20]; - std::vector CorrdgmeanSliceTot[20],ErrCorrdgmeanSliceTot[20],CorrdgsgmSliceTot[20],ErrCorrdgsgmSliceTot[20]; - - //__________________________________________// - TF1 *Ft; - TF1 *FitMPVExtra; - TF1 *FitSgmExtra; - TF1 *FitMeanExtra1; - TF1 *FitSgmExtra1; - TF1 *FitSlopeExtra1; - TF1 *FitFracExtra1; - - TFormula *FtFormula; - - std::vector Corrmean; - std::vector Corrsgm; - - std::vector Corrdgmean; - std::vector Corrdgsgm; - - std::vector Corrdglfrac; - std::vector Corrdglmeang; - std::vector Corrdglsgmg; - std::vector Corrdglmpv; - std::vector Corrdglsgml; - - std::vector ClSzCorrpmean; - std::vector ClSzCorrpsgm; - std::vector ClSzCorrdgmean; - std::vector ClSzCorrdgsgm; - std::vector ClSzCorrdglfrac; - std::vector ClSzCorrdglmpvl; - std::vector ClSzCorrdglsgml; - std::vector ClSzCorrdglmeang; - std::vector ClSzCorrdglsgmg; - - double maxExSlp; - double ExSgmlep; - double ExSgmhad; - inline void Calc_maxExSlp(); - inline void Calc_ExSgmlep(); - inline void Calc_ExSgmhad(); - - - - //___________________________________________// - - std::map itpm; - ROOT::Math::Interpolator* itp1; - + TString DataAlg; + TChain* data1; + TChain* datalep; + TChain* datahad; + TFile* file; + + double bgT, bgTlep, bgThad, tmpmaxEx0Tot, tmpErrmaxEx0Tot, tmpmaxExSlpTot, tmpErrmaxExSlpTot, tmpExSgmTotlep, + tmpExSgmTothad; + double tmptFfracTot, tmptFerrfracTot, tmptFmpv1Tot, tmptFerrmpv1Tot, tmptFsgm1Tot, tmptFerrsgm1Tot, tmptFmpv2Tot, + tmptFerrmpv2Tot, tmptFsgm2Tot, tmptFerrsgm2Tot; + double tmpCorrIntTot, tmpErrCorrIntTot, tmpCorrSlpTot, tmpErrCorrSlpTot; + + std::vector* tmpCorrmeanSliceTot = nullptr; + std::vector* tmpErrCorrmeanSliceTot = nullptr; + std::vector* tmpCorrsgmSliceTot = nullptr; + std::vector* tmpErrCorrsgmSliceTot = nullptr; + + std::vector* tmpCorrdglfracSliceTot = nullptr; + std::vector* tmpErrCorrdglfracSliceTot = nullptr; + std::vector* tmpCorrdglmeangSliceTot = nullptr; + std::vector* tmpErrCorrdglmeangSliceTot = nullptr; + std::vector* tmpCorrdglsgmgSliceTot = nullptr; + std::vector* tmpErrCorrdglsgmgSliceTot = nullptr; + std::vector* tmpCorrdglmpvSliceTot = nullptr; + std::vector* tmpErrCorrdglmpvSliceTot = nullptr; + std::vector* tmpCorrdglsgmlSliceTot = nullptr; + std::vector* tmpErrCorrdglsgmlSliceTot = nullptr; + + std::vector* tmpCorrdgmeanSliceTot = nullptr; + std::vector* tmpErrCorrdgmeanSliceTot = nullptr; + std::vector* tmpCorrdgsgmSliceTot = nullptr; + std::vector* tmpErrCorrdgsgmSliceTot = nullptr; + + //____________________________________________________// + std::vector bgv, bglep, bghad; + std::vector maxEx0Tot, ErrmaxEx0Tot, maxExSlpTot, ErrmaxExSlpTot, ExSgmTotlep, ExSgmTothad; + std::vector tFfracTot, tFerrfracTot, tFmpv1Tot, tFerrmpv1Tot, tFsgm1Tot, tFerrsgm1Tot, tFmpv2Tot, + tFerrmpv2Tot, tFsgm2Tot, tFerrsgm2Tot; + std::vector CorrmeanSliceTot[20], ErrCorrmeanSliceTot[20], CorrsgmSliceTot[20], ErrCorrsgmSliceTot[20]; + std::vector CorrIntTot, ErrCorrIntTot, CorrSlpTot, ErrCorrSlpTot; + std::vector CorrdglfracSliceTot[20], ErrCorrdglfracSliceTot[20], CorrdglmeangSliceTot[20], + ErrCorrdglmeangSliceTot[20], CorrdglsgmgSliceTot[20], ErrCorrdglsgmgSliceTot[20], CorrdglmpvSliceTot[20], + ErrCorrdglmpvSliceTot[20], CorrdglsgmlSliceTot[20], ErrCorrdglsgmlSliceTot[20]; + std::vector CorrdgmeanSliceTot[20], ErrCorrdgmeanSliceTot[20], CorrdgsgmSliceTot[20], + ErrCorrdgsgmSliceTot[20]; + + //__________________________________________// + TF1* Ft; + TF1* FitMPVExtra; + TF1* FitSgmExtra; + TF1* FitMeanExtra1; + TF1* FitSgmExtra1; + TF1* FitSlopeExtra1; + TF1* FitFracExtra1; + + TFormula* FtFormula; + + std::vector Corrmean; + std::vector Corrsgm; + + std::vector Corrdgmean; + std::vector Corrdgsgm; + + std::vector Corrdglfrac; + std::vector Corrdglmeang; + std::vector Corrdglsgmg; + std::vector Corrdglmpv; + std::vector Corrdglsgml; + + std::vector ClSzCorrpmean; + std::vector ClSzCorrpsgm; + std::vector ClSzCorrdgmean; + std::vector ClSzCorrdgsgm; + std::vector ClSzCorrdglfrac; + std::vector ClSzCorrdglmpvl; + std::vector ClSzCorrdglsgml; + std::vector ClSzCorrdglmeang; + std::vector ClSzCorrdglsgmg; + + double maxExSlp; + double ExSgmlep; + double ExSgmhad; + inline void Calc_maxExSlp(); + inline void Calc_ExSgmlep(); + inline void Calc_ExSgmhad(); + + //___________________________________________// + + std::map itpm; + ROOT::Math::Interpolator* itp1; }; - //_______________________________________________// +//_______________________________________________// + +AlgData::AlgData() { + DataAlg = ""; + data1 = nullptr; + datalep = nullptr; + datahad = nullptr; + file = nullptr; + + tmpCorrmeanSliceTot = nullptr; + tmpErrCorrmeanSliceTot = nullptr; + tmpCorrsgmSliceTot = nullptr; + tmpErrCorrsgmSliceTot = nullptr; + + tmpCorrdglfracSliceTot = nullptr; + tmpErrCorrdglfracSliceTot = nullptr; + tmpCorrdglmeangSliceTot = nullptr; + tmpErrCorrdglmeangSliceTot = nullptr; + tmpCorrdglsgmgSliceTot = nullptr; + tmpErrCorrdglsgmgSliceTot = nullptr; + tmpCorrdglmpvSliceTot = nullptr; + tmpErrCorrdglmpvSliceTot = nullptr; + tmpCorrdglsgmlSliceTot = nullptr; + tmpErrCorrdglsgmlSliceTot = nullptr; + + tmpCorrdgmeanSliceTot = nullptr; + tmpErrCorrdgmeanSliceTot = nullptr; + tmpCorrdgsgmSliceTot = nullptr; + tmpErrCorrdgsgmSliceTot = nullptr; + + Ft = nullptr; + FtFormula = nullptr; + FitMPVExtra = nullptr; + FitSgmExtra = nullptr; + FitMeanExtra1 = nullptr; + FitSgmExtra1 = nullptr; + FitSlopeExtra1 = nullptr; + FitFracExtra1 = nullptr; +} + +AlgData::~AlgData() { + if (file->IsOpen()) + file->Close(); +} + +void AlgData::read_file(TString dataAlg) { + DataAlg = dataAlg; + + std::cout << " ---reading of---" << DataAlg << std::endl; + + data1 = new TChain("DataAlg"); + datalep = new TChain("DataAlglep"); + datahad = new TChain("DataAlghad"); + + if (dataAlg.Contains(".root")) { + data1->Add(dataAlg.Data()); + datalep->Add(dataAlg.Data()); + datahad->Add(dataAlg.Data()); + } + + data1->SetBranchAddress("bgT", &bgT); + data1->SetBranchAddress("tmpmaxEx0Tot", &tmpmaxEx0Tot); + data1->SetBranchAddress("tmpErrmaxEx0Tot", &tmpErrmaxEx0Tot); + data1->SetBranchAddress("tmpmaxExSlpTot", &tmpmaxExSlpTot); + data1->SetBranchAddress("tmpErrmaxExSlpTot", &tmpErrmaxExSlpTot); + data1->SetBranchAddress("tmptFfracTot", &tmptFfracTot); + data1->SetBranchAddress("tFerrfracTot", &tmptFerrfracTot); + data1->SetBranchAddress("tmptFerrfracTot", &tmptFmpv1Tot); + data1->SetBranchAddress("tmptFerrmpv1Tot", &tmptFerrmpv1Tot); + data1->SetBranchAddress("tmptFerrmpv1Tot", &tmptFerrmpv1Tot); + data1->SetBranchAddress("tmptFsgm1Tot", &tmptFsgm1Tot); + data1->SetBranchAddress("tmptFerrsgm1Tot", &tmptFerrsgm1Tot); + data1->SetBranchAddress("tmptFmpv2Tot", &tmptFmpv2Tot); + data1->SetBranchAddress("tmptFerrmpv2Tot", &tmptFerrmpv2Tot); + data1->SetBranchAddress("tmptFsgm2Tot", &tmptFsgm2Tot); + data1->SetBranchAddress("tmptFerrsgm2Tot", &tmptFerrsgm2Tot); + + data1->SetBranchAddress("tmpCorrIntTot", &tmpCorrIntTot); + data1->SetBranchAddress("tmpErrCorrIntTot", &tmpErrCorrIntTot); + data1->SetBranchAddress("tmpCorrSlpTot", &tmpCorrSlpTot); + data1->SetBranchAddress("tmpErrCorrSlpTot", &tmpErrCorrSlpTot); + + data1->SetBranchAddress("tmpCorrmeanSliceTot", &tmpCorrmeanSliceTot); + data1->SetBranchAddress("tmpErrCorrmeanSliceTot", &tmpErrCorrmeanSliceTot); + data1->SetBranchAddress("tmpCorrsgmSliceTot", &tmpCorrsgmSliceTot); + data1->SetBranchAddress("tmpErrCorrsgmSliceTot", &tmpErrCorrsgmSliceTot); + + data1->SetBranchAddress("tmpCorrdgmeanSliceTot", &tmpCorrdgmeanSliceTot); + data1->SetBranchAddress("tmpErrCorrdgmeanSliceTot", &tmpErrCorrdgmeanSliceTot); + data1->SetBranchAddress("tmpCorrdgsgmSliceTot", &tmpCorrdgsgmSliceTot); + data1->SetBranchAddress("tmpErrCorrdgsgmSliceTot", &tmpErrCorrdgsgmSliceTot); + + data1->SetBranchAddress("tmpCorrdglfracSliceTot", &tmpCorrdglfracSliceTot); + data1->SetBranchAddress("tmpErrCorrdglfracSliceTot", &tmpErrCorrdglfracSliceTot); + data1->SetBranchAddress("tmpCorrdglmeangSliceTot", &tmpCorrdglmeangSliceTot); + data1->SetBranchAddress("tmpErrCorrdglmeangSliceTot", &tmpErrCorrdglmeangSliceTot); + data1->SetBranchAddress("tmpCorrdglsgmgSliceTot", &tmpCorrdglsgmgSliceTot); + data1->SetBranchAddress("tmpErrCorrdglsgmgSliceTot", &tmpErrCorrdglsgmgSliceTot); + data1->SetBranchAddress("tmpCorrdglmpvSliceTot", &tmpCorrdglmpvSliceTot); + data1->SetBranchAddress("tmpErrCorrdglmpvSliceTot", &tmpErrCorrdglmpvSliceTot); + data1->SetBranchAddress("tmpCorrdglsgmlSliceTot", &tmpCorrdglsgmlSliceTot); + data1->SetBranchAddress("tmpErrCorrdglsgmlSliceTot", &tmpErrCorrdglsgmlSliceTot); + + datalep->SetBranchAddress("bgTlep", &bgTlep); + datalep->SetBranchAddress("tmpExSgmTotlep", &tmpExSgmTotlep); + + datahad->SetBranchAddress("bgThad", &bgThad); + datahad->SetBranchAddress("tmpExSgmTothad", &tmpExSgmTothad); + + for (int i = 0; i < data1->GetEntries(); i++) { + data1->GetEntry(i); + bgv.push_back(bgT); + // std::cout<size(); ++j) { + CorrmeanSliceTot[j].push_back(tmpCorrmeanSliceTot->at(j)); + ErrCorrmeanSliceTot[j].push_back(tmpErrCorrmeanSliceTot->at(j)); + CorrsgmSliceTot[j].push_back(tmpCorrsgmSliceTot->at(j)); + ErrCorrsgmSliceTot[j].push_back(tmpErrCorrsgmSliceTot->at(j)); + } + + for (unsigned int n = 0; n < tmpCorrdgmeanSliceTot->size(); ++n) { + CorrdgmeanSliceTot[n].push_back(tmpCorrdgmeanSliceTot->at(n)); + ErrCorrdgmeanSliceTot[n].push_back(tmpErrCorrdgmeanSliceTot->at(n)); + CorrdgsgmSliceTot[n].push_back(tmpCorrdgsgmSliceTot->at(n)); + ErrCorrdgsgmSliceTot[n].push_back(tmpErrCorrdgsgmSliceTot->at(n)); + } + + for (unsigned int k = 0; k < tmpCorrdglfracSliceTot->size(); ++k) { + CorrdglfracSliceTot[k].push_back(tmpCorrdglfracSliceTot->at(k)); + ErrCorrdglfracSliceTot[k].push_back(tmpErrCorrdglfracSliceTot->at(k)); + CorrdglmeangSliceTot[k].push_back(tmpCorrdglmeangSliceTot->at(k)); + ErrCorrdglmeangSliceTot[k].push_back(tmpErrCorrdglmeangSliceTot->at(k)); + CorrdglsgmgSliceTot[k].push_back(tmpCorrdglsgmgSliceTot->at(k)); + ErrCorrdglsgmgSliceTot[k].push_back(tmpErrCorrdglsgmgSliceTot->at(k)); + CorrdglmpvSliceTot[k].push_back(tmpCorrdglmpvSliceTot->at(k)); + ErrCorrdglmpvSliceTot[k].push_back(tmpErrCorrdglmpvSliceTot->at(k)); + CorrdglsgmlSliceTot[k].push_back(tmpCorrdglsgmlSliceTot->at(k)); + ErrCorrdglsgmlSliceTot[k].push_back(tmpErrCorrdglsgmlSliceTot->at(k)); + // std::cout<at(k)<GetEntries(); i++) { + datalep->GetEntry(i); + bglep.push_back(bgTlep); + ExSgmTotlep.push_back(tmpExSgmTotlep); + // std::cout<< " bgTlep "<GetEntries(); i++) { + datahad->GetEntry(i); + bghad.push_back(bgThad); + ExSgmTothad.push_back(tmpExSgmTothad); + } +} + +TF1* AlgData::read_graph(TString dataAlg, TString cvName, TString fitName) { + DataAlg = dataAlg; + file = TFile::Open(dataAlg.Data(), "read"); + TCanvas* cv = (TCanvas*)file->Get(cvName.Data()); + TGraph* gr = (TGraph*)cv->GetListOfPrimitives()->FindObject("Graph"); + TF1* ft = (TF1*)gr->GetListOfFunctions()->FindObject(fitName.Data()); + TFormula* ftFormula = (TFormula*)ft->GetFormula(); + FtFormula = ftFormula; + return ft; +} + +void AlgData::Load_file(TString dataAlg) { + read_file(dataAlg.Data()); + FitMPVExtra = read_graph(dataAlg.Data(), "cMPVExtrabgTot", "fit_MPV"); + FitSgmExtra = read_graph(dataAlg.Data(), "cSgmExtrabgTot", "fit_sgmEx"); + FitMeanExtra1 = read_graph(dataAlg.Data(), "cMeanExtra1bgTot", "expeff"); + FitSgmExtra1 = read_graph(dataAlg.Data(), "cSgmExtra1bgTot", "expeffNeg"); + FitSlopeExtra1 = read_graph(dataAlg.Data(), "cSlopeExtra1bgTot", "fit_slp"); + FitFracExtra1 = read_graph(dataAlg.Data(), "cfracbgTot", "fit_frEx"); + Calc_ExSgmhad(); + Calc_ExSgmlep(); + Calc_maxExSlp(); +} + +TF1* AlgData::get_fit() { return Ft; } + +TFormula* AlgData::get_formula() { return FtFormula; } + +double AlgData::get_MPVExtra(double betagamma) const { return FitMPVExtra->Eval(betagamma); } + +double AlgData::get_SgmExtra(double betagamma) const { return FitSgmExtra->Eval(betagamma); } -AlgData::AlgData(){ +double AlgData::get_MeanExtra1(double betagamma) const { return FitMeanExtra1->Eval(betagamma); } - DataAlg=""; - data1=nullptr; - datalep=nullptr; - datahad=nullptr; - file=nullptr; +double AlgData::get_SgmExtra1(double betagamma) const { return FitSgmExtra1->Eval(betagamma); } +double AlgData::get_FracExtra1(double betagamma) const { return FitFracExtra1->Eval(betagamma); } - tmpCorrmeanSliceTot=nullptr; - tmpErrCorrmeanSliceTot=nullptr; - tmpCorrsgmSliceTot=nullptr; - tmpErrCorrsgmSliceTot=nullptr; +double AlgData::get_SlopeExtra1(double betagamma) const { return FitSlopeExtra1->Eval(betagamma); } - tmpCorrdglfracSliceTot=nullptr; - tmpErrCorrdglfracSliceTot=nullptr; - tmpCorrdglmeangSliceTot=nullptr; - tmpErrCorrdglmeangSliceTot=nullptr; - tmpCorrdglsgmgSliceTot=nullptr; - tmpErrCorrdglsgmgSliceTot=nullptr; - tmpCorrdglmpvSliceTot=nullptr; - tmpErrCorrdglmpvSliceTot=nullptr; - tmpCorrdglsgmlSliceTot=nullptr; - tmpErrCorrdglsgmlSliceTot=nullptr; +ROOT::Math::Interpolator* AlgData::Interpvalue(std::vector bg, std::vector Ydata, int Intertype) { + //tipo di interpolazione 0=linear;1=pol;2=cspline;4=akima + return itp1 = new ROOT::Math::Interpolator(bg, Ydata, (ROOT::Math::Interpolation::Type)Intertype); +} + +void AlgData::Load_interp() { + itpm["maxEx0Tot"] = Interpvalue(bgv, maxEx0Tot, 4); + itpm["tFfracTot"] = Interpvalue(bgv, tFfracTot, 4); + itpm["tFmpv1Tot"] = Interpvalue(bgv, tFmpv1Tot, 4); + itpm["tFsgm1Tot"] = Interpvalue(bgv, tFsgm1Tot, 4); + itpm["tFmpv2Tot"] = Interpvalue(bgv, tFmpv2Tot, 4); + itpm["tFsgm2Tot"] = Interpvalue(bgv, tFsgm2Tot, 4); + + itpm["CorrIntTot"] = Interpvalue(bgv, CorrIntTot, 4); + itpm["CorrSlpTot"] = Interpvalue(bgv, CorrSlpTot, 4); + for (unsigned int i = 0; i < tmpCorrmeanSliceTot->size(); ++i) { + TString nameCm = "CorrmeanSliceTot"; + TString nameCs = "CorrsgmSliceTot"; + Corrmean.push_back(nameCm + i); + Corrsgm.push_back(nameCs + i); + itpm[Corrmean.at(i)] = Interpvalue(bgv, CorrmeanSliceTot[i], 4); + itpm[Corrsgm.at(i)] = Interpvalue(bgv, CorrsgmSliceTot[i], 4); + } + + for (unsigned int j = 0; j < tmpCorrdglfracSliceTot->size(); ++j) { + TString nameCdf = "Corrdglfrac"; + TString nameCdm = "Corrdglmeang"; + TString nameCds = "Corrdglsgmg"; + TString nameCdmpv = "Corrdglmpv"; + TString nameCdsgl = "Corrdglsgmgl"; + Corrdglfrac.push_back(nameCdf + j); + Corrdglmeang.push_back(nameCdm + j); + Corrdglsgmg.push_back(nameCds + j); + Corrdglmpv.push_back(nameCdmpv + j); + Corrdglsgml.push_back(nameCdsgl + j); + itpm[Corrdglfrac.at(j)] = Interpvalue(bgv, CorrdglfracSliceTot[j], 4); + itpm[Corrdglmeang.at(j)] = Interpvalue(bgv, CorrdglmeangSliceTot[j], 4); + itpm[Corrdglsgmg.at(j)] = Interpvalue(bgv, CorrdglsgmgSliceTot[j], 4); + itpm[Corrdglmpv.at(j)] = Interpvalue(bgv, CorrdglmpvSliceTot[j], 4); + itpm[Corrdglsgml.at(j)] = Interpvalue(bgv, CorrdglsgmlSliceTot[j], 4); + } + + for (unsigned int n = 0; n < tmpCorrdgmeanSliceTot->size(); ++n) { + TString nameCdmg = "Corrdgmean"; + TString nameCdsgg = "Corrdgsgm"; + Corrdgmean.push_back(nameCdmg + n); + Corrdgsgm.push_back(nameCdsgg + n); + itpm[Corrdgmean.at(n)] = Interpvalue(bgv, CorrdgmeanSliceTot[n], 4); + itpm[Corrdgsgm.at(n)] = Interpvalue(bgv, CorrdgsgmSliceTot[n], 4); + } +} - tmpCorrdgmeanSliceTot=nullptr; - tmpErrCorrdgmeanSliceTot=nullptr; - tmpCorrdgsgmSliceTot=nullptr; - tmpErrCorrdgsgmSliceTot=nullptr; +double AlgData::get_maxEx0(double betagamma) { return itpm["maxEx0Tot"]->Eval(betagamma); } +double AlgData::get_Ffrac(double betagamma) { return itpm["tFfracTot"]->Eval(betagamma); } +double AlgData::get_Fmpv1(double betagamma) { return itpm["tFmpv1Tot"]->Eval(betagamma); } +double AlgData::get_Fsgm1(double betagamma) { return itpm["tFsgm1Tot"]->Eval(betagamma); } +double AlgData::get_Fmpv2(double betagamma) { return itpm["tFmpv2Tot"]->Eval(betagamma); } +double AlgData::get_Fsgm2(double betagamma) { return itpm["tFsgm2Tot"]->Eval(betagamma); } +double AlgData::get_ClSzCorrInt(double betagamma) { return itpm["CorrIntTot"]->Eval(betagamma); } +double AlgData::get_ClSzCorrSlp(double betagamma) { return itpm["CorrSlpTot"]->Eval(betagamma); } + +std::vector AlgData::get_ClSzCorrpmean(double betagamma) { + ClSzCorrpmean.clear(); + // std::cout<<"size "<size()<size(); ++i) { + ClSzCorrpmean.push_back(itpm[Corrmean.at(i)]->Eval(betagamma)); + } + + return ClSzCorrpmean; +} - Ft=nullptr; - FtFormula=nullptr; - FitMPVExtra=nullptr; - FitSgmExtra=nullptr; - FitMeanExtra1=nullptr; - FitSgmExtra1=nullptr; - FitSlopeExtra1=nullptr; - FitFracExtra1=nullptr; +std::vector AlgData::get_ClSzCorrpsgm(double betagamma) { + ClSzCorrpsgm.clear(); + for (unsigned int i = 0; i < tmpCorrmeanSliceTot->size(); ++i) { + ClSzCorrpsgm.push_back(itpm[Corrsgm.at(i)]->Eval(betagamma)); + } + return ClSzCorrpsgm; +} + +std::vector AlgData::get_ClSzCorrdglfrac(double betagamma) { + ClSzCorrdglfrac.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglfrac.push_back(itpm[Corrdglfrac.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglfrac; +} + +std::vector AlgData::get_ClSzCorrdglmeang(double betagamma) { + ClSzCorrdglmeang.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglmeang.push_back(itpm[Corrdglmeang.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglmeang; +} + +std::vector AlgData::get_ClSzCorrdglsgmg(double betagamma) { + ClSzCorrdglsgmg.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglsgmg.push_back(itpm[Corrdglsgmg.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglsgmg; +} + +std::vector AlgData::get_ClSzCorrdglmpvl(double betagamma) { + ClSzCorrdglmpvl.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglmpvl.push_back(itpm[Corrdglmpv.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglmpvl; +} + +std::vector AlgData::get_ClSzCorrdglsgml(double betagamma) { + ClSzCorrdglsgml.clear(); + for (unsigned int i = 0; i < tmpCorrdglfracSliceTot->size(); ++i) { + ClSzCorrdglsgml.push_back(itpm[Corrdglsgml.at(i)]->Eval(betagamma)); + } + return ClSzCorrdglsgml; +} + +std::vector AlgData::get_ClSzCorrdgmean(double betagamma) { + ClSzCorrdgmean.clear(); + for (unsigned int i = 0; i < tmpCorrdgmeanSliceTot->size(); ++i) { + ClSzCorrdgmean.push_back(itpm[Corrdgmean.at(i)]->Eval(betagamma)); + } + return ClSzCorrdgmean; +} + +std::vector AlgData::get_ClSzCorrdgsgm(double betagamma) { + ClSzCorrdgsgm.clear(); + for (unsigned int i = 0; i < tmpCorrdgmeanSliceTot->size(); ++i) { + ClSzCorrdgsgm.push_back(itpm[Corrdgsgm.at(i)]->Eval(betagamma)); + } + return ClSzCorrdgsgm; +} + +double AlgData::get_maxExSlp() { return maxExSlp; } + +void AlgData::Calc_maxExSlp() { + double sum = 0.0; + for (unsigned int i = 0; i < maxExSlpTot.size(); ++i) { + sum += maxExSlpTot[i]; + } + maxExSlp = sum / (double)maxExSlpTot.size(); +} + +double AlgData::get_ExSgmlep() { return ExSgmlep; } + +void AlgData::Calc_ExSgmlep() { + double sum = 0.0; + for (unsigned int i = 0; i < ExSgmTotlep.size(); ++i) { + sum += ExSgmTotlep[i]; + } + ExSgmlep = sum / (double)ExSgmTotlep.size(); +} +double AlgData::get_ExSgmhad() { return ExSgmhad; } +void AlgData::Calc_ExSgmhad() { + double sum = 0.0; + for (unsigned int i = 0; i < ExSgmTothad.size(); ++i) { + sum += ExSgmTothad[i]; + } + ExSgmhad = sum / (double)ExSgmTothad.size(); } - AlgData::~AlgData(){ - if(file->IsOpen())file->Close(); - } - - - void AlgData::read_file(TString dataAlg){ - - DataAlg=dataAlg; - - std::cout<<" ---reading of---"<Add(dataAlg.Data()); - datalep->Add(dataAlg.Data()); - datahad->Add(dataAlg.Data()); - - } - - data1->SetBranchAddress("bgT",&bgT); - data1->SetBranchAddress("tmpmaxEx0Tot",&tmpmaxEx0Tot); - data1->SetBranchAddress("tmpErrmaxEx0Tot",&tmpErrmaxEx0Tot); - data1->SetBranchAddress("tmpmaxExSlpTot",&tmpmaxExSlpTot); - data1->SetBranchAddress("tmpErrmaxExSlpTot",&tmpErrmaxExSlpTot); - data1->SetBranchAddress("tmptFfracTot",&tmptFfracTot); - data1->SetBranchAddress("tFerrfracTot",&tmptFerrfracTot); - data1->SetBranchAddress("tmptFerrfracTot",&tmptFmpv1Tot); - data1->SetBranchAddress("tmptFerrmpv1Tot",&tmptFerrmpv1Tot); - data1->SetBranchAddress("tmptFerrmpv1Tot",&tmptFerrmpv1Tot); - data1->SetBranchAddress("tmptFsgm1Tot",&tmptFsgm1Tot); - data1->SetBranchAddress("tmptFerrsgm1Tot",&tmptFerrsgm1Tot); - data1->SetBranchAddress("tmptFmpv2Tot",&tmptFmpv2Tot); - data1->SetBranchAddress("tmptFerrmpv2Tot",&tmptFerrmpv2Tot); - data1->SetBranchAddress("tmptFsgm2Tot",&tmptFsgm2Tot); - data1->SetBranchAddress("tmptFerrsgm2Tot",&tmptFerrsgm2Tot); - - data1->SetBranchAddress("tmpCorrIntTot",&tmpCorrIntTot); - data1->SetBranchAddress("tmpErrCorrIntTot",&tmpErrCorrIntTot); - data1->SetBranchAddress("tmpCorrSlpTot",&tmpCorrSlpTot); - data1->SetBranchAddress("tmpErrCorrSlpTot",&tmpErrCorrSlpTot); - - data1->SetBranchAddress("tmpCorrmeanSliceTot",&tmpCorrmeanSliceTot); - data1->SetBranchAddress("tmpErrCorrmeanSliceTot",&tmpErrCorrmeanSliceTot); - data1->SetBranchAddress("tmpCorrsgmSliceTot",&tmpCorrsgmSliceTot); - data1->SetBranchAddress("tmpErrCorrsgmSliceTot",&tmpErrCorrsgmSliceTot); - - data1->SetBranchAddress("tmpCorrdgmeanSliceTot",&tmpCorrdgmeanSliceTot); - data1->SetBranchAddress("tmpErrCorrdgmeanSliceTot",&tmpErrCorrdgmeanSliceTot); - data1->SetBranchAddress("tmpCorrdgsgmSliceTot",&tmpCorrdgsgmSliceTot); - data1->SetBranchAddress("tmpErrCorrdgsgmSliceTot",&tmpErrCorrdgsgmSliceTot); - - data1->SetBranchAddress("tmpCorrdglfracSliceTot",&tmpCorrdglfracSliceTot); - data1->SetBranchAddress("tmpErrCorrdglfracSliceTot",&tmpErrCorrdglfracSliceTot); - data1->SetBranchAddress("tmpCorrdglmeangSliceTot",&tmpCorrdglmeangSliceTot); - data1->SetBranchAddress("tmpErrCorrdglmeangSliceTot",&tmpErrCorrdglmeangSliceTot); - data1->SetBranchAddress("tmpCorrdglsgmgSliceTot",&tmpCorrdglsgmgSliceTot); - data1->SetBranchAddress("tmpErrCorrdglsgmgSliceTot",&tmpErrCorrdglsgmgSliceTot); - data1->SetBranchAddress("tmpCorrdglmpvSliceTot",&tmpCorrdglmpvSliceTot); - data1->SetBranchAddress("tmpErrCorrdglmpvSliceTot",&tmpErrCorrdglmpvSliceTot); - data1->SetBranchAddress("tmpCorrdglsgmlSliceTot",&tmpCorrdglsgmlSliceTot); - data1->SetBranchAddress("tmpErrCorrdglsgmlSliceTot",&tmpErrCorrdglsgmlSliceTot); - - datalep->SetBranchAddress("bgTlep",&bgTlep); - datalep->SetBranchAddress("tmpExSgmTotlep",&tmpExSgmTotlep); - - datahad->SetBranchAddress("bgThad",&bgThad); - datahad->SetBranchAddress("tmpExSgmTothad",&tmpExSgmTothad); - - - for (int i=0;iGetEntries();i++) { - data1->GetEntry(i); - bgv.push_back(bgT); -// std::cout<size();++j){ - CorrmeanSliceTot[j].push_back(tmpCorrmeanSliceTot->at(j)); - ErrCorrmeanSliceTot[j].push_back(tmpErrCorrmeanSliceTot->at(j)); - CorrsgmSliceTot[j].push_back(tmpCorrsgmSliceTot->at(j)); - ErrCorrsgmSliceTot[j].push_back(tmpErrCorrsgmSliceTot->at(j)); - } - - for( unsigned int n=0;nsize();++n){ - CorrdgmeanSliceTot[n].push_back(tmpCorrdgmeanSliceTot->at(n)); - ErrCorrdgmeanSliceTot[n].push_back(tmpErrCorrdgmeanSliceTot->at(n)); - CorrdgsgmSliceTot[n].push_back(tmpCorrdgsgmSliceTot->at(n)); - ErrCorrdgsgmSliceTot[n].push_back(tmpErrCorrdgsgmSliceTot->at(n)); - } - - for(unsigned int k=0;ksize();++k){ - CorrdglfracSliceTot[k].push_back(tmpCorrdglfracSliceTot->at(k)); - ErrCorrdglfracSliceTot[k].push_back(tmpErrCorrdglfracSliceTot->at(k)); - CorrdglmeangSliceTot[k].push_back(tmpCorrdglmeangSliceTot->at(k)); - ErrCorrdglmeangSliceTot[k].push_back(tmpErrCorrdglmeangSliceTot->at(k)); - CorrdglsgmgSliceTot[k].push_back(tmpCorrdglsgmgSliceTot->at(k)); - ErrCorrdglsgmgSliceTot[k].push_back(tmpErrCorrdglsgmgSliceTot->at(k)); - CorrdglmpvSliceTot[k].push_back(tmpCorrdglmpvSliceTot->at(k)); - ErrCorrdglmpvSliceTot[k].push_back(tmpErrCorrdglmpvSliceTot->at(k)); - CorrdglsgmlSliceTot[k].push_back(tmpCorrdglsgmlSliceTot->at(k)); - ErrCorrdglsgmlSliceTot[k].push_back(tmpErrCorrdglsgmlSliceTot->at(k)); -// std::cout<at(k)<GetEntries();i++) { - datalep->GetEntry(i); - bglep.push_back(bgTlep); - ExSgmTotlep.push_back(tmpExSgmTotlep); -// std::cout<< " bgTlep "<GetEntries();i++) { - datahad->GetEntry(i); - bghad.push_back(bgThad); - ExSgmTothad.push_back(tmpExSgmTothad); - } - - } - - - - TF1* AlgData::read_graph(TString dataAlg, TString cvName, TString fitName){ - DataAlg=dataAlg; - file = TFile::Open(dataAlg.Data(),"read"); - TCanvas *cv=(TCanvas *)file->Get(cvName.Data()); - TGraph *gr=(TGraph *)cv->GetListOfPrimitives()->FindObject("Graph"); - TF1* ft=(TF1*)gr->GetListOfFunctions()->FindObject(fitName.Data()); - TFormula *ftFormula=(TFormula *)ft->GetFormula(); - FtFormula=ftFormula; - return ft; - } - - void AlgData::Load_file(TString dataAlg){ - read_file(dataAlg.Data()); - FitMPVExtra=read_graph(dataAlg.Data(), "cMPVExtrabgTot", "fit_MPV"); - FitSgmExtra=read_graph(dataAlg.Data(), "cSgmExtrabgTot", "fit_sgmEx"); - FitMeanExtra1=read_graph(dataAlg.Data(), "cMeanExtra1bgTot", "expeff"); - FitSgmExtra1=read_graph(dataAlg.Data(), "cSgmExtra1bgTot", "expeffNeg"); - FitSlopeExtra1=read_graph(dataAlg.Data(), "cSlopeExtra1bgTot", "fit_slp"); - FitFracExtra1=read_graph(dataAlg.Data(), "cfracbgTot", "fit_frEx"); - Calc_ExSgmhad(); - Calc_ExSgmlep(); - Calc_maxExSlp(); - } - - TF1* AlgData::get_fit(){ - return Ft; - } - - TFormula * AlgData::get_formula(){ - return FtFormula; - } - - double AlgData::get_MPVExtra (double betagamma) const{ - return FitMPVExtra->Eval(betagamma); - } - - double AlgData::get_SgmExtra(double betagamma) const{ - return FitSgmExtra->Eval(betagamma); - } - - double AlgData::get_MeanExtra1(double betagamma) const{ - return FitMeanExtra1->Eval(betagamma); - } - - double AlgData::get_SgmExtra1(double betagamma) const{ - return FitSgmExtra1->Eval(betagamma); - } - - double AlgData::get_FracExtra1(double betagamma) const{ - return FitFracExtra1->Eval(betagamma); - } - - double AlgData::get_SlopeExtra1(double betagamma) const{ - return FitSlopeExtra1->Eval(betagamma); - } - - - ROOT::Math::Interpolator* AlgData::Interpvalue(std::vectorbg,std::vectorYdata,int Intertype) { - //tipo di interpolazione 0=linear;1=pol;2=cspline;4=akima - return itp1=new ROOT::Math::Interpolator(bg,Ydata,(ROOT::Math::Interpolation::Type)Intertype); - } - - void AlgData::Load_interp(){ - itpm["maxEx0Tot"]=Interpvalue(bgv,maxEx0Tot,4); - itpm["tFfracTot"]=Interpvalue(bgv,tFfracTot,4); - itpm["tFmpv1Tot"]=Interpvalue(bgv,tFmpv1Tot,4); - itpm["tFsgm1Tot"]=Interpvalue(bgv,tFsgm1Tot,4); - itpm["tFmpv2Tot"]=Interpvalue(bgv,tFmpv2Tot,4); - itpm["tFsgm2Tot"]=Interpvalue(bgv,tFsgm2Tot,4); - - itpm["CorrIntTot"]=Interpvalue(bgv,CorrIntTot,4); - itpm["CorrSlpTot"]=Interpvalue(bgv,CorrSlpTot,4); - for(unsigned int i=0;isize();++i){ - TString nameCm="CorrmeanSliceTot"; - TString nameCs="CorrsgmSliceTot"; - Corrmean.push_back(nameCm+i); - Corrsgm.push_back(nameCs+i); - itpm[Corrmean.at(i)]=Interpvalue(bgv, CorrmeanSliceTot[i], 4); - itpm[Corrsgm.at(i)]=Interpvalue(bgv, CorrsgmSliceTot[i], 4); - - } - - for(unsigned int j=0;jsize();++j){ - TString nameCdf="Corrdglfrac"; - TString nameCdm="Corrdglmeang"; - TString nameCds="Corrdglsgmg"; - TString nameCdmpv="Corrdglmpv"; - TString nameCdsgl="Corrdglsgmgl"; - Corrdglfrac.push_back(nameCdf+j); - Corrdglmeang.push_back(nameCdm+j); - Corrdglsgmg.push_back(nameCds+j); - Corrdglmpv.push_back(nameCdmpv+j); - Corrdglsgml.push_back(nameCdsgl+j); - itpm[Corrdglfrac.at(j)]=Interpvalue(bgv, CorrdglfracSliceTot[j], 4); - itpm[Corrdglmeang.at(j)]=Interpvalue(bgv, CorrdglmeangSliceTot[j], 4); - itpm[Corrdglsgmg.at(j)]=Interpvalue(bgv, CorrdglsgmgSliceTot[j], 4); - itpm[Corrdglmpv.at(j)]=Interpvalue(bgv, CorrdglmpvSliceTot[j], 4); - itpm[Corrdglsgml.at(j)]=Interpvalue(bgv, CorrdglsgmlSliceTot[j], 4); - } - - for(unsigned int n=0;nsize();++n){ - TString nameCdmg="Corrdgmean"; - TString nameCdsgg="Corrdgsgm"; - Corrdgmean.push_back(nameCdmg+n); - Corrdgsgm.push_back(nameCdsgg+n); - itpm[Corrdgmean.at(n)]=Interpvalue(bgv, CorrdgmeanSliceTot[n], 4); - itpm[Corrdgsgm.at(n)]=Interpvalue(bgv, CorrdgsgmSliceTot[n], 4); - } - } - - double AlgData::get_maxEx0(double betagamma) { - return itpm["maxEx0Tot"]->Eval(betagamma); - } - double AlgData::get_Ffrac(double betagamma) { - return itpm["tFfracTot"]->Eval(betagamma); - } - double AlgData::get_Fmpv1(double betagamma) { - return itpm["tFmpv1Tot"]->Eval(betagamma); - } - double AlgData::get_Fsgm1(double betagamma) { - return itpm["tFsgm1Tot"]->Eval(betagamma); - } - double AlgData::get_Fmpv2(double betagamma) { - return itpm["tFmpv2Tot"]->Eval(betagamma); - } - double AlgData::get_Fsgm2(double betagamma) { - return itpm["tFsgm2Tot"]->Eval(betagamma); - } - double AlgData::get_ClSzCorrInt(double betagamma){ - return itpm["CorrIntTot"]->Eval(betagamma); - } - double AlgData::get_ClSzCorrSlp(double betagamma){ - return itpm["CorrSlpTot"]->Eval(betagamma); - } - - std::vector AlgData::get_ClSzCorrpmean(double betagamma){ - ClSzCorrpmean.clear(); -// std::cout<<"size "<size()<size();++i){ - ClSzCorrpmean.push_back(itpm[Corrmean.at(i)]->Eval(betagamma)); - } - - return ClSzCorrpmean; - } - - std::vector AlgData::get_ClSzCorrpsgm(double betagamma){ - ClSzCorrpsgm.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrpsgm.push_back(itpm[Corrsgm.at(i)]->Eval(betagamma)); - } - return ClSzCorrpsgm; - } - - std::vector AlgData::get_ClSzCorrdglfrac(double betagamma){ - ClSzCorrdglfrac.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdglfrac.push_back(itpm[Corrdglfrac.at(i)]->Eval(betagamma)); - } - return ClSzCorrdglfrac; - } - - std::vector AlgData::get_ClSzCorrdglmeang(double betagamma){ - ClSzCorrdglmeang.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdglmeang.push_back(itpm[Corrdglmeang.at(i)]->Eval(betagamma)); - } - return ClSzCorrdglmeang; - } - - std::vector AlgData::get_ClSzCorrdglsgmg(double betagamma){ - ClSzCorrdglsgmg.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdglsgmg.push_back(itpm[Corrdglsgmg.at(i)]->Eval(betagamma)); - } - return ClSzCorrdglsgmg; - } - - std::vector AlgData::get_ClSzCorrdglmpvl(double betagamma){ - ClSzCorrdglmpvl.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdglmpvl.push_back(itpm[Corrdglmpv.at(i)]->Eval(betagamma)); - } - return ClSzCorrdglmpvl; - } - - std::vector AlgData::get_ClSzCorrdglsgml(double betagamma){ - ClSzCorrdglsgml.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdglsgml.push_back(itpm[Corrdglsgml.at(i)]->Eval(betagamma)); - } - return ClSzCorrdglsgml; - } - - std::vector AlgData::get_ClSzCorrdgmean(double betagamma){ - ClSzCorrdgmean.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdgmean.push_back(itpm[Corrdgmean.at(i)]->Eval(betagamma)); - } - return ClSzCorrdgmean; - } - - std::vector AlgData::get_ClSzCorrdgsgm(double betagamma){ - ClSzCorrdgsgm.clear(); - for(unsigned int i=0;isize();++i){ - ClSzCorrdgsgm.push_back(itpm[Corrdgsgm.at(i)]->Eval(betagamma)); - } - return ClSzCorrdgsgm; - } - - - double AlgData::get_maxExSlp(){ - return maxExSlp; - } - - void AlgData::Calc_maxExSlp(){ - double sum=0.0; - for(unsigned int i=0;i #include +#include // data extension for the DCH v2 #include "DDRec/DCH_info.h" // ROOT headers -#include "TVector3.h" -#include "TRandom3.h" #include "TFile.h" #include "TH1D.h" +#include "TRandom3.h" +#include "TVector3.h" // Class developed by Walaa for the CLS #include "AlgData.h" @@ -78,20 +76,17 @@ constexpr double MM_TO_CM = 0.1; struct DCHdigi final : k4FWCore::MultiTransformer< - std::tuple( + std::tuple( const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&)> { DCHdigi(const std::string& name, ISvcLocator* svcLoc); StatusCode initialize() override; StatusCode finalize() override; - std::tuple operator()( - const edm4hep::SimTrackerHitCollection& , - const edm4hep::EventHeaderCollection& ) const override; + std::tuple + operator()(const edm4hep::SimTrackerHitCollection&, const edm4hep::EventHeaderCollection&) const override; private: - - //------------------------------------------------------------------ // machinery for geometry @@ -103,29 +98,30 @@ struct DCHdigi final Gaudi::Property m_DCH_name{this, "DCH_name", "DCH_v2", "Name of the Drift Chamber detector"}; /// Pointer to the geometry service - SmartIF m_geoSvc; + SmartIF m_geoSvc; /// Decoder for the cellID dd4hep::DDSegmentation::BitFieldCoder* m_decoder; /// Pointer to drift chamber data extension - dd4hep::rec::DCH_info * dch_data = {nullptr}; + dd4hep::rec::DCH_info* dch_data = {nullptr}; //------------------------------------------------------------------ // machinery for smearing the position /// along the sense wire position resolution in mm - Gaudi::Property m_z_resolution{this, "zResolution_mm", 1.0, - "Spatial resolution in the z direction (from reading out the wires at both sides) in mm. Default 1 mm."}; + Gaudi::Property m_z_resolution{ + this, "zResolution_mm", 1.0, + "Spatial resolution in the z direction (from reading out the wires at both sides) in mm. Default 1 mm."}; /// xy resolution in mm - Gaudi::Property m_xy_resolution{this, "xyResolution_mm", 0.1, "Spatial resolution in the xy direction in mm. Default 0.1 mm."}; + Gaudi::Property m_xy_resolution{this, "xyResolution_mm", 0.1, + "Spatial resolution in the xy direction in mm. Default 0.1 mm."}; /// create seed using the uid - SmartIF m_uidSvc; + SmartIF m_uidSvc; /// use thread local engine from C++ standard inline static thread_local std::mt19937_64 m_engine; - void PrepareRandomEngine(const edm4hep::EventHeaderCollection& headers) const; - + void PrepareRandomEngine(const edm4hep::EventHeaderCollection& headers) const; // Operator std::normal_distribution::operator()(Generator& g) is a non-const member function and thus cannot be called for a constant object. So we defined the distribution as mutable. // Gaussian random number generator used for the smearing of the z position, in cm! @@ -137,7 +133,7 @@ struct DCHdigi final //------------------------------------------------------------------ // ancillary functions - bool IsFileGood(std::string & ifilename) const {return std::ifstream(ifilename).good(); } + bool IsFileGood(std::string& ifilename) const { return std::ifstream(ifilename).good(); } /// Print algorithm configuration void PrintConfiguration(std::ostream& io); @@ -146,60 +142,60 @@ struct DCHdigi final void ThrowException(std::string s) const; int CalculateLayerFromCellID(dd4hep::DDSegmentation::CellID id) const { - return m_decoder->get ( id,"layer" ) + dch_data->nlayersPerSuperlayer*m_decoder->get ( id,"superlayer" ) + 1; + return m_decoder->get(id, "layer") + dch_data->nlayersPerSuperlayer * m_decoder->get(id, "superlayer") + 1; } - int CalculateNphiFromCellID(dd4hep::DDSegmentation::CellID id) const { - return m_decoder->get ( id,"nphi" ); - } + int CalculateNphiFromCellID(dd4hep::DDSegmentation::CellID id) const { return m_decoder->get(id, "nphi"); } - TVector3 Convert_EDM4hepVector_to_TVector3(const edm4hep::Vector3d & v, double scale) const { - return {v[0]*scale,v[1]*scale,v[2]*scale}; + TVector3 Convert_EDM4hepVector_to_TVector3(const edm4hep::Vector3d& v, double scale) const { + return {v[0] * scale, v[1] * scale, v[2] * scale}; }; - edm4hep::Vector3d Convert_TVector3_to_EDM4hepVector(const TVector3 & v, double scale) const{ - return {v.x()*scale,v.y()*scale,v.z()*scale}; + edm4hep::Vector3d Convert_TVector3_to_EDM4hepVector(const TVector3& v, double scale) const { + return {v.x() * scale, v.y() * scale, v.z() * scale}; }; // the following functions should be upstreamed to the data extension at DD4hep // to avoid code duplication and keep it centralized - TVector3 Calculate_hitpos_to_wire_vector(int ilayer, int nphi, const TVector3 & hit_position /*in cm*/) const; - TVector3 Calculate_wire_vector_ez (int ilayer, int nphi) const; - TVector3 Calculate_wire_z0_point (int ilayer, int nphi) const; - double Calculate_wire_phi_z0 (int ilayer, int nphi) const; - + TVector3 Calculate_hitpos_to_wire_vector(int ilayer, int nphi, const TVector3& hit_position /*in cm*/) const; + TVector3 Calculate_wire_vector_ez(int ilayer, int nphi) const; + TVector3 Calculate_wire_z0_point(int ilayer, int nphi) const; + double Calculate_wire_phi_z0(int ilayer, int nphi) const; //------------------------------------------------------------------ // cluster calculation, developed by Walaa /// file with distributions to be sampled - Gaudi::Property m_fileDataAlg{this, "fileDataAlg", "/eos/project/f/fccsw-web/www/filesForSimDigiReco/IDEA/DataAlgFORGEANT.root", "ROOT file with cluster size distributions"}; + Gaudi::Property m_fileDataAlg{ + this, "fileDataAlg", "/eos/project/f/fccsw-web/www/filesForSimDigiReco/IDEA/DataAlgFORGEANT.root", + "ROOT file with cluster size distributions"}; /// pointer to wrapper class, which contains the cluster size and number distributions - AlgData * flData; + AlgData* flData; /// code developed by Walaa for calculating number of clusters and cluster size - std::pair CalculateClusters(const edm4hep::SimTrackerHit & input_sim_hit) const; + std::pair CalculateClusters(const edm4hep::SimTrackerHit& input_sim_hit) const; //------------------------------------------------------------------ // debug information /// Flag to create output file with debug histgrams - Gaudi::Property m_create_debug_histos{this, "create_debug_histograms", false, "Create output file with histograms for debugging"}; - + Gaudi::Property m_create_debug_histos{this, "create_debug_histograms", false, + "Create output file with histograms for debugging"}; + /// name for the file that will contain the histograms for debugging - Gaudi::Property m_out_debug_filename{this, "out_debug_filename", "dch_digi_alg_debug.root", "name for the file that will contain the histograms for debugging"}; + Gaudi::Property m_out_debug_filename{this, "out_debug_filename", "dch_digi_alg_debug.root", + "name for the file that will contain the histograms for debugging"}; /// histogram to store distance from hit position to the wire - TH1D * hDpw; + TH1D* hDpw; /// histogram to store distance from hit projection to the wire (should be zero) - TH1D * hDww; + TH1D* hDww; /// histogram to store smearing along the wire - TH1D * hSz; + TH1D* hSz; /// histogram to store smearing perpendicular the wire - TH1D * hSxy; - + TH1D* hSxy; }; DECLARE_COMPONENT(DCHdigi); diff --git a/DCHdigi/src/DCHdigi.cpp b/DCHdigi/src/DCHdigi.cpp index a83e611..e69de29 100644 --- a/DCHdigi/src/DCHdigi.cpp +++ b/DCHdigi/src/DCHdigi.cpp @@ -1,664 +0,0 @@ -#include "DCHdigi.h" - -// STL -#include -#include - -#include "extension/MutableDriftChamberDigiV2.h" - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// DCHdigi constructor //////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -// -- KeyValues("name of the variable that holds the name of the collection exposed in the python steering file", {"default name for the collection"}), -DCHdigi::DCHdigi(const std::string& name, ISvcLocator* svcLoc) -: MultiTransformer(name, svcLoc, - { - KeyValues("DCH_simhits", {""}), - KeyValues("HeaderName", {"EventHeader"}), - }, - { - KeyValues("DCH_DigiCollection", {"DCH_DigiCollection"}), - KeyValues("DCH_DigiSimAssociationCollection", {"DCH_DigiSimAssociationCollection"}) - } - ) -{ - m_geoSvc = serviceLocator()->service(m_geoSvcName); - m_uidSvc = serviceLocator()->service(m_uidSvcName); -} - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// initialize //////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -StatusCode DCHdigi::initialize() { - - if (!m_uidSvc) - ThrowException( "Unable to get UniqueIDGenSvc" ); - - m_gauss_z_cm = std::normal_distribution(0., m_z_resolution.value()*MM_TO_CM ); - m_gauss_xy_cm = std::normal_distribution(0., m_xy_resolution.value()*MM_TO_CM); - - //----------------- - // Retrieve the subdetector - std::string DCH_name(m_DCH_name.value()); - if ( 0 == m_geoSvc->getDetector()->detectors().count(DCH_name) ) - { - ThrowException( "Detector <<" + DCH_name + ">> does not exist." ); - } - - dd4hep::DetElement DCH_DE = m_geoSvc->getDetector()->detectors().at(DCH_name); - - /////////////////////////////////////////////////////////////////////////////////// - /////////////////////////// retrieve data extension ////////////////////////// - /////////////////////////////////////////////////////////////////////////////////// - this->dch_data = DCH_DE.extension(); - - /////////////////////////////////////////////////////////////////////////////////// - - //----------------- - // Retrieve the readout associated with the detector element (subdetector) - dd4hep::SensitiveDetector dch_sd = m_geoSvc->getDetector()->sensitiveDetector(DCH_name); - if(not dch_sd.isValid() ) - ThrowException("No valid Sensitive Detector was found for detector <<" + DCH_name + ">>."); - - dd4hep::Readout dch_readout = dch_sd.readout(); - // set the cellID decoder - m_decoder = dch_readout.idSpec().decoder(); - - /////////////////////////////////////////////////////////////////////////////////// - ////////////////// initialize Walaa's code for Cluster counting ///////////////// - /////////////////////////////////////////////////////////////////////////////////// - debug() << Form("Opening %s...",m_fileDataAlg.value().c_str())<> not found."); - flData = new AlgData(); - flData->Load_file(m_fileDataAlg.value().c_str()); - flData->Load_interp(); - debug() << Form("Opening %s... Done",m_fileDataAlg.value().c_str())<SetDirectory(0); - hDww = new TH1D("hDww", "Distance hit projection to the wire, in cm. Should be zero", 100,0,1); - hDww->SetDirectory(0); - hSz = new TH1D("hSz", "Smearing along the wire, in cm", 100,0,5*m_z_resolution.value()); - hSz->SetDirectory(0); - hSxy = new TH1D("hSxy", "Smearing perpendicular the wire, in cm", 100,0,5*m_xy_resolution.value()); - hSxy->SetDirectory(0); - } - return StatusCode::SUCCESS; - } - - - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// operator() //////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -std::tuple -DCHdigi::operator()(const edm4hep::SimTrackerHitCollection& input_sim_hits, - const edm4hep::EventHeaderCollection& headers) const { - - // initialize seed for random engine - this->PrepareRandomEngine( headers ); - - debug() << "Input Sim Hit collection size: " << input_sim_hits.size() << endmsg; - - // Create the collections we are going to return - extension::DriftChamberDigiV2Collection output_digi_hits; - extension::MCRecoDriftChamberDigiV2AssociationCollection output_digi_sim_association; - - //loop over hit collection - for (const auto& input_sim_hit : input_sim_hits) - { - dd4hep::DDSegmentation::CellID cellid = input_sim_hit.getCellID(); - int ilayer = this->CalculateLayerFromCellID(cellid ); - int nphi = this->CalculateNphiFromCellID(cellid ); - auto hit_position = Convert_EDM4hepVector_to_TVector3( input_sim_hit.getPosition(), MM_TO_CM ); - - // ------------------------------------------------------------------------- - // calculate hit position projection into the wire - TVector3 hit_to_wire_vector = this->Calculate_hitpos_to_wire_vector(ilayer, nphi,hit_position); - TVector3 hit_projection_on_the_wire = hit_position + hit_to_wire_vector; - if( m_create_debug_histos.value() ) - { - double distance_hit_wire = hit_to_wire_vector.Mag(); - hDpw->Fill(distance_hit_wire); - } - TVector3 wire_direction_ez = this->Calculate_wire_vector_ez(ilayer, nphi); - - // ------------------------------------------------------------------------- - // smear the position - - // smear position along the wire - double smearing_z = m_gauss_z_cm( m_engine ); - if( m_create_debug_histos.value() ) hSz->Fill( smearing_z ); - - hit_projection_on_the_wire += smearing_z*(wire_direction_ez.Unit()); - if( m_create_debug_histos.value() ) - { - // the distance from the hit projection and the wire should be zero - TVector3 dummy_vector = this->Calculate_hitpos_to_wire_vector(ilayer, nphi,hit_projection_on_the_wire); - hDww->Fill( dummy_vector.Mag() ); - } - - // smear position perpendicular to the wire - double smearing_xy = m_gauss_xy_cm( m_engine ); - if( m_create_debug_histos.value() ) hSxy->Fill( smearing_xy ); - float distanceToWire_real = hit_to_wire_vector.Mag(); - - // protect against negative values - float distanceToWire_smeared = std::max(0.0, distanceToWire_real +smearing_xy ); - - std::int32_t type = 0; - std::int32_t quality = 0; - float eDepError =0; - // length units back to mm - auto positionSW = Convert_TVector3_to_EDM4hepVector(hit_projection_on_the_wire, 1./MM_TO_CM ); - auto directionSW = Convert_TVector3_to_EDM4hepVector(wire_direction_ez , 1./MM_TO_CM ); - float distanceToWire = distanceToWire_smeared/MM_TO_CM; - - auto [nCluster,nElectronPerCluster] = CalculateClusters(input_sim_hit); - - extension::MutableDriftChamberDigiV2 oDCHdigihit; - oDCHdigihit.setCellID(input_sim_hit.getCellID()); - oDCHdigihit.setType(type); - oDCHdigihit.setQuality(quality); - oDCHdigihit.setTime(input_sim_hit.getTime()); - oDCHdigihit.setEDep(input_sim_hit.getEDep()); - oDCHdigihit.setEDepError(eDepError); - oDCHdigihit.setPosition(positionSW); - oDCHdigihit.setDirectionSW(directionSW); - oDCHdigihit.setDistanceToWire(distanceToWire); - oDCHdigihit.setNCluster(nCluster); - oDCHdigihit.setNElectronPerCluster(nElectronPerCluster); - - output_digi_hits.push_back(oDCHdigihit); - - extension::MutableMCRecoDriftChamberDigiV2Association oDCHsimdigi_association; - oDCHsimdigi_association.setDigi( oDCHdigihit ); - oDCHsimdigi_association.setSim( input_sim_hit ); - output_digi_sim_association.push_back(oDCHsimdigi_association); - - }// end loop over hit collection - - ///////////////////////////////////////////////////////////////// - return std::make_tuple(std::move(output_digi_hits),std::move(output_digi_sim_association)); -} - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// finalize ////////////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -StatusCode DCHdigi::finalize() -{ - if( m_create_debug_histos.value() ) - { - std::unique_ptr ofile{TFile::Open ( m_out_debug_filename.value().c_str() , "recreate" ) }; - ofile->cd(); - hDpw->Write(); - hDww->Write(); - hSxy->Write(); - hSz->Write(); - } - - return StatusCode::SUCCESS; -} - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// ThrowException //////////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -void DCHdigi::ThrowException(std::string s) const { - error() << s.c_str() << endmsg; - throw std::runtime_error(s); - } - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// PrintConfiguration //////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -void DCHdigi::PrintConfiguration(std::ostream& io) -{ - io << "DCHdigi will use the following components:\n"; - io << "\tGeometry Service: " << m_geoSvcName.value().c_str() << "\n"; - io << "\tUID Service: " << m_uidSvcName.value().c_str() << "\n"; - io << "\tDetector name: " << m_DCH_name.value().c_str() << "\n"; - io << "\t\t|--Volume bitfield: " << m_decoder->fieldDescription().c_str() << "\n"; - io << "\t\t|--Number of layers: " << dch_data->database.size() << "\n"; - io << "\tCluster distributions taken from: " << m_fileDataAlg.value().c_str() << "\n"; - io << "\tResolution along the wire (mm): " << m_z_resolution.value() << "\n"; - io << "\tResolution perp. to the wire (mm): " << m_xy_resolution.value() << "\n"; - return; -} - -void DCHdigi::PrepareRandomEngine(const edm4hep::EventHeaderCollection& headers) const -{ - uint32_t evt_n = headers[0].getEventNumber(); - uint32_t run_n = headers[0].getRunNumber(); - size_t seed = m_uidSvc->getUniqueID(evt_n, run_n, this->name() ); - m_engine.seed(seed); - // advance internal state to minimize possibility of creating correlations - m_engine.discard(10); - for(int i=0; i<10; ++i) myRandom.Rndm(); -} - -/////////////////////////////////////////////////////////////////////////////////////// -///// Ancillary functions for calculating the distance to the wire //////// -/////////////////////////////////////////////////////////////////////////////////////// - - -TVector3 DCHdigi::Calculate_wire_vector_ez(int ilayer, int nphi) const -{ - auto & l = this->dch_data->database.at(ilayer); - - // See original paper Hoshina et al, Computer Physics Communications 153 (2003) 3 - // eq. 2.9, for the definition of ez, vector along the wire - - // initialize some variables - int stereosign = l.StereoSign(); - double rz0 = l.radius_sw_z0; - double dphi = dch_data->twist_angle; - // kappa is the same as in eq. 2.9 - double kappa = (1./dch_data->Lhalf)*tan(dphi/2); - - //--- calculating wire position - // the points p1 and p2 correspond to the ends of the wire - - // point 1 - // double x1 = rz0; // m - // double y1 = 0.; // m - // double z1 = 0.; // m - double x1 = rz0; // m - double y1 = -stereosign*rz0*kappa*dch_data->Lhalf; // m - double z1 = -dch_data->Lhalf; // m - - TVector3 p1 (x1,y1,z1); - - - // point 2 - double x2 = rz0; // m - double y2 = stereosign*rz0*kappa*dch_data->Lhalf; // m - double z2 = dch_data->Lhalf; // m - - TVector3 p2 (x2,y2,z2); - - // calculate phi rotation of whole twisted tube, ie, rotation at z=0 - double phi_z0 = Calculate_wire_phi_z0(ilayer,nphi); - p1.RotateZ(phi_z0); - p2.RotateZ(phi_z0); - - //--- end calculating wire position - - return (p2-p1).Unit(); - -} - -TVector3 DCHdigi::Calculate_wire_z0_point(int ilayer, int nphi) const -{ - auto & l = this->dch_data->database.at(ilayer); - double rz0 = l.radius_sw_z0; - TVector3 p1 (rz0,0,0); - double phi_z0 = Calculate_wire_phi_z0(ilayer,nphi); - p1.RotateZ(phi_z0); - return p1; -} - -// calculate phi rotation of whole twisted tube, ie, rotation at z=0 -double DCHdigi::Calculate_wire_phi_z0(int ilayer, int nphi) const -{ - auto & l = this->dch_data->database.at(ilayer); - int ncells = l.nwires/2; - double phistep = TMath::TwoPi()/ncells; - double phi_z0 = (nphi + 0.25*(l.layer%2))*phistep; - return phi_z0; -} - - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// Calculate vector from hit position to wire ///////////////// -/////////////////////////////////////////////////////////////////////////////////////// -TVector3 DCHdigi::Calculate_hitpos_to_wire_vector(int ilayer, int nphi, const TVector3 & hit_position /*in cm*/) const -{ - // Solution distance from a point to a line given here: - // https://en.wikipedia.org/wiki/Distance_from_a_point_to_a_line#Vector_formulation - TVector3 n = this->Calculate_wire_vector_ez(ilayer, nphi); - TVector3 a = this->Calculate_wire_z0_point (ilayer, nphi); - // Remember using cm as natural units of DD4hep consistently! - // TVector3 p {hit_position.x()*MM_TO_CM,hit_position.y()*MM_TO_CM,hit_position.z()*MM_TO_CM}; - - TVector3 a_minus_p = a - hit_position; - double a_minus_p_dot_n = a_minus_p.Dot( n ); - TVector3 scaled_n = a_minus_p_dot_n * n; - //hit_to_wire_vector = a_minus_p - scaled_n; - return (a_minus_p - scaled_n); -} - -/////////////////////////////////////////////////////////////////////////////////////// -/////////////////////// CalculateNClusters //////////////////////////////// -/////////////////////////////////////////////////////////////////////////////////////// -std::pair DCHdigi::CalculateClusters(const edm4hep::SimTrackerHit & input_sim_hit) const -{ - std::pair return_values = {0,0}; - uint32_t & Ncl = return_values.first; - uint32_t & ClSz = return_values.second; - //_________________SET NECESSARY PARAMETERS FOR THE CLS ALGORITHM-----WALAA_________________// - - auto thisParticle = input_sim_hit.getParticle(); - // Parameters needed for the CL Algo //Added by Walaa/////// - ///from Createclusters.cpp//// - // int fNClusters; - double me = 0.511; - // float rnd; - float Eloss= 0.0; - // float DeltaE_per_track= 0.0; - // float EtotCell= 0.0; - int NCl; - int NCl1,NClp; - // float ClSz,ClSzP; - // float Rt=0.87; - float EIzp=15.8; - float ExECl1; - float ExECl; - float cut= 1000;//controlla - // float prc=0.83; - float EIzs=25.6; - // float ratio=prc/EIzs; - // int NCltot=0; - // float Edep=0.0; - // float DeltaE=0.0; - // int count; - // bool loop; - // float prvDiff; - // float tmpDiff[10]; - // float vecExECl[10]; - // float mintmpDiff; - // int iloop; - float ExECl1totRec; - // float SgmaxExECl; - // int Ncl; - float rndCorr; - float maxExECl=0.0; - // float ExSgm=0.0; - float totExECl; -// int ParentID; - const int nhEp=10;//10 - float hEpcut[10]={100,200,300,400,500,600,700,800,900,1000}; - int minE=1000; - int maxE=10000; - int binE=1000; - int nhE=(maxE-minE)/binE; - TString parClass; - // int Ncltot=0; - int NEltot=0; - // float Etot=0.0; - float maxEx0len,ExSgmlen; - int choice= 2; - // int Nev=500000; - // int val; - // float EtCut=0.7; - ///////end from Createclusters.cpp//// - -///////from Createclusters.h///////// - // TF1Convolution **cnvl =new TF1Convolution*[100];//100 for 1 cm gas box //Added by Walaa - // TF1 **fcnvl=new TF1 *[100]; //Added by Walaa - // float fnorm[100]; //Added by Walaa - double vecpar[5]; //Added by Walaa - std::vector vecEtr,vecEtrTot; - float sumVecTot,meanVecTot,sumVec,meanVec; -/////////end from Createclusters.h///////// - - // get rid of spurious compiler error, Werror=unused-but-set-variable - (void)NCl; - (void)vecpar; - (void)sumVecTot; - - int NCld; - TString parName; - double bg, Momentum, mass=-1e300; - std::vector vecExtraD; - - double Ffrac, Fmpv1, Fsgm1, Fmpv2, Fsgm2; - std::vector CorrpMean,CorrpSgm,Corrdgmean,Corrdgsgm,Corrdglfrac,Corrdglmpvl,Corrdglsgml,Corrdglmeang,Corrdglsgmg; - float Ekdelta; - float maxEx0, maxExSlp, ExSgmlep, ExSgmhad; - float MPVEx, SgmEx, MeanEx1, SgmEx1, frac, Slp, CorrSlp, CorrInt; - double LengthTrack, Etot_per_track; - //////end parameters - bool IsSecondaryWithinDCH = false; - { - auto vertex = thisParticle.getVertex(); // in mm - auto vertexRsquared = vertex[0]*vertex[0] + vertex[1]*vertex[1]; - auto vertexZabs = std::fabs(vertex[2]); - float DCH_halflengh = 2000; - float DCH_rin_squared = 350*350; - float DCH_rout_squared = 2000*2000; - IsSecondaryWithinDCH = (vertexZabs < DCH_halflengh) && (vertexRsquared>DCH_rin_squared) && (vertexRsquaredget_Ffrac(bg);//0.9; - Fmpv1=flData->get_Fmpv1(bg);//43; - Fsgm1=flData->get_Fsgm1(bg);//9.64; - Fmpv2=flData->get_Fmpv2(bg);//311.3; - Fsgm2=flData->get_Fsgm2(bg);//35; - CorrpMean=flData->get_ClSzCorrpmean(bg); - CorrpSgm=flData->get_ClSzCorrpsgm(bg); - Corrdgmean=flData->get_ClSzCorrdgmean(bg); - Corrdgsgm=flData->get_ClSzCorrdgsgm(bg); - Corrdglfrac=flData->get_ClSzCorrdglfrac(bg); - Corrdglmpvl=flData->get_ClSzCorrdglmpvl(bg); - Corrdglsgml=flData->get_ClSzCorrdglsgml(bg); - Corrdglmeang=flData->get_ClSzCorrdglmeang(bg); - Corrdglsgmg=flData->get_ClSzCorrdglsgmg(bg); - maxEx0=flData->get_maxEx0(bg);//379.4 electron 100 gev - maxExSlp=flData->get_maxExSlp(); - // lepton PDG id goes from 11 to 16 (antiparticles have negative id) - bool IsLepton = (11 <= abs(pdgid) ) && (16 >= abs(pdgid) ); - if( IsLepton ){ - ExSgmlep=flData->get_ExSgmlep(); - }else{ - ExSgmhad=flData->get_ExSgmhad(); - } - - MPVEx=flData->get_MPVExtra(bg);//103.2 - SgmEx=flData->get_SgmExtra(bg);//28.5;// - MeanEx1=flData->get_MeanExtra1(bg);//13.8;// - SgmEx1=flData->get_SgmExtra1(bg);//9.72;// - frac=flData->get_FracExtra1(bg);//0.13;// - Slp=flData->get_SlopeExtra1(bg);//7.95;// - debug() << " maxEx0 "<< maxEx0<< " ExSgmhad "<get_ClSzCorrSlp(bg); - CorrInt=flData->get_ClSzCorrInt(bg); - - vecpar[0]=Ffrac; - vecpar[1]=Fmpv1; - vecpar[2]=Fsgm1; - vecpar[3]=Fmpv2; - vecpar[4]=Fsgm2; - - double Tmax = (2.0*me*pow(bg,2)/(1+(2.0*(1+pow(bg,2))*me/mass)+pow(me/mass,2)))*1e+6; - float maxEcut=cut; - if(TmaxSetRange(0,maxEcut); - * fcnvl[0]->SetParameters(vecpar[0],vecpar[1],vecpar[2],vecpar[3],vecpar[4]); - * cnvl[0]=nullptr;*/ - - - TF1 *land= new TF1("land","landaun"); - land->SetParameter(0,1); - land->SetParameter(1,MPVEx); - land->SetParameter(2,SgmEx); - land->SetRange(0,maxEcut); - - TF1 *exGauss=new TF1("exGauss","[0]*([1]*TMath::Gaus(x,[2],[3],true)+(1.0-[1])*TMath::Exp(-x/[4])/[4])"); - exGauss->SetParameter(0,1); - exGauss->SetParameter(1,frac); - exGauss->SetParameter(2,MeanEx1); - exGauss->SetParameter(3,SgmEx1); - exGauss->SetParameter(4,Slp); - exGauss->SetRange(0,90); - - - // int MaxEv = Entries > Nev ? Nev : Entries; - // cout<<" MAXEV"<myRandom.Gaus(0,ExSgmlen))/maxExSlp; - // cout<<"maxExECl= "<(EIzp+EIzs)&&maxExECl>totExECl && while1counter<1e6){ - while1counter++; - // cout<<"start again"<Eloss) ExECl=Eloss; - totExECl+=ExECl; - // cout<<"totExECl= " <myRandom.Uniform(0,1); - if(rndCorrmyRandom.Gaus(Corrdglmeang[tmphE],Corrdglsgmg[tmphE]); - }else{ - tmpCl=this->myRandom.Landau(Corrdglmpvl[tmphE],Corrdglsgml[tmphE]); - } - }else{ - tmpCl=this->myRandom.Gaus(Corrdgmean[tmphE],Corrdgsgm[tmphE]); - } - - ClSz=TMath::Nint(vecExtraD[i]*CorrSlp+CorrInt-tmpCl); - if(ClSz<2) {ClSz=2;} - // TODO Alvaro: check what is ClSz??? - // output_digi_hit.setClusterSize(ClSz); - // hClSzRec->Fill(ClSz); - // hClSzDRec->Fill(ClSz); - // fClusterCharge.push_back(ClSz); - NEltot+=ClSz; - - } - - Ncl=NCl1+NClp+NCld; - debug() <<"Ncl= "<< Ncl<<" NCl1= "<size() << endmsg; - - // output_digi_hit.setClusterCount(Ncl); - - /* for (int icl=0;icl