diff --git a/DCHdigi/include/AlgData.h b/DCHdigi/include/AlgData.h new file mode 100644 index 0000000..b1a7e99 --- /dev/null +++ b/DCHdigi/include/AlgData.h @@ -0,0 +1,589 @@ +#include + +#include "TTree.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 + +using namespace std; + + +class AlgData{ +public: + AlgData(); + ~AlgData(); + +public: + void read_file(TString dataAlg); + TF1* read_graph(TString dataAlg, TString cvName, TString fitName); + + + TF1* get_fit(); + TFormula* get_formula(); + double get_Fitvalue(double betagamma) const; + ROOT::Math::Interpolator* Interpvalue(vector bg, vector Ydata,int Intertype); + void Load_file(TString dataAlg); + void Load_interp(); + double get_MPVExtra (double betagamma) const; + double get_SgmExtra (double betagamma) const; + double get_MeanExtra1 (double betagamma) const; + double get_SgmExtra1 (double betagamma) const; + double get_SlopeExtra1 (double betagamma)const; + double get_FracExtra1(double betagamma) const; + double get_FfracExtra(double betagamma) ; + double get_maxEx0(double betagamma); + double get_maxExSlp(); + double get_ExSgmlep(); + double get_ExSgmhad(); + double get_Ffrac(double betagamma); + double get_Fmpv1(double betagamma); + double get_Fsgm1(double betagamma); + double get_Fmpv2(double betagamma); + double get_Fsgm2(double betagamma); + double get_ClSzCorrInt(double betagamma); + double get_ClSzCorrSlp(double betagamma); + vector get_ClSzCorrpmean(double betagamma); + vector get_ClSzCorrpsgm(double betagamma); + + vector get_ClSzCorrdgmean(double betagamma); + vector get_ClSzCorrdgsgm(double betagamma); + + vector get_ClSzCorrdglfrac(double betagamma); + vector get_ClSzCorrdglmpvl(double betagamma); + vector get_ClSzCorrdglsgml(double betagamma); + vector get_ClSzCorrdglmeang(double betagamma); + 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; + + vector *tmpCorrmeanSliceTot=nullptr; + vector *tmpErrCorrmeanSliceTot=nullptr; + vector *tmpCorrsgmSliceTot=nullptr; + vector *tmpErrCorrsgmSliceTot=nullptr; + + vector *tmpCorrdglfracSliceTot=nullptr; + vector *tmpErrCorrdglfracSliceTot=nullptr; + vector *tmpCorrdglmeangSliceTot=nullptr; + vector *tmpErrCorrdglmeangSliceTot=nullptr; + vector *tmpCorrdglsgmgSliceTot=nullptr; + vector *tmpErrCorrdglsgmgSliceTot=nullptr; + vector *tmpCorrdglmpvSliceTot=nullptr; + vector *tmpErrCorrdglmpvSliceTot=nullptr; + vector *tmpCorrdglsgmlSliceTot=nullptr; + vector *tmpErrCorrdglsgmlSliceTot=nullptr; + + vector *tmpCorrdgmeanSliceTot=nullptr; + vector *tmpErrCorrdgmeanSliceTot=nullptr; + vector *tmpCorrdgsgmSliceTot=nullptr; + vector *tmpErrCorrdgsgmSliceTot=nullptr; + + //____________________________________________________// + vector bgv,bglep,bghad; + vector maxEx0Tot,ErrmaxEx0Tot,maxExSlpTot,ErrmaxExSlpTot,ExSgmTotlep,ExSgmTothad; + vector tFfracTot, tFerrfracTot, tFmpv1Tot, tFerrmpv1Tot, tFsgm1Tot, tFerrsgm1Tot, tFmpv2Tot, tFerrmpv2Tot, tFsgm2Tot, tFerrsgm2Tot; + vector CorrmeanSliceTot[20],ErrCorrmeanSliceTot[20],CorrsgmSliceTot[20],ErrCorrsgmSliceTot[20]; + vector CorrIntTot,ErrCorrIntTot,CorrSlpTot,ErrCorrSlpTot; + vector CorrdglfracSliceTot[20],ErrCorrdglfracSliceTot[20],CorrdglmeangSliceTot[20],ErrCorrdglmeangSliceTot[20],CorrdglsgmgSliceTot[20],ErrCorrdglsgmgSliceTot[20],CorrdglmpvSliceTot[20],ErrCorrdglmpvSliceTot[20],CorrdglsgmlSliceTot[20],ErrCorrdglsgmlSliceTot[20]; + 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; + + vector Corrmean; + vector Corrsgm; + + vector Corrdgmean; + vector Corrdgsgm; + + vector Corrdglfrac; + vector Corrdglmeang; + vector Corrdglsgmg; + vector Corrdglmpv; + vector Corrdglsgml; + + vector ClSzCorrpmean; + vector ClSzCorrpsgm; + vector ClSzCorrdgmean; + vector ClSzCorrdgsgm; + vector ClSzCorrdglfrac; + vector ClSzCorrdglmpvl; + vector ClSzCorrdglsgml; + vector ClSzCorrdglmeang; + vector ClSzCorrdglsgmg; + + double maxExSlp; + double ExSgmlep; + double ExSgmhad; + void Calc_maxExSlp(); + void Calc_ExSgmlep(); + void Calc_ExSgmhad(); + + + + //___________________________________________// + + 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; + + 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); +// 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( 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)); +// cout<at(k)<GetEntries();i++) { + datalep->GetEntry(i); + bglep.push_back(bgTlep); + ExSgmTotlep.push_back(tmpExSgmTotlep); +// 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(vectorbg,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); + } + + vector AlgData::get_ClSzCorrpmean(double betagamma){ + ClSzCorrpmean.clear(); +// cout<<"size "<size()<size();++i){ + ClSzCorrpmean.push_back(itpm[Corrmean.at(i)]->Eval(betagamma)); + } + + return ClSzCorrpmean; + } + + 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; + } + + 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; + } + + 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; + } + + 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; + } + + 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; + } + + 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; + } + + 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; + } + + 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 m_input_sim_hits{"inputSimHits", Gaudi::DataHandle::Reader, this}; + // Output digitized tracker hit collection name + DataHandle m_output_digi_hits{"outputDigiHits", Gaudi::DataHandle::Writer, this}; + + // Detector readout name + Gaudi::Property m_readoutName{this, "readoutName", "CDCHHits", "Name of the detector readout"}; + // Pointer to the geometry service + ServiceHandle m_geoSvc; + // Decoder for the cellID + dd4hep::DDSegmentation::BitFieldCoder* m_decoder; + // Volume manager to get the physical cell sensitive volume + dd4hep::VolumeManager m_volman; + + // z position resolution in mm + FloatProperty m_z_resolution{this, "zResolution", 1.0, + "Spatial resolution in the z direction (from reading out the wires at both sides) [mm]"}; + // xy resolution in mm + FloatProperty m_xy_resolution{this, "xyResolution", 0.1, "Spatial resolution in the xy direction [mm]"}; + + // Random Number Service + IRndmGenSvc* m_randSvc; + // Gaussian random number generator used for the smearing of the z position + Rndm::Numbers m_gauss_z; + // Gaussian random number generator used for the smearing of the xy position + Rndm::Numbers m_gauss_xy; + + // Parameters needed for the CL Algo //Added by Walaa/////// + ///from Createclusters.cpp//// + Int_t fNClusters; + double me = 0.511; + float rnd; + float Eloss= 0.0; float DeltaE_per_track= 0.0; float EtotCell= 0.0; + int NCl,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_t 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///////// + int NCld; + TString parName; + double bg, Momentum, mass; + AlgData *flData; + 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_t LengthTrack, Etot_per_track; + //////end parameters + +}; diff --git a/DCHdigi/src/DCHsimpleDigitizerWClustercount.cpp b/DCHdigi/src/DCHsimpleDigitizerWClustercount.cpp new file mode 100644 index 0000000..72be0f2 --- /dev/null +++ b/DCHdigi/src/DCHsimpleDigitizerWClustercount.cpp @@ -0,0 +1,344 @@ +#include "DCHsimpleDigitizerWClustercount.h" + +// DD4hep +#include "DD4hep/Detector.h" +#include "DDRec/Vector3D.h" + +// ROOT +#include "Math/Cylindrical3D.h" +#include "AlgData.h" +#include "TParticlePDG.h" +#include "TDatabasePDG.h" +#include "TMath.h" +#include "edm4hep/ParticleIDData.h" +#include "TRandom.h" +DECLARE_COMPONENT(DCHsimpleDigitizerWClustercount) + +DCHsimpleDigitizerWClustercount::DCHsimpleDigitizerWClustercount(const std::string& aName, ISvcLocator* aSvcLoc) + : GaudiAlgorithm(aName, aSvcLoc), m_geoSvc("GeoSvc", "DCHsimpleDigitizerWClustercount") { + declareProperty("inputSimHits", m_input_sim_hits, "Input sim tracker hit collection name"); + declareProperty("outputDigiHits", m_output_digi_hits, "Output digitized tracker hit collection name"); +} + +DCHsimpleDigitizerWClustercount::~DCHsimpleDigitizerWClustercount() {} + +StatusCode DCHsimpleDigitizerWClustercount::initialize() { + // Initialize random services + if (service("RndmGenSvc", m_randSvc).isFailure()) { + error() << "Couldn't get RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } + if (m_gauss_z.initialize(m_randSvc, Rndm::Gauss(0., m_z_resolution)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } + if (m_gauss_xy.initialize(m_randSvc, Rndm::Gauss(0., m_xy_resolution)).isFailure()) { + error() << "Couldn't initialize RndmGenSvc!" << endmsg; + return StatusCode::FAILURE; + } + + // check if readout exists + if (m_geoSvc->lcdd()->readouts().find(m_readoutName) == m_geoSvc->lcdd()->readouts().end()) { + error() << "Readout <<" << m_readoutName << ">> does not exist." << endmsg; + return StatusCode::FAILURE; + } + // set the cellID decoder + m_decoder = m_geoSvc->lcdd()->readout(m_readoutName).idSpec().decoder(); + // retrieve the volume manager + m_volman = m_geoSvc->lcdd()->volumeManager(); + +////Added by Walaa for the CLS/////////// + TString DataAlg= "/afs/cern.ch/work/w/welmeten/public/IDEA_Key4HEP_Jun22/Cluster_sim_code/DataAlgFORGEANT.root"; + debug() << "root file opened"<Load_file(DataAlg.Data()); + flData->Load_interp(); +//////////////////////////////////////// + return StatusCode::SUCCESS; +} + +StatusCode DCHsimpleDigitizerWClustercount::execute() { + // Get the input collection with Geant4 hits + const edm4hep::SimTrackerHitCollection* input_sim_hits = m_input_sim_hits.get(); + debug() << "Input Sim Hit collection size: " << input_sim_hits->size() << endmsg; + + // Digitize the sim hits + extension::DriftChamberDigiCollection* output_digi_hits = m_output_digi_hits.createAndPut(); + for (const auto& input_sim_hit : *input_sim_hits) { + auto output_digi_hit = output_digi_hits->create(); + // smear the hit position: need to go in the wire local frame to smear in the direction aligned/perpendicular with the wire for z/distanceToWire, taking e.g. stereo angle into account + // retrieve the cell detElement + dd4hep::DDSegmentation::CellID cellID = input_sim_hit.getCellID(); + auto cellDetElement = m_volman.lookupDetElement(cellID); + // retrieve the wire (in DD4hep 1.23 there is no easy way to access the volume daughters we have to pass by detElements, in later versions volumes can be used) + const std::string& wireDetElementName = + Form("superLayer_%d_layer_%d_phi_%d_wire", m_decoder->get(cellID, "superLayer"), + m_decoder->get(cellID, "layer"), m_decoder->get(cellID, "phi")); + dd4hep::DetElement wireDetElement = cellDetElement.child(wireDetElementName); + // get the transformation matrix used to place the wire + const auto& wireTransformMatrix = wireDetElement.nominal().worldTransformation(); + // Retrieve global position in mm and apply unit transformation (translation matrix is stored in cm) + double simHitGlobalPosition[3] = {input_sim_hit.getPosition().x * dd4hep::mm, + input_sim_hit.getPosition().y * dd4hep::mm, + input_sim_hit.getPosition().z * dd4hep::mm}; + double simHitLocalPosition[3] = {0, 0, 0}; + // get the simHit coordinate in cm in the wire reference frame to be able to apply smearing of radius perpendicular to the wire + wireTransformMatrix.MasterToLocal(simHitGlobalPosition, simHitLocalPosition); + debug() << "Cell ID string: " << m_decoder->valueString(cellID) << endmsg; + debug() << "Global simHit x " << simHitGlobalPosition[0] << " --> Local simHit x " << simHitLocalPosition[0] + << " in cm" << endmsg; + debug() << "Global simHit y " << simHitGlobalPosition[1] << " --> Local simHit y " << simHitLocalPosition[1] + << " in cm" << endmsg; + debug() << "Global simHit z " << simHitGlobalPosition[2] << " --> Local simHit z " << simHitLocalPosition[2] + << " in cm" << endmsg; + // build a vector to easily apply smearing of distance to the wire + dd4hep::rec::Vector3D simHitLocalPositionVector(simHitLocalPosition[0], simHitLocalPosition[1], + simHitLocalPosition[2]); + // get the smeared distance to the wire (cylindrical coordinate as the smearing should be perpendicular to the wire) + debug() << "Original distance to wire: " << simHitLocalPositionVector.rho() << endmsg; + double smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; + while(smearedDistanceToWire < 0){ + debug() << "Negative smearedDistanceToWire (" << smearedDistanceToWire << ") shooting another random number" << endmsg; + smearedDistanceToWire = simHitLocalPositionVector.rho() + m_gauss_xy.shoot() * dd4hep::mm; + } + debug() << "Smeared distance to wire: " << smearedDistanceToWire << endmsg; + // smear the z position (in local coordinate the z axis is aligned with the wire i.e. it take the stereo angle into account); + double smearedZ = simHitLocalPositionVector.z() + m_gauss_z.shoot() * dd4hep::mm; + // fill the output DriftChamberDigi (making sure we are back in mm) + output_digi_hit.setCellID(cellID); + output_digi_hit.setDistanceToWire(smearedDistanceToWire / dd4hep::mm); + output_digi_hit.setZPositionAlongWire(smearedZ / dd4hep::mm); + + //_________________SET NECESSARY PARAMETERS FOR THE CLS ALGORITHM-----WALAA_________________// + + Momentum = sqrt((input_sim_hit.getMomentum().x * input_sim_hit.getMomentum().x) + (input_sim_hit.getMomentum().y * input_sim_hit.getMomentum().y) + (input_sim_hit.getMomentum().z * input_sim_hit.getMomentum().z)) ; + + NCld =0; + int ParentID = input_sim_hit.getCellID(); + int pdgid = static_cast(input_sim_hit.getTime()); + + if(ParentID==1&&pdgid==11){ + NCld++; + Ekdelta=(TMath::Sqrt(Momentum*Momentum+me*me)-me)*1e6; + vecExtraD.push_back(Ekdelta); + } + debug() << "NCld= "<0) debug() << "There is delta rays"<GetParticle(pdgid); + if(partRoot!=nullptr){ + parName=partRoot->GetName(); + parClass=partRoot->ParticleClass(); + debug() << "parName "<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"<Gaus(0,ExSgmlen))/maxExSlp; + // cout<<"maxExECl= "<(EIzp+EIzs)&&maxExECl>totExECl){ + // cout<<"start again"<Eloss) ExECl=Eloss; + totExECl+=ExECl; + // cout<<"totExECl= " <Uniform(0,1); + if(rndCorrGaus(Corrdglmeang[tmphE],Corrdglsgmg[tmphE]); + }else{ + tmpCl=gRandom->Landau(Corrdglmpvl[tmphE],Corrdglsgml[tmphE]); + } + }else{ + tmpCl=gRandom->Gaus(Corrdgmean[tmphE],Corrdgsgm[tmphE]); + } + + ClSz=TMath::Nint(vecExtraD[i]*CorrSlp+CorrInt-tmpCl); + if(ClSz<2) {ClSz=2;} + 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