21 #include "FairLogger.h"
22 #include "FairMCPoint.h"
23 #include "FairTrackParam.h"
33 #include "KFParticleTopoReconstructor.h"
38 #include "TClonesArray.h"
44 #include <boost/assign/list_of.hpp>
63 #define M2E 2.6112004954086e-7
66 using boost::assign::list_of;
69 : FairTask(
"CbmAnaConversion")
82 , fhElectronSources(NULL)
83 , fhNofPi0_perEvent(NULL)
84 , fhNofPi0_perEvent_cut(NULL)
85 , fhNofPi0_perEvent_cut2(NULL)
86 , fhNofPi0_perEvent_cut3(NULL)
87 , fhNofEta_perEvent(NULL)
88 , fhNofEta_perEvent_cut(NULL)
89 , fhNofEta_perEvent_cut2(NULL)
93 , fhPi0_pt_vs_rap(NULL)
95 , fhPi0_theta_vs_rap(NULL)
97 , fhEta_pt_vs_rap(NULL)
99 , fhEta_theta_vs_rap(NULL)
101 , fhRho_pt_vs_rap(NULL)
103 , fhRho_theta_vs_rap(NULL)
105 , fhOmega_pt_vs_rap(NULL)
106 , fhOmega_theta(NULL)
107 , fhOmega_theta_vs_rap(NULL)
108 , fhElectronsFromPi0_z(NULL)
109 , fhNofTracks_mctrack(NULL)
110 , fhNofTracks_globaltrack(NULL)
111 , fhInvariantMass_test(NULL)
112 , fhInvariantMass_test2(NULL)
113 , fhInvariantMass_test3(NULL)
114 , fhInvariantMassReco_test(NULL)
115 , fhInvariantMassReco_test2(NULL)
116 , fhInvariantMassReco_test3(NULL)
117 , fhInvariantMassReco_pi0(NULL)
118 , fhMomentum_MCvsReco(NULL)
119 , fhMomentum_MCvsReco_diff(NULL)
120 , fhSearchGammas(NULL)
125 , fRichRingMatches(NULL)
128 , fStsTrackMatches(NULL)
129 , fGlobalTracks(NULL)
130 , fhANN_output_electrons(NULL)
131 , fhANN_output_electrons2(NULL)
132 , fhANN_output_electrons_chiCut(NULL)
133 , fhANN_output_electrons_vs_p(NULL)
134 , fhANN_output_else(NULL)
135 , fhANN_output_else2(NULL)
136 , fhANN_output_else_chiCut(NULL)
137 , fhANN_output_else_vs_p(NULL)
143 , fKFparticleFinderQA(NULL)
147 , particlecounter_2daughters(0)
148 , particlecounter_3daughters(0)
149 , particlecounter_4daughters(0)
150 , particlecounter_all(0)
151 , fNofGeneratedPi0_allEvents(0)
152 , fNofPi0_kfparticle_allEvents(0)
153 , fNofGeneratedPi0(0)
154 , fNofPi0_kfparticle(0)
157 , fhPi0_NDaughters(NULL)
160 , fHistoList_tomography()
162 , fHistoList_reco_mom()
163 , fHistoList_kfparticle()
164 , fHistoList_richrings()
165 , fHistoList_furtherAnalyses()
169 , fRecoTracklistEPEM()
170 , fRecoTracklistEPEM_id()
171 , fRecoTracklistEPEM_chi()
172 , fRecoTracklistEPEM_gtid()
174 , fTestTracklist_momentum()
175 , fTestTracklist_chi()
176 , fTestTracklist_richInd()
177 , fTestTracklist_ndf()
178 , fTestTracklist_nofhits()
179 , fTestTracklist_noRichInd()
180 , fTestTracklist_noRichInd_MCindex()
181 , fTestTracklist_noRichInd_momentum()
182 , fTestTracklist_noRichInd_chi()
183 , fTestTracklist_noRichInd_richInd()
184 , fTestTracklist_noRichInd_gTrackId()
185 , fTestTracklist_noRichInd_ndf()
186 , fTestTracklist_noRichInd_nofhits()
188 , fRecoRefittedMomentum()
189 , fhNofElectrons_4epem(NULL)
190 , fhPi0_MC_occurence(NULL)
191 , fhPi0_MC_occurence2(NULL)
192 , fhPi0_Reco_occurence(NULL)
193 , fhPi0_Reco_occurence2(NULL)
194 , fhPi0_Reco_angles(NULL)
195 , fhPi0_Reco_chi(NULL)
196 , fhPi0_Reco_ndf(NULL)
197 , fhPi0_Reco_ndf_vs_chi(NULL)
198 , fhPi0_Reco_ndf_vs_startvertex(NULL)
199 , fhPi0_Reco_startvertex_vs_chi(NULL)
200 , fhPi0_Reco_startvertex_vs_nofhits(NULL)
201 , fhPi0_noRichInd_diffPhi(NULL)
202 , fhPi0_noRichInd_diffTheta(NULL)
203 , fhPi0_Reco_invmass_cases(NULL)
204 , fhPi0_Reco_noRichInd_invmass_cases(NULL)
205 , fhPi0_Reco_invmass(NULL)
206 , fhPi0_Reco_invmass_mc(NULL)
207 , fhPi0_Reco_noRichInd_invmass(NULL)
208 , fhPi0_Reco_noRichInd_invmass_mc(NULL)
209 , fhPi0_Reco_noRichInd_ndf_vs_nofhits(NULL)
210 , fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi(NULL)
211 , fhPi0_Reco_ndf_vs_nofhits(NULL)
212 , fhPi0_Reco_ndf_vs_nofhits_withChi(NULL)
213 , fhPi0_Reco_noRichInd_chi_vs_momentum(NULL)
214 , fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0(NULL)
215 , fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target(NULL)
216 , fhPi0_Reco_noRichInd_chi_vs_momentum_eRest(NULL)
217 , fhPi0_Reco_noRichInd_chi_vs_pt(NULL)
218 , fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0(NULL)
219 , fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target(NULL)
220 , fhPi0_Reco_noRichInd_chi_vs_pt_eRest(NULL)
221 , fhPi0_Reco_chi_vs_momentum(NULL)
222 , fhPi0_Reco_chi_vs_momentum_eFromPi0(NULL)
223 , fhPi0_Reco_chi_vs_momentum_eFromPi0_Target(NULL)
224 , fhPi0_Reco_chi_vs_pt(NULL)
225 , fhPi0_Reco_chi_vs_pt_eFromPi0(NULL)
226 , fhPi0_Reco_chi_vs_pt_eFromPi0_Target(NULL)
235 , fAnaTomography(NULL)
251 cout <<
"CbmAnaConversion::Init" << endl;
252 FairRootManager* ioman = FairRootManager::Instance();
254 Fatal(
"CbmAnaConversion::Init",
"RootManager not instantised!");
257 fRichPoints = (TClonesArray*) ioman->GetObject(
"RichPoint");
259 Fatal(
"CbmAnaConversion::Init",
"No RichPoint array!");
262 fMcTracks = (TClonesArray*) ioman->GetObject(
"MCTrack");
264 Fatal(
"CbmAnaConversion::Init",
"No MCTrack array!");
267 fStsTracks = (TClonesArray*) ioman->GetObject(
"StsTrack");
269 Fatal(
"CbmAnaConversion::Init",
"No StsTrack array!");
274 Fatal(
"CbmAnaConversion::Init",
"No StsTrackMatch array!");
277 fGlobalTracks = (TClonesArray*) ioman->GetObject(
"GlobalTrack");
279 Fatal(
"CbmAnaConversion::Init",
"No GlobalTrack array!");
289 if (
nullptr ==
fPrimVertex) { LOG(fatal) <<
"No PrimaryVertex array!"; }
291 fRichRings = (TClonesArray*) ioman->GetObject(
"RichRing");
293 Fatal(
"CbmAnaConversion::Init",
"No RichRing array!");
298 Fatal(
"CbmAnaConversion::Init",
"No RichRingMatch array!");
369 fhPdgCodes =
new TH1D(
"fhPdfCodes",
"fhPdgCodes;pdg code;#", 1000, 0, 1000);
371 new TH1D(
"fhNofElPrim",
"fhNofElPrim;Nof prim El;Entries", 10., -0.5, 9.5);
373 new TH1D(
"fhNofElSec",
"fhNofElSec;Nof Sec El;Entries", 20., -0.5, 19.5);
375 new TH1D(
"fhNofElAll",
"fhNofElAll;Nof All El;Entries", 30., -0.5, 29.5);
377 "fhNofPi0_perEvent;Nof pi0;Entries",
382 new TH1D(
"fhNofPi0_perEvent_cut",
383 "fhNofPi0_perEvent_cut (Z<10cm);Nof pi0;Entries",
388 new TH1D(
"fhNofPi0_perEvent_cut2",
389 "fhNofPi0_perEvent_cut2 (motherId = -1);Nof pi0;Entries",
394 new TH1D(
"fhNofPi0_perEvent_cut3",
395 "fhNofPi0_perEvent_cut3 (conversion before 70cm);Nof pi0;Entries",
400 "fhNofEta_perEvent",
"fhNofEta_perEvent;Nof eta;Entries", 100., -0.5, 99.5);
402 new TH1D(
"fhNofEta_perEvent_cut",
403 "fhNofEta_perEvent_cut (Z<4cm);Nof eta;Entries",
408 new TH1D(
"fhNofEta_perEvent_cut2",
409 "fhNofEta_perEvent_cut2 (motherId = -1);Nof eta;Entries",
413 fhPi0_z =
new TH1D(
"fhPi0_z",
"fhPi0_z;z in cm;Entries", 600., -0.5, 599.5);
415 new TH1D(
"fhPi0_z_cut",
"fhPi0_z_cut;z in cm;Entries", 600., -0.5, 599.5);
417 new TH1D(
"fhPi0_pt",
"fhPi0_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
419 "fhPi0_pt_vs_rap;p_{t} in GeV/c; rapidity y",
427 "fhPi0_theta;emission angle #theta in deg;Entries",
432 "fhPi0_theta_vs_rap;theta in deg;rapidity y",
440 new TH1D(
"fhEta_pt",
"fhEta_pt;p_{t} in GeV;Entries", 200., 0., 10.);
442 "fhEta_pt_vs_rap;p_{t} in GeV; rapidity y",
450 "fhEta_theta;emission angle #theta in deg;Entries",
455 new TH2D(
"fhEta_theta_vs_rap",
456 "fhEta_theta_vs_rap;emission angle #theta in deg;rapidity y",
464 new TH1D(
"fhRho_pt",
"fhRho_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
466 "fhRho_pt_vs_rap;p_{t} in GeV/c; rapidity y",
474 "fhRho_theta;emission angle #theta in deg;Entries",
479 new TH2D(
"fhRho_theta_vs_rap",
480 "fhRho_theta_vs_rap;emission angle #theta in deg;rapidity y",
488 new TH1D(
"fhOmega_pt",
"fhOmega_pt;p_{t} in GeV/c;Entries", 200., 0., 10.);
490 "fhOmega_pt_vs_rap;p_{t} in GeV; rapidity y",
498 "fhOmega_theta;emission angle #theta in deg;Entries",
503 new TH2D(
"fhOmega_theta_vs_rap",
504 "fhOmega_theta_vs_rap;emission angle #theta in deg;rapidity y",
514 "fhElectronSources",
"fhElectronSources;Source;Entries", 6., 0., 6.);
516 new TH1D(
"fhElectronsFromPi0_z",
517 "fhElectronsFromPi0_z (= pos. of gamma conversion);z [cm];Entries",
559 "fhNofTracks_mctrack",
"fhNofTracks_mctrack;nof;#", 1000., 0., 1000.);
562 "fhNofTracks_globaltrack;nof;#",
571 new TH1D(
"fhInvariant",
"fhInvariant;mass in GeV/c^{2}];#", 2000, 0., 2.);
573 new TH1D(
"fhInvariant2",
"fhInvariant2;mass in GeV/c^{2}];#", 2000, 0., 2.);
575 new TH1D(
"fhInvariant3",
"fhInvariant3;mass in GeV/c^{2}];#", 2000, 0., 2.);
581 "fhInvariantReco",
"fhInvariantReco;mass [GeV/c^2];#", 2000, 0., 2.);
583 "fhInvariantReco2",
"fhInvariantReco2;mass [GeV/c^2];#", 2000, 0., 2.);
585 "fhInvariantReco3",
"fhInvariantReco3;mass [GeV/c^2];#", 2000, 0., 2.);
591 "fhInvariantReco_pi0;mass [GeV/c^2];#",
600 "fhMomentum_MCvsReco;MC;Reco",
608 "fhMomentum_MCvsReco_diff;(MC-Reco)/MC",
617 new TH1D(
"fhSearchGammas",
"fhSearchGammas;mass;#", 100, -0.005, 0.995);
622 "fhANN_output_electrons",
"fhANN_output_electrons;ann output", 400, -2, 2.);
624 "fhANN_output_electrons2;ann output",
629 new TH1D(
"fhANN_output_electrons_chiCut",
630 "fhANN_output_electrons_chiCut;ann output",
635 new TH2D(
"fhANN_output_electrons_vs_p",
636 "fhANN_output_electrons_vs_p;momentum in GeV/c; ann output",
644 new TH1D(
"fhANN_output_else",
"fhANN_output_else;ann output", 400, -2, 2.);
646 "fhANN_output_else2",
"fhANN_output_else2;ann output", 400, -2, 2.);
648 "fhANN_output_else_chiCut;ann output",
653 new TH2D(
"fhANN_output_else_vs_p",
654 "fhANN_output_else_vs_p;momentum in GeV/c; ann output",
674 "fhPi0_NDaughters",
"fhPi0_NDaughters;number of daughers;#", 4, 0.5, 4.5);
675 fhPi0Ratio =
new TH1D(
"fhPi0Ratio",
"fhPi0Ratio; ratio;#", 1000, 0., 0.02);
676 fhPi0_mass =
new TH1D(
"fhPi0_mass",
"fhPi0_mass;mass;#", 500, 0., 0.5);
683 new TH1D(
"fhNofElectrons_4epem",
684 "fhNofElectrons_4epem;number of electrons per event;#",
691 new TH1D(
"fhPi0_MC_occurence",
"fhPi0_MC_occurence;;#", 20, 0, 20);
700 "!= -1: all pi0 not from target");
714 new TH1D(
"fhPi0_MC_occurence2",
"fhPi0_MC_occurence2;;#", 20, 0, 20);
717 1,
"!= -1: all pi0 not from target");
725 new TH1D(
"fhPi0_Reco_occurence",
"fhPi0_Reco_occurence;;#", 16, 0, 16);
745 new TH1D(
"fhPi0_Reco_occurence2",
"fhPi0_Reco_occurence2;;#", 16, 0, 16);
766 new TH1D(
"fhPi0_Reco_angles",
"fhPi0_Reco_angles;angle;#", 500, 0, 50);
769 new TH1D(
"fhPi0_Reco_chi",
"fhPi0_Reco_chi;chi;#", 500, 0, 500);
772 new TH1D(
"fhPi0_Reco_ndf",
"fhPi0_Reco_ndf;ndf;#", 500, 0, 50);
775 "fhPi0_Reco_ndf_vs_chi;ndf;chi",
784 new TH2D(
"fhPi0_Reco_ndf_vs_startvertex",
785 "fhPi0_Reco_ndf_vs_startvertex;ndf;startvertex",
794 new TH2D(
"fhPi0_Reco_startvertex_vs_chi",
795 "fhPi0_Reco_startvertex_vs_chi;startvertex;chi",
804 new TH2D(
"fhPi0_Reco_startvertex_vs_nofhits",
805 "fhPi0_Reco_startvertex_vs_nofhits;startvertex;nofhits",
814 "fhPi0_noRichInd_diffPhi;phi difference;#",
820 new TH1D(
"fhPi0_noRichInd_diffTheta",
821 "fhPi0_noRichInd_diffTheta;theta difference;#",
828 new TH2D(
"fhPi0_Reco_invmass_cases",
829 "fhPi0_Reco_invmass_cases;cases;invmass [GeV]",
844 new TH2D(
"fhPi0_Reco_noRichInd_invmass_cases",
845 "fhPi0_Reco_noRichInd_invmass_cases;cases;invmass in GeV^{2}",
860 7,
"richInd-no cuts");
870 "fhPi0_Reco_invmass",
"fhPi0_Reco_invmass;invmass [GeV];#", 300, 0, 3);
873 "fhPi0_Reco_invmass_mc;invmass_mc [GeV];#",
879 new TH1D(
"fhPi0_Reco_noRichInd_invmass",
880 "fhPi0_Reco_noRichInd_invmass;invmass_noRichInd [GeV];#",
886 new TH1D(
"fhPi0_Reco_noRichInd_invmass_mc",
887 "fhPi0_Reco_noRichInd_invmass_mc;invmass_noRichInd_mc [GeV];#",
894 new TH2D(
"fhPi0_Reco_noRichInd_ndf_vs_nofhits",
895 "fhPi0_Reco_noRichInd_ndf_vs_nofhits;ndf;nofhits",
904 "fhPi0_Reco_ndf_vs_nofhits;ndf;nofhits",
913 new TH2D(
"fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi",
914 "fhPi0_Reco_noRichInd_ndf_vs_nofhits_withChi;ndf;nofhits",
924 new TH2D(
"fhPi0_Reco_ndf_vs_nofhits_withChi",
925 "fhPi0_Reco_ndf_vs_nofhits_withChi;ndf;nofhits",
938 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_momentum",
939 "fhPi0_Reco_noRichInd_chi_vs_momentum;momentum [GeV];#chi^{2} of "
949 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0",
950 "fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0;momentum "
951 "[GeV];#chi^{2} of momentum fit",
961 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target",
962 "fhPi0_Reco_noRichInd_chi_vs_momentum_eFromPi0_Target;momentum "
963 "[GeV];#chi^{2} of momentum fit",
973 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_momentum_eRest",
974 "fhPi0_Reco_noRichInd_chi_vs_momentum_eRest;momentum "
975 "[GeV];#chi^{2} of momentum fit",
986 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_pt",
987 "fhPi0_Reco_noRichInd_chi_vs_pt;pt [GeV];#chi^{2} of momentum fit",
996 "fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0",
997 "fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0;pt [GeV];#chi^{2} of momentum fit",
1006 new TH2D(
"fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target",
1007 "fhPi0_Reco_noRichInd_chi_vs_pt_eFromPi0_Target;pt [GeV];#chi^{2} "
1018 "fhPi0_Reco_noRichInd_chi_vs_pt_eRest",
1019 "fhPi0_Reco_noRichInd_chi_vs_pt_eRest;pt [GeV];#chi^{2} of momentum fit",
1032 "fhPi0_Reco_chi_vs_momentum",
1033 "fhPi0_Reco_chi_vs_momentum;momentum [GeV];#chi^{2} of momentum fit",
1042 new TH2D(
"fhPi0_Reco_chi_vs_momentum_eFromPi0",
1043 "fhPi0_Reco_chi_vs_momentum_eFromPi0;momentum [GeV];#chi^{2} of "
1053 new TH2D(
"fhPi0_Reco_chi_vs_momentum_eFromPi0_Target",
1054 "fhPi0_Reco_chi_vs_momentum_eFromPi0_Target;momentum "
1055 "[GeV];#chi^{2} of momentum fit",
1066 new TH2D(
"fhPi0_Reco_chi_vs_pt",
1067 "fhPi0_Reco_chi_vs_pt;pt [GeV];#chi^{2} of momentum fit",
1076 new TH2D(
"fhPi0_Reco_chi_vs_pt_eFromPi0",
1077 "fhPi0_Reco_chi_vs_pt_eFromPi0;pt [GeV];#chi^{2} of momentum fit",
1086 "fhPi0_Reco_chi_vs_pt_eFromPi0_Target",
1087 "fhPi0_Reco_chi_vs_pt_eFromPi0_Target;pt [GeV];#chi^{2} of momentum fit",
1104 <<
"======================================================================="
1106 cout <<
"========== CbmAnaConversion, event No. " <<
fEventNum << endl;
1108 <<
"======================================================================="
1145 int countPrimEl = 0;
1150 int countPrimPart = 0;
1152 int countPi0MC_cut = 0;
1153 int countPi0MC_fromPrimary = 0;
1154 int countPi0MC_reconstructible = 0;
1156 int countEtaMC_cut = 0;
1157 int countEtaMC_fromPrimary = 0;
1162 Fatal(
"CbmAnaConversion::Exec",
"No PrimaryVertex array!");
1192 Int_t nofMcTracks =
fMcTracks->GetEntriesFast();
1194 for (
int i = 0;
i < nofMcTracks;
i++) {
1196 if (mctrack == NULL)
continue;
1202 if (mctrack->
GetMotherId() == -1) { countPrimPart++; }
1204 && TMath::Abs(mctrack->
GetPdgCode()) == 11) {
1208 && TMath::Abs(mctrack->
GetPdgCode()) == 11) {
1211 if (TMath::Abs(mctrack->
GetPdgCode()) == 11) { countAllEl++; }
1218 if (
v.Z() <= 4) { countPi0MC_cut++; }
1230 if (motherId == -1) {
1231 countPi0MC_fromPrimary++;
1241 if (reconstructible) countPi0MC_reconstructible++;
1245 if (TMath::Abs(mctrack->
GetPdgCode()) == 221) {
1247 TVector3
v, momentum;
1251 if (
v.Z() <= 4) { countEtaMC_cut++; }
1253 if (motherId == -1) {
1254 countEtaMC_fromPrimary++;
1267 TVector3
v, momentum;
1284 TVector3
v, momentum;
1297 if (TMath::Abs(mctrack->
GetPdgCode()) == 11) {
1303 cout <<
"CbmAnaConversion::Exec - Number of pi0 in MC sample: " << countPi0MC
1305 cout <<
"CbmAnaConversion::Exec - Number of pi0 from primary: "
1306 << countPi0MC_fromPrimary << endl;
1345 Int_t nofElectrons4epem = 0;
1349 for (Int_t iGTrack = 0; iGTrack < ngTracks; iGTrack++) {
1351 if (NULL == gTrack)
continue;
1357 if (stsInd < 0)
continue;
1359 if (stsTrack == NULL)
continue;
1362 if (stsMatch == NULL)
continue;
1365 if (stsMcTrackId < 0)
continue;
1367 if (mcTrack1 == NULL)
continue;
1371 TVector3 refittedMomentum;
1373 vector<CbmStsTrack> stsTracks;
1374 stsTracks.resize(1);
1375 stsTracks[0] = *stsTrack;
1376 vector<L1FieldRegion> vField;
1377 vector<float> chiPrim;
1381 const FairTrackParam* vtxTrack = stsTracks[0].GetParamFirst();
1382 vtxTrack->Momentum(refittedMomentum);
1383 float result_chi = chiPrim[0];
1388 vector<CbmStsTrack> stsTracks_electron;
1389 stsTracks_electron.resize(1);
1390 stsTracks_electron[0] = *stsTrack;
1391 vector<L1FieldRegion> vField_electron;
1392 vector<float> chiPrim_electron;
1393 vector<int> pidHypo_electron;
1394 pidHypo_electron.push_back(11);
1395 fPFFitter_electron.
Fit(stsTracks_electron, pidHypo_electron);
1397 stsTracks_electron, vField_electron, chiPrim_electron,
fKFVertex, 3e6);
1399 TVector3 refittedMomentum_electron;
1400 const FairTrackParam* vtxTrack_electron =
1401 stsTracks_electron[0].GetParamFirst();
1402 vtxTrack_electron->Momentum(refittedMomentum_electron);
1403 float result_chi_electron = chiPrim_electron[0];
1404 float result_ndf_electron = stsTracks_electron[0].GetNDF();
1407 Double_t startvertexZ = vtxTrack_electron->GetZ();
1414 Double_t nofhits_sts = stsTrack->
GetNofHits();
1418 result_chi_electron);
1420 result_chi_electron);
1432 if (richInd < 0)
continue;
1439 iGTrack, refittedMomentum_electron.Mag());
1441 if (TMath::Abs(pdg) == 11) {
1444 refittedMomentum_electron.Perp())) {
1448 if (TMath::Abs(pdg) != 11) {
1451 refittedMomentum_electron.Perp())) {
1458 TVector3 stsMomentumVec;
1460 Double_t stsMomentum = stsMomentumVec.Mag();
1462 TVector3 mcMomentumVec;
1464 Double_t mcMomentum = mcMomentumVec.Mag();
1471 bothtogether.SetX(mcMomentumVec.X());
1472 bothtogether.SetY(stsMomentumVec.Y());
1473 bothtogether.SetZ(stsMomentumVec.Z());
1483 if (isFilled) nofElectrons4epem++;
1488 if (richMatch == NULL)
continue;
1490 if (richMcTrackId < 0)
continue;
1492 if (mcTrack2 == NULL)
continue;
1497 if (stsMcTrackId == richMcTrackId) {
1501 iGTrack, refittedMomentum_electron.Mag());
1502 if (TMath::Abs(pdg) == 11) {
1507 if (TMath::Abs(pdg) != 11) {
1533 result_chi_electron);
1535 result_chi_electron);
1577 cout <<
"\n\n############### CALLING FINISH ROUTINES... ############" << endl;
1581 gDirectory->mkdir(
"analysis-conversion");
1582 gDirectory->cd(
"analysis-conversion");
1585 gDirectory->mkdir(
"KFParticle");
1586 gDirectory->cd(
"KFParticle");
1590 gDirectory->cd(
"..");
1604 gDirectory->mkdir(
"further analyses");
1605 gDirectory->cd(
"further analyses");
1609 gDirectory->cd(
"..");
1614 gDirectory->cd(
"..");
1621 cout <<
"############### FINISHED MAIN TASK ##############" << endl;
1630 cout <<
"#####################################" << endl;
1633 cout <<
"Number of reconstructed pi0 (all events): "
1637 cout <<
"#####################################" << endl;
1638 cout <<
"############### OVERALL TIMERS ###############" << endl;
1640 cout << std::setprecision(1);
1641 cout <<
"Complete time: " <<
fTime_all << endl;
1644 cout <<
"############### ############## ###############" << endl;
1645 cout <<
"Number of events in fhNofPi0_perEvent histogram: "
1647 cout <<
"############### ############## ###############" << endl;
1656 if (motherId == -1)
return;
1658 int mcMotherPdg = -1;
1659 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1662 if (mcMotherPdg == 22) {
1665 if (grandmotherId == -1)
return;
1667 int mcGrandmotherPdg = -1;
1668 if (NULL != grandmother) mcGrandmotherPdg = grandmother->
GetPdgCode();
1672 if (mcGrandmotherPdg == 111) {
1680 if (mcMotherPdg != 22 && mcMotherPdg != 111 && mcMotherPdg != 221)
1683 if (mcMotherPdg == 22) {
1699 Double_t energyP = TMath::Sqrt(momP.Mag2() +
M2E);
1700 TLorentzVector lorVecP(momP, energyP);
1704 Double_t energyM = TMath::Sqrt(momM.Mag2() +
M2E);
1705 TLorentzVector lorVecM(momM, energyM);
1707 TVector3 momPair = momP + momM;
1708 Double_t energyPair = energyP + energyM;
1709 Double_t ptPair = momPair.Perp();
1710 Double_t pzPair = momPair.Pz();
1712 0.5 * TMath::Log((energyPair + pzPair) / (energyPair - pzPair));
1713 Double_t anglePair = lorVecM.Angle(lorVecP.Vect());
1714 Double_t theta = 180. * anglePair /
TMath::Pi();
1716 2. * TMath::Sin(anglePair / 2.) * TMath::Sqrt(momM.Mag() * momP.Mag());
1719 params.
fPt = ptPair;
1721 params.
fMinv = minv;
1732 Double_t energy1 = TMath::Sqrt(mom1.Mag2());
1733 TLorentzVector lorVec1(mom1, energy1);
1737 Double_t energy2 = TMath::Sqrt(mom2.Mag2());
1738 TLorentzVector lorVec2(mom2, energy2);
1740 TLorentzVector sum = lorVec1 + lorVec2;
1751 Double_t energy1 = TMath::Sqrt(mom1.Mag2() +
M2E);
1752 TLorentzVector lorVec1(mom1, energy1);
1756 Double_t energy2 = TMath::Sqrt(mom2.Mag2() +
M2E);
1757 TLorentzVector lorVec2(mom2, energy2);
1759 TLorentzVector sum = lorVec1 + lorVec2;
1806 TLorentzVector lorVec1;
1809 TLorentzVector lorVec2;
1812 TLorentzVector lorVec3;
1815 TLorentzVector lorVec4;
1820 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1821 cout <<
"mc: \t" << sum.Px() <<
" / " << sum.Py() <<
" / " << sum.Pz()
1822 <<
" / " << sum.E() <<
"\t => mag = " << sum.Mag() << endl;
1828 const TVector3 part2,
1829 const TVector3 part3,
1830 const TVector3 part4)
1833 Double_t energy1 = TMath::Sqrt(part1.Mag2() +
M2E);
1834 TLorentzVector lorVec1(part1, energy1);
1836 Double_t energy2 = TMath::Sqrt(part2.Mag2() +
M2E);
1837 TLorentzVector lorVec2(part2, energy2);
1839 Double_t energy3 = TMath::Sqrt(part3.Mag2() +
M2E);
1840 TLorentzVector lorVec3(part3, energy3);
1842 Double_t energy4 = TMath::Sqrt(part4.Mag2() +
M2E);
1843 TLorentzVector lorVec4(part4, energy4);
1846 sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
1860 Bool_t electrons =
true;
1861 Bool_t gammas =
true;
1864 if (TMath::Abs(mctrack->
GetPdgCode()) == 22 && gammas) {
1866 if (motherId != -1) {
1868 int mcMotherPdg = -1;
1869 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1870 if (mcMotherPdg == 111) {
1877 if (TMath::Abs(mctrack->
GetPdgCode()) == 11 && electrons) {
1882 if (motherId != -1 || motherId == -1) {
1884 int mcMotherPdg = -1;
1885 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1886 if (mcMotherPdg == 22 || mcMotherPdg == 111 || mcMotherPdg == 221) {
1888 if (grandmotherId != -1 || grandmotherId == -1) {
1904 if (TMath::Abs(mctrack->
GetPdgCode()) == 11) {
1906 if (motherId != -1) {
1908 int mcMotherPdg = -1;
1909 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1910 if (mcMotherPdg == 111
1911 || mcMotherPdg == 22) {
1922 TVector3 stsMomentum,
1923 TVector3 refittedMom,
1926 Int_t GlobalTrackId) {
1927 Bool_t isFilled =
false;
1928 if (TMath::Abs(mctrack->
GetPdgCode()) == 11) {
1930 if (motherId != -1) {
1953 int mcMotherPdg_i = -1;
1954 if (NULL != mother_i) mcMotherPdg_i = mother_i->
GetPdgCode();
1960 if (motherId_i == motherId_j
1968 if (mcMotherPdg_i == 111) {
1970 cout <<
"-- mother particle of decay: pi0" << endl;
1982 for (
unsigned int j =
i + 1; j <
fMCTracklist.size(); j++) {
1995 if (motherId_i == motherId_j) {
1998 if (invmass < 0.001) {
2003 if (vi.Z() == vj.Z() && vi.Y() == vj.Y() && vi.X() == vj.X()) {
2034 int grandmotherId1 = -1;
2035 int grandmotherId2 = -1;
2036 int grandmotherId3 = -1;
2037 int grandmotherId4 = -1;
2039 int mcMotherPdg1 = -1;
2040 int mcMotherPdg2 = -1;
2041 int mcMotherPdg3 = -1;
2050 if (motherId1 != -1) {
2052 if (NULL != mother1) mcMotherPdg1 = mother1->
GetPdgCode();
2054 if (grandmotherId1 != -1) {
2059 if (motherId2 != -1) {
2061 if (NULL != mother2) mcMotherPdg2 = mother2->
GetPdgCode();
2063 if (grandmotherId2 != -1) {
2068 if (motherId3 != -1) {
2070 if (NULL != mother3) mcMotherPdg3 = mother3->
GetPdgCode();
2072 if (grandmotherId3 != -1) {
2077 if (motherId4 != -1) {
2081 if (grandmotherId4 != -1) {
2088 if (motherId1 == motherId2 && motherId2 == motherId3
2089 && motherId3 == motherId4) {
2098 if (motherId1 == motherId2 && motherId1 == motherId3) {}
2099 if (motherId1 == motherId2 && motherId1 == motherId4) {}
2100 if (motherId1 == motherId3 && motherId1 == motherId4) {}
2101 if (motherId2 == motherId3 && motherId2 == motherId4) {}
2106 if (((motherId1 == motherId2 && motherId3 == motherId4)
2107 && (mcMotherPdg1 == 22) && (mcMotherPdg3 == 22)
2108 && grandmotherId1 == grandmotherId3)
2109 || ((motherId1 == motherId3 && motherId2 == motherId4)
2110 && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
2111 && grandmotherId1 == grandmotherId2)
2112 || ((motherId1 == motherId4 && motherId2 == motherId3)
2113 && (mcMotherPdg1 == 22) && (mcMotherPdg2 == 22)
2114 && grandmotherId1 == grandmotherId2)) {
2159 if (motherId1 == motherId2) {
2176 cout <<
"kf works" << endl;
2178 cout <<
"kf works not" << endl;
2184 Bool_t reconstructible =
false;
2188 Int_t daughters = 0;
2189 Int_t gammaDaughters = 0;
2190 Int_t electronDaughters = 0;
2191 Int_t electronDaughtersBeforeRICH = 0;
2192 Int_t electronDaughtersInTarget = 0;
2193 vector<int> gammaDaughterIDs;
2194 gammaDaughterIDs.clear();
2195 vector<int> electronDaughterIDs;
2196 electronDaughterIDs.clear();
2198 Int_t nofMcTracks =
fMcTracks->GetEntriesFast();
2199 for (
int i = 0;
i < nofMcTracks;
i++) {
2201 if (mctrack_switch == NULL)
continue;
2204 if (motherID == trackid && pdg == 22) {
2206 gammaDaughterIDs.push_back(
i);
2208 if (motherID == trackid) { daughters++; }
2220 if (gammaDaughters == 2) {
2223 for (
int i = 0;
i < nofMcTracks;
i++) {
2225 if (mctrack_switch == NULL)
continue;
2228 if (TMath::Abs(pdg) == 11
2229 && (motherID == gammaDaughterIDs[0]
2230 || motherID == gammaDaughterIDs[1])) {
2231 electronDaughters++;
2232 electronDaughterIDs.push_back(
i);
2234 TVector3 startvertex;
2236 if (startvertex.Z() <= 70) { electronDaughtersBeforeRICH++; }
2237 if (startvertex.Z() <= 4) { electronDaughtersInTarget++; }
2242 if (electronDaughtersBeforeRICH == 4) {
2244 reconstructible =
true;
2250 Int_t daughters = 0;
2251 Int_t gammaDaughters = 0;
2252 Int_t electronDaughters = 0;
2253 Int_t electronDaughtersBeforeRICH = 0;
2254 Int_t electronDaughtersInTarget = 0;
2255 vector<int> gammaDaughterIDs;
2256 gammaDaughterIDs.clear();
2257 vector<int> electronDaughterIDs;
2258 electronDaughterIDs.clear();
2260 Int_t nofMcTracks =
fMcTracks->GetEntriesFast();
2261 for (
int i = 0;
i < nofMcTracks;
i++) {
2263 if (mctrack_switch == NULL)
continue;
2266 if (motherID == trackid && pdg == 22) {
2268 gammaDaughterIDs.push_back(
i);
2270 if (motherID == trackid) { daughters++; }
2273 if (gammaDaughters == 2) {
2275 for (
int i = 0;
i < nofMcTracks;
i++) {
2277 if (mctrack_switch == NULL)
continue;
2280 if (TMath::Abs(pdg) == 11
2281 && (motherID == gammaDaughterIDs[0]
2282 || motherID == gammaDaughterIDs[1])) {
2283 electronDaughters++;
2284 electronDaughterIDs.push_back(
i);
2286 TVector3 startvertex;
2288 if (startvertex.Z() <= 70) { electronDaughtersBeforeRICH++; }
2289 if (startvertex.Z() <= 4) { electronDaughtersInTarget++; }
2298 return reconstructible;
2303 Int_t electrons = 0;
2304 std::multimap<int, int>
2306 std::multimap<int, int> electronPi0ID_map_richInd;
2309 for (
int i = 0;
i < nof;
i++) {
2313 if (TMath::Abs(pdg) == 11 && motherID != -1) {
2315 if (mctrack_mother == NULL)
continue;
2316 Int_t grandmotherID = mctrack_mother->
GetMotherId();
2317 Int_t motherPdg = mctrack_mother->
GetPdgCode();
2319 if (motherPdg == 22 && grandmotherID != -1) {
2322 if (mctrack_grandmother == NULL)
continue;
2323 Int_t grandmotherPdg = mctrack_grandmother->
GetPdgCode();
2325 if (grandmotherPdg == 111) {
2327 electronPi0ID_map.insert(std::pair<int, int>(grandmotherID,
i));
2329 electronPi0ID_map_richInd.insert(
2330 std::pair<int, int>(grandmotherID,
i));
2337 TVector3 startvertex;
2339 if (startvertex.Z() < 1) {
2357 int samePi0counter = 0;
2359 for (std::map<int, int>::iterator it = electronPi0ID_map.begin();
2360 it != electronPi0ID_map.end();
2362 if (it == electronPi0ID_map.begin()) check = 1;
2363 if (it != electronPi0ID_map.begin()) {
2364 std::map<int, int>::iterator zwischen = it;
2367 int id_old = zwischen->first;
2382 std::map<int, int>::iterator alt3 = zwischen;
2384 std::map<int, int>::iterator alt4 = alt3;
2386 std::map<int, int>::iterator alt5 = alt4;
2388 cout <<
"CbmAnaConversion: AnalysePi0_Reco: " << alt5->first <<
"/"
2389 << zwischen->first <<
"/" << alt3->first <<
"/" << alt4->first
2391 cout <<
"CbmAnaConversion: AnalysePi0_Reco: " << alt5->second <<
"/"
2392 << zwischen->second <<
"/" << alt3->second <<
"/" << alt4->second
2395 alt5->second, zwischen->second, alt3->second, alt4->second);
2410 Double_t invmass_mc =
2419 std::map<int, int>::iterator temp = zwischen;
2420 for (
int i = 0;
i < 4;
i++) {
2428 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
2436 std::map<int, int>::iterator temp2 = zwischen;
2437 for (
int i = 0;
i < 4;
i++) {
2457 if (samePi0counter >= 4) {
2462 int samePi0counter2 = 0;
2464 cout <<
"CbmAnaConversion: RecoPi0: electronmapsize: "
2465 << electronPi0ID_map_richInd.size() << endl;
2466 for (std::map<int, int>::iterator it = electronPi0ID_map_richInd.begin();
2467 it != electronPi0ID_map_richInd.end();
2469 if (it == electronPi0ID_map_richInd.begin()) check = 1;
2470 if (it != electronPi0ID_map_richInd.begin()) {
2471 std::map<int, int>::iterator zwischen = it;
2474 int id_old = zwischen->first;
2533 Bool_t IsWithinCuts =
false;
2534 Double_t OpeningAngleCut = 1.5;
2536 if (motherID_e1 == motherID_e2 && motherID_e3 == motherID_e4) {
2537 Double_t anglePair1 = lorVec1.Angle(lorVec2.Vect());
2538 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2540 Double_t anglePair2 = lorVec3.Angle(lorVec4.Vect());
2541 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2543 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2544 IsWithinCuts =
true;
2547 cout <<
"CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 <<
"/"
2551 if (motherID_e1 == motherID_e3 && motherID_e2 == motherID_e4) {
2552 Double_t anglePair1 = lorVec1.Angle(lorVec3.Vect());
2553 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2555 Double_t anglePair2 = lorVec2.Angle(lorVec4.Vect());
2556 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2558 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2559 IsWithinCuts =
true;
2562 cout <<
"CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 <<
"/"
2566 if (motherID_e1 == motherID_e4 && motherID_e2 == motherID_e3) {
2567 Double_t anglePair1 = lorVec1.Angle(lorVec4.Vect());
2568 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2570 Double_t anglePair2 = lorVec2.Angle(lorVec3.Vect());
2571 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2573 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2574 IsWithinCuts =
true;
2577 cout <<
"CbmAnaConversion: AnalysePi0_Reco_calc: " << theta1 <<
"/"
2583 return IsWithinCuts;
2588 Int_t electrons = 0;
2589 std::multimap<int, int>
2591 std::multimap<int, int>
2592 electronPi0ID_map_richInd;
2595 for (
int i = 0;
i < nof;
i++) {
2599 if (TMath::Abs(pdg) == 11 && motherID != -1) {
2601 if (mctrack_mother == NULL)
continue;
2602 Int_t grandmotherID = mctrack_mother->
GetMotherId();
2603 Int_t motherPdg = mctrack_mother->
GetPdgCode();
2605 if (motherPdg == 22 && grandmotherID != -1) {
2608 if (mctrack_grandmother == NULL)
continue;
2609 Int_t grandmotherPdg = mctrack_grandmother->
GetPdgCode();
2611 if (grandmotherPdg == 111) {
2613 electronPi0ID_map.insert(std::pair<int, int>(grandmotherID,
i));
2615 electronPi0ID_map_richInd.insert(
2616 std::pair<int, int>(grandmotherID,
i));
2625 TVector3 startvertex;
2627 if (startvertex.Z() < 1) {
2652 if (TMath::Abs(pdg) != 11)
continue;
2669 int samePi0counter = 0;
2671 int check_withRichInd = 0;
2672 for (std::map<int, int>::iterator it = electronPi0ID_map.begin();
2673 it != electronPi0ID_map.end();
2675 if (it == electronPi0ID_map.begin()) check = 1;
2676 if (it != electronPi0ID_map.begin()) {
2677 std::map<int, int>::iterator zwischen = it;
2680 int id_old = zwischen->first;
2688 check_withRichInd++;
2698 std::map<int, int>::iterator alt3 = zwischen;
2700 std::map<int, int>::iterator alt4 = alt3;
2702 std::map<int, int>::iterator alt5 = alt4;
2704 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd: " << alt5->first
2705 <<
"/" << zwischen->first <<
"/" << alt3->first <<
"/"
2706 << alt4->first << endl;
2707 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd: "
2708 << alt5->second <<
"/" << zwischen->second <<
"/" << alt3->second
2709 <<
"/" << alt4->second << endl;
2711 alt5->second, zwischen->second, alt3->second, alt4->second);
2731 Double_t invmass_mc =
2740 std::map<int, int>::iterator temp = zwischen;
2741 for (
int i = 0;
i < 4;
i++) {
2748 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
2756 std::map<int, int>::iterator temp2 = zwischen;
2757 for (
int i = 0;
i < 4;
i++) {
2769 if (richInd_e1 >= 0 && richInd_e2 >= 0 && richInd_e3 >= 0
2770 && richInd_e4 >= 0) {
2776 if (chi_e1 <= 3 && chi_e2 <= 3 && chi_e3 <= 3 && chi_e4 <= 3) {
2786 std::map<int, int>::iterator temp = zwischen;
2787 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd, check>4: ";
2788 for (
int i = 0;
i < check;
i++) {
2789 TVector3 momentum_mc;
2791 Double_t theta_mc = 180.0 * momentum_mc.Theta() /
TMath::Pi();
2792 Double_t phi_mc = 180.0 * momentum_mc.Phi() /
TMath::Pi();
2793 TVector3 momentum_reco =
2795 Double_t theta_reco = 180.0 * momentum_reco.Theta() /
TMath::Pi();
2796 Double_t phi_reco = 180.0 * momentum_reco.Phi() /
TMath::Pi();
2797 cout <<
"(" << temp->first <<
"/" << temp->second <<
"/"
2799 << momentum_reco.Mag() <<
"/" << theta_mc <<
"-" << phi_mc
2800 <<
"/" << theta_reco <<
"-" << phi_reco <<
") \t";
2819 check_withRichInd = 1;
2821 check_withRichInd = 0;
2827 if (samePi0counter >= 4) {
2865 Bool_t IsWithinCuts =
false;
2866 Double_t OpeningAngleCut = 1.5;
2868 if (motherID_e1 == motherID_e2 && motherID_e3 == motherID_e4) {
2869 Double_t anglePair1 = lorVec1.Angle(lorVec2.Vect());
2870 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2872 Double_t anglePair2 = lorVec3.Angle(lorVec4.Vect());
2873 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2875 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2876 IsWithinCuts =
true;
2877 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1
2878 <<
"/" << theta2 << endl;
2881 if (motherID_e1 == motherID_e3 && motherID_e2 == motherID_e4) {
2882 Double_t anglePair1 = lorVec1.Angle(lorVec3.Vect());
2883 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2885 Double_t anglePair2 = lorVec2.Angle(lorVec4.Vect());
2886 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2888 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2889 IsWithinCuts =
true;
2890 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1
2891 <<
"/" << theta2 << endl;
2894 if (motherID_e1 == motherID_e4 && motherID_e2 == motherID_e3) {
2895 Double_t anglePair1 = lorVec1.Angle(lorVec4.Vect());
2896 Double_t theta1 = 180. * anglePair1 /
TMath::Pi();
2898 Double_t anglePair2 = lorVec2.Angle(lorVec3.Vect());
2899 Double_t theta2 = 180. * anglePair2 /
TMath::Pi();
2901 if (theta1 <= OpeningAngleCut && theta2 <= OpeningAngleCut)
2902 IsWithinCuts =
true;
2903 cout <<
"CbmAnaConversion: AnalysePi0_Reco_noRichInd_calc: " << theta1
2904 <<
"/" << theta2 << endl;
2909 return IsWithinCuts;