13 #include "FairEventHeader.h"
14 #include "FairMCPoint.h"
15 #include "FairRootManager.h"
23 #include "CbmTofHit.h"
27 #include "FairTrackParam.h"
34 #include "TClonesArray.h"
39 #include "TLorentzVector.h"
41 #include "TObjArray.h"
47 #include "FairBaseParSet.h"
48 #include "FairGeoMedium.h"
49 #include "FairGeoNode.h"
50 #include "FairGeoTransform.h"
51 #include "FairGeoVector.h"
52 #include "FairGeoVolume.h"
53 #include "FairRunAna.h"
54 #include "FairRuntimeDb.h"
55 #include "TDatabasePDG.h"
64 #include "TStopwatch.h"
90 hist[
i] =
new TH1D(hname.c_str(), hname.c_str(), nBins,
min,
max);
91 hist[
i]->GetXaxis()->SetTitle(axisX.c_str());
92 hist[
i]->GetYaxis()->SetTitle(axisY.c_str());
93 fHistoList.push_back(hist[
i]);
113 hname.c_str(), hname.c_str(), nBinsX, minX, maxX, nBinsY, minY, maxY);
114 hist[
i]->GetXaxis()->SetTitle(axisX.c_str());
115 hist[
i]->GetYaxis()->SetTitle(axisY.c_str());
116 hist[
i]->GetZaxis()->SetTitle(axisZ.c_str());
117 fHistoList.push_back(hist[
i]);
132 hist[
i] =
new TH1D(hname.c_str(), hname.c_str(), nBins,
min,
max);
133 hist[
i]->GetXaxis()->SetTitle(axisX.c_str());
134 hist[
i]->GetYaxis()->SetTitle(axisY.c_str());
135 fHistoList.push_back(hist[
i]);
155 hname.c_str(), hname.c_str(), nBinsX, minX, maxX, nBinsY, minY, maxY);
156 hist[
i]->GetXaxis()->SetTitle(axisX.c_str());
157 hist[
i]->GetYaxis()->SetTitle(axisY.c_str());
158 hist[
i]->GetZaxis()->SetTitle(axisZ.c_str());
159 fHistoList.push_back(hist[
i]);
164 : FairTask(
"CbmAnaDielectronTask")
165 , fMCEventHeader(NULL)
170 , fRichRingMatches(NULL)
172 , fGlobalTracks(NULL)
174 , fStsTrackMatches(NULL)
178 , fMvdHitMatches(NULL)
181 , fTrdTrackMatches(NULL)
183 , fTofHitsMatches(NULL)
199 , fPionMisidLevel(-1.)
200 , fRandom3(new TRandom3(0))
203 , fNofHitsInRingMap()
204 , fh_mc_signal_mom_angle()
205 , fh_nof_charged_particles()
206 , fh_nof_charged_particles_acc()
207 , fh_mc_mother_pdg(NULL)
208 , fh_acc_mother_pdg(NULL)
209 , fh_signal_pmtXY(NULL)
211 , fh_gamma_pmtXY(NULL)
212 , fh_vertex_el_gamma_xz()
213 , fh_vertex_el_gamma_yz()
214 , fh_vertex_el_gamma_xy()
215 , fh_vertex_el_gamma_rz()
223 , fh_signal_minv_pt()
226 , fh_bg_truematch_minv()
227 , fh_bg_truematch_el_minv()
228 , fh_bg_truematch_notel_minv()
229 , fh_bg_mismatch_minv()
230 , fh_source_bg_minv()
244 , fh_ttcut_truepair()
246 , fh_stcut_truepair()
248 , fh_rtcut_truepair()
255 , fh_mvd1cut_mc_dist_gamma(NULL)
256 , fh_mvd1cut_mc_dist_pi0(NULL)
257 , fh_mvd2cut_mc_dist_gamma(NULL)
258 , fh_mvd2cut_mc_dist_pi0(NULL)
261 , fh_source_pairs_epem()
262 , fh_source_pairs(NULL)
263 , fh_event_number(NULL)
264 , fh_nof_bg_tracks(NULL)
265 , fh_nof_el_tracks(NULL)
266 , fh_source_tracks(NULL)
267 , fh_nof_topology_pairs_gamma(NULL)
268 , fh_nof_topology_pairs_pi0(NULL)
269 , fh_nof_rec_pairs_gamma(NULL)
270 , fh_nof_rec_pairs_pi0(NULL)
271 , fh_nof_rec_gamma(NULL)
272 , fh_nof_rec_pi0(NULL)
273 , fh_nof_mismatches(NULL)
274 , fh_nof_mismatches_rich(NULL)
275 , fh_nof_mismatches_trd(NULL)
276 , fh_nof_mismatches_tof(NULL)
277 , fh_nof_ghosts(NULL)
282 , fh_pi_mom_acc(NULL)
283 , fh_pi_mom_rec(NULL)
284 , fh_pi_mom_rec_only_sts(NULL)
285 , fh_pi_mom_rec_sts_rich_trd(NULL)
286 , fh_pi_mom_rec_sts_rich_trd_tof(NULL)
287 , fh_pi_rapidity_mc(NULL)
288 , fh_piprim_mom_mc(NULL)
289 , fh_piprim_mom_acc(NULL)
290 , fh_piprim_mom_rec(NULL)
291 , fh_piprim_mom_rec_only_sts(NULL)
292 , fh_piprim_mom_rec_sts_rich_trd(NULL)
293 , fh_piprim_mom_rec_sts_rich_trd_tof(NULL)
294 , fh_piprim_plus_rapidity_mc(NULL)
295 , fh_piprim_minus_rapidity_mc(NULL)
296 , fh_pi0prim_rapidity_mc(NULL)
297 , fh_etaprim_rapidity_mc(NULL) {
317 new TH2D(
"fh_mc_signal_mom_angle",
318 "fh_mc_signal_mom_angle; #sqrt{p_{e^{#pm}} p_{e^{#mp}}} [GeV/c]; "
319 "#theta_{e^{+},e^{-}} [deg] ;Counter",
330 new TH1D(
"fh_nof_charged_particles",
331 "fh_nof_charged_particles; nof charged particles; Yield",
337 new TH1D(
"fh_nof_charged_particles_acc",
338 "fh_nof_charged_particles_acc; nof charged particles; Yield",
346 "fh_mc_mother_pdg; Pdg code; Particles per event",
352 new TH1D(
"fh_acc_mother_pdg",
353 "fh_acc_mother_pdg; Pdg code; Particles per event",
361 "fh_signal_pmtXY;X [cm];Y [cm];Counter",
370 "fh_pi0_pmtXY;X [cm];Y [cm];Counter",
379 "fh_gamma_pmtXY;X [cm];Y [cm];Counter",
390 "fh_vertex_el_gamma_xz",
401 "fh_vertex_el_gamma_yz",
412 "fh_vertex_el_gamma_xy",
423 "fh_vertex_el_gamma_rz",
425 "#sqrt{X^{2}+Y^{2}} [cm]",
436 "fh_nof_bg_tracks;Analysis steps;Tracks/event",
442 "fh_nof_el_tracks;Analysis steps;Tracks/event",
448 "fh_source_tracks;Analysis steps;Particle",
458 new TH1D(
"fh_nof_topology_pairs_gamma",
459 "fh_nof_topology_pairs_gamma;Pair type;Pairs/event",
466 new TH1D(
"fh_nof_topology_pairs_pi0",
467 "fh_nof_topology_pairs_pi0;Pair type;Pairs/event",
475 "fh_nof_mismatches;Analysis steps;Tracks/event",
481 new TH1D(
"fh_nof_mismatches_rich",
482 "fh_nof_mismatches_rich;Analysis steps;Tracks/event",
488 new TH1D(
"fh_nof_mismatches_trd",
489 "fh_nof_mismatches_trd;Analysis steps;Tracks/event",
495 new TH1D(
"fh_nof_mismatches_tof",
496 "fh_nof_mismatches_tof;Analysis steps;Tracks/event",
502 "fh_nof_ghosts;Analysis steps;Tracks/event",
510 "fh_source_pairs;Analysis steps;Pair",
520 fh_event_number =
new TH1D(
"fh_event_number",
"fh_event_number", 1, 0, 1.);
524 fh_richann,
"fh_richann",
"ANN output",
"Yield", 100, -1.1, 1.1);
526 fh_trdann,
"fh_trdann",
"ANN output",
"Yield", 100, -1.1, 1.1);
530 "m^{2} [GeV/c^{2}]^{2}",
546 fh_chi2sts,
"fh_chi2sts",
"#chi^{2}",
"Yield", 200, 0., 20.);
549 fh_chi2prim,
"fh_chi2prim",
"#chi^{2}_{prim}",
"Yield", 200, 0., 20.);
553 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
554 "#theta_{e^{+},e^{-}} [deg]",
565 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
566 "#theta_{e^{#pm},rec} [deg]",
577 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
578 "#theta_{e^{#pm},rec} [deg]",
613 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
614 "#theta_{e^{+},e^{-}} [deg]",
624 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
625 "#theta_{e^{+},e^{-}} [deg]",
635 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
636 "#theta_{e^{+},e^{-}} [deg]",
646 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
647 "#theta_{e^{+},e^{-}} [deg]",
657 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
658 "#theta_{e^{+},e^{-}} [deg]",
668 "#sqrt{p_{e^{#pm}} p_{rec}} [GeV/c]",
669 "#theta_{e^{+},e^{-}} [deg]",
680 "Number of hits in MVD",
687 "Number of hits in STS",
704 fh_mvd1r,
"fh_mvd1r",
"#sqrt{X^{2}+Y^{2}} [cm]",
"Yield", 60, 0., 3.);
717 fh_mvd2r,
"fh_mvd2r",
"#sqrt{X^{2}+Y^{2}} [cm]",
"Yield", 60, 0., 6.);
721 fh_mvd1cut_qa,
"fh_mvd1cut_qa",
"MVD hit assignment",
"Yield", 2, 0., 2.);
723 fh_mvd2cut_qa,
"fh_mvd2cut_qa",
"MVD hit assignment",
"Yield", 2, 0., 2.);
728 "M_{ee} [GeV/c^{2}]",
734 fh_bg_minv,
"fh_bg_minv",
"M_{ee} [GeV/c^{2}]",
"Yield", 4000, 0., 4.);
736 fh_pi0_minv,
"fh_pi0_minv",
"M_{ee} [GeV/c^{2}]",
"Yield", 4000, 0., 4.);
738 fh_eta_minv,
"fh_eta_minv",
"M_{ee} [GeV/c^{2}]",
"Yield", 4000, 0., 4.);
741 "M_{ee} [GeV/c^{2}]",
748 "fh_bg_truematch_minv",
749 "M_{ee} [GeV/c^{2}]",
755 "fh_bg_truematch_el_minv",
756 "M_{ee} [GeV/c^{2}]",
762 "fh_bg_truematch_notel_minv",
763 "M_{ee} [GeV/c^{2}]",
769 "fh_bg_mismatch_minv",
770 "M_{ee} [GeV/c^{2}]",
779 ss <<
"fh_source_bg_minv_" <<
i;
782 "M_{ee} [GeV/c^{2}]",
791 "M_{ee} [GeV/c^{2}]",
802 "M_{ee} [GeV/c^{2}]",
813 "M_{ee} [GeV/c^{2}]",
825 fh_signal_mom,
"fh_signal_mom",
"P [GeV/c]",
"Yield", 100, 0., 15.);
840 "fh_source_pairs_epem",
841 "mother particle e+",
842 "mother particle e-",
862 string hname =
"", htitle =
"";
865 htitle = hname +
";#Theta_{1,2} [deg];Yield";
867 new TH1D(hname.c_str(), htitle.c_str(), 160, 0., 80.);
872 htitle = hname +
";P [GeV/c];Yield";
874 new TH1D(hname.c_str(), htitle.c_str(), 300, 0., 15.);
879 htitle = hname +
";P_{t} [GeV/c];Yield";
881 new TH1D(hname.c_str(), htitle.c_str(), 100, 0., 5.);
888 "fh_pi_mom_mc",
"fh_pi_mom_mc;p [GeV/c];dN/dP [1/GeV/c]", 30, 0., 3.);
891 "fh_pi_mom_acc",
"fh_pi_mom_acc;p [GeV/c];dN/dP [1/GeV/c]", 30, 0., 3.);
894 "fh_pi_mom_rec",
"fh_pi_mom_rec;p [GeV/c];dN/dP [1/GeV/c]", 30, 0., 3.);
897 new TH1D(
"fh_pi_mom_rec_only_sts",
898 "fh_pi_mom_rec_only_sts;p [GeV/c];dN/dP [1/GeV/c]",
904 new TH1D(
"fh_pi_mom_rec_sts_rich_trd",
905 "fh_pi_mom_rec_sts_rich_trd;p [GeV/c];dN/dP [1/GeV/c]",
911 new TH1D(
"fh_pi_mom_rec_sts_rich_trd_tof",
912 "fh_pi_mom_rec_sts_rich_trd_tof;p [GeV/c];dN/dP [1/GeV/c]",
918 "fh_pi_rapidity_mc",
"fh_pi_rapidity_mc;Rapidity;dN/dY", 400, 0., 4.);
924 "fh_piprim_mom_mc;p [GeV/c];dN/dP [1/GeV/c]",
930 "fh_piprim_mom_acc;p [GeV/c];dN/dP [1/GeV/c]",
936 "fh_piprim_mom_rec;p [GeV/c];dN/dP [1/GeV/c]",
942 new TH1D(
"fh_piprim_mom_rec_only_sts",
943 "fh_piprim_mom_rec_only_sts;p [GeV/c];dN/dP [1/GeV/c]",
949 new TH1D(
"fh_piprim_mom_rec_sts_rich_trd",
950 "fh_piprim_mom_rec_sts_rich_trd;p [GeV/c];dN/dP [1/GeV/c]",
956 new TH1D(
"fh_piprim_mom_rec_sts_rich_trd_tof",
957 "fh_piprim_mom_rec_sts_rich_trd_tof;p [GeV/c];dN/dP [1/GeV/c]",
964 new TH1D(
"fh_piprim_plus_rapidity_mc",
965 "fh_piprim_plus_rapidity_mc;Rapidity;dN/dY",
971 new TH1D(
"fh_piprim_minus_rapidity_mc",
972 "fh_piprim_minus_rapidity_mc;Rapidity;dN/dY",
978 "fh_pi0prim_rapidity_mc;Rapidity;dN/dY",
984 "fh_etaprim_rapidity_mc;Rapidity;dN/dY",
992 new TH1D(
"fh_nof_rec_pairs_gamma",
993 "fh_nof_rec_pairs_gamma;Pair category; Number per event",
999 new TH1D(
"fh_nof_rec_pairs_pi0",
1000 "fh_nof_rec_pairs_pi0;Pair category; Number per event",
1007 new TH1D(
"fh_nof_rec_gamma",
1008 "fh_nof_rec_gamma;Track category; Number per event",
1014 "fh_nof_rec_pi0;Track category; Number per event",
1022 cout <<
"InitStatus CbmAnaDielectronTask::Init" << endl;
1024 FairRootManager* ioman = FairRootManager::Instance();
1025 if (NULL == ioman) {
1026 Fatal(
"CbmAnaDielectronTask::Init",
"No FairRootManager!");
1029 fMCEventHeader = (FairMCEventHeader*) ioman->GetObject(
"MCEventHeader.");
1031 Fatal(
"CbmAnaDielectronTask::Init",
"No MCEventHeader array!");
1034 fMCTracks = (TClonesArray*) ioman->GetObject(
"MCTrack");
1036 Fatal(
"CbmAnaDielectronTask::Init",
"No MCTrack array!");
1040 fRichHits = (TClonesArray*) ioman->GetObject(
"RichHit");
1042 Fatal(
"CbmAnaDielectronTask::Init",
"No RichHit array!");
1045 fRichRings = (TClonesArray*) ioman->GetObject(
"RichRing");
1047 Fatal(
"CbmAnaDielectronTask::Init",
"No RichRing array!");
1050 fRichPoints = (TClonesArray*) ioman->GetObject(
"RichPoint");
1052 Fatal(
"CbmAnaDielectronTask::Init",
"No RichPoint array!");
1057 Fatal(
"CbmAnaDielectronTask::Init",
"No RichRingMatch array!");
1060 fRichProj = (TClonesArray*) ioman->GetObject(
"RichProjection");
1062 Fatal(
"CbmAnaDielectronTask::Init",
"No RichProjection array!");
1068 Fatal(
"CbmAnaDielectronTask::Init",
"No StsTrackMatch array!");
1071 fStsTracks = (TClonesArray*) ioman->GetObject(
"StsTrack");
1073 Fatal(
"CbmAnaDielectronTask::Init",
"No StsTrack array!");
1076 fStsHits = (TClonesArray*) ioman->GetObject(
"StsHit");
1078 Fatal(
"CbmAnaDielectronTask::Init",
"No StsHit array!");
1082 fMvdHits = (TClonesArray*) ioman->GetObject(
"MvdHit");
1084 Fatal(
"CbmAnaDielectronTask::Init",
"No MvdHit array!");
1087 fMvdPoints = (TClonesArray*) ioman->GetObject(
"MvdPoint");
1089 Fatal(
"CbmAnaDielectronTask::Init",
": No MvdPoint array!");
1092 fMvdHitMatches = (TClonesArray*) ioman->GetObject(
"MvdHitMatch");
1094 Fatal(
"CbmAnaDielectronTask::Init",
": No MvdHitMatch array!");
1098 fGlobalTracks = (TClonesArray*) ioman->GetObject(
"GlobalTrack");
1100 Fatal(
"CbmAnaDielectronTask::Init",
"No GlobalTrack array!");
1104 fTrdTracks = (TClonesArray*) ioman->GetObject(
"TrdTrack");
1106 Fatal(
"CbmAnaDielectronTask::Init",
"No TrdTrack array!");
1111 Fatal(
"CbmAnaDielectronTask::Init",
"No TrdTrackMatch array!");
1116 fTofPoints = (TClonesArray*) ioman->GetObject(
"TofPoint");
1118 Fatal(
"CbmAnaDielectronTask::Init",
"No TofPoint array!");
1121 fTofHits = (TClonesArray*) ioman->GetObject(
"TofHit");
1123 Fatal(
"CbmAnaDielectronTask::Init",
"No TofHit array!");
1128 Fatal(
"CbmAnaDielectronTask::Init",
"No TofHitMatch Array! ");
1140 if (
nullptr ==
fPrimVertex) { LOG(fatal) <<
"No PrimaryVertex array!"; }
1159 Bool_t useMbias =
false;
1161 Bool_t isCentralCollision =
false;
1165 if (impactPar <= 7.7) isCentralCollision =
true;
1168 cout <<
"-I- CbmAnaDielectronTask, event number "
1172 cout <<
"fWeight = " <<
fWeight << endl;
1177 Fatal(
"CbmAnaDielectronTask::Exec",
"No PrimaryVertex array!");
1182 if (useMbias || (!useMbias && isCentralCollision)) {
1201 Int_t nofRichHits =
fRichHits->GetEntriesFast();
1202 for (Int_t iHit = 0; iHit < nofRichHits; iHit++) {
1204 if (NULL == hit)
continue;
1207 if (iPoint < 0)
continue;
1209 FairMCPoint* point =
static_cast<FairMCPoint*
>(
fRichPoints->At(iPoint));
1210 if (NULL == point)
continue;
1212 Int_t iMCTrack = point->GetTrackID();
1214 if (NULL == track)
continue;
1217 if (iMother == -1)
continue;
1225 Int_t nMcTracks =
fMCTracks->GetEntries();
1226 for (Int_t
i = 0;
i < nMcTracks;
i++) {
1229 Int_t pdg = TMath::Abs(mctrack->
GetPdgCode());
1230 Double_t mom = mctrack->
GetP();
1233 Bool_t isMcGammaElectron =
1235 if (isMcSignalElectron) {
1237 for (Int_t iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
1238 if (
i == iMc2)
continue;
1241 if (motherId == motherIdMc2
1246 Double_t angle = pKin.
fAngle;
1247 Double_t pMc = mctrack->
GetP();
1248 Double_t pMc2 = mctrack2->
GetP();
1249 Double_t sqrtPMc = TMath::Sqrt(pMc * pMc2);
1253 for (Int_t iMc2 = 0; iMc2 < nMcTracks; iMc2++) {
1254 if (
i == iMc2)
continue;
1257 if (motherId == motherIdMc2
1262 Double_t angle = pKin.
fAngle;
1263 Double_t pMc = mctrack->
GetP();
1264 Double_t pMc2 = mctrack2->
GetP();
1265 Double_t sqrtPMc = TMath::Sqrt(pMc * pMc2);
1270 if (isMcGammaElectron) {
1277 sqrt(
v.X() *
v.X() +
v.Y() *
v.Y()));
1281 Int_t mcMotherPdg = 0;
1283 if (motherId != -1) {
1285 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1295 Int_t nofRichHits =
fRichHits->GetEntriesFast();
1296 for (Int_t iH = 0; iH < nofRichHits; iH++) {
1298 if (richHit == NULL)
continue;
1299 Int_t pointInd = richHit->
GetRefId();
1300 if (pointInd < 0)
continue;
1302 FairMCPoint* pointPhoton =
1303 static_cast<FairMCPoint*
>(
fRichPoints->At(pointInd));
1304 if (NULL == pointPhoton)
continue;
1306 Int_t iMCTrackPhoton = pointPhoton->GetTrackID();
1309 if (NULL == trackPhoton)
continue;
1312 if (iMCTrack == -1)
continue;
1315 if (NULL == mctrack)
continue;
1319 Bool_t isPrim = (
v.Z() < 2.);
1322 Bool_t isMcGammaElectron =
1326 if (isMcSignalElectron) {
1329 if (isMcGammaElectron && isPrim) {
1332 if (isMcPi0Electron && isPrim) {
1339 Int_t nofMcTracks =
fMCTracks->GetEntries();
1340 Int_t nofChargedUrqmdParticles = 0;
1341 Int_t nofChargedUrqmdParticlesAcc = 0;
1342 for (Int_t
i = 0;
i < nofMcTracks;
i++) {
1344 Bool_t isMcTrackCharged =
false;
1345 if (mcTrack->
GetCharge() != 0) isMcTrackCharged =
true;
1347 Bool_t isMcTrackAccepted =
false;
1349 Bool_t isPrimary =
false;
1351 if (motherId == -1) isPrimary =
true;
1352 if (!isMcSignalElectron && isMcTrackCharged && isPrimary)
1353 nofChargedUrqmdParticles++;
1354 if (!isMcSignalElectron && isMcTrackCharged && isMcTrackAccepted
1356 nofChargedUrqmdParticlesAcc++;
1365 if (tr == NULL)
return false;
1374 Int_t nMcTracks =
fMCTracks->GetEntries();
1375 for (Int_t
i = 0;
i < nMcTracks;
i++) {
1378 Int_t pdg = TMath::Abs(mctrack->
GetPdgCode());
1383 Bool_t isAcc = (nMvdPoints + nStsPoints >= 4 && nRichPoints >= 7);
1384 Bool_t isMcGammaElectron =
1387 if (isMcGammaElectron) {
1394 sqrt(
v.X() *
v.X() +
v.Y() *
v.Y()));
1397 Int_t mcMotherPdg = 0;
1398 if (pdg == 11 && isAcc) {
1399 if (motherId != -1) {
1401 if (NULL != mother) mcMotherPdg = mother->
GetPdgCode();
1411 Int_t nMcTracks =
fMCTracks->GetEntries();
1412 for (Int_t iP = 0; iP < nMcTracks; iP++) {
1416 if (pdgP != 11)
continue;
1418 for (Int_t iM = 0; iM < nMcTracks; iM++) {
1419 if (iP == iM)
continue;
1423 if (pdgM != -11)
continue;
1441 if (isAccP && isAccM) {
1463 Int_t nMcTracks =
fMCTracks->GetEntries();
1464 for (Int_t
i = 0;
i < nMcTracks;
i++) {
1467 Int_t pdg = TMath::Abs(mctrack->
GetPdgCode());
1468 double momentum = mctrack->
GetP();
1472 Bool_t isAcc = (nMvdPoints + nStsPoints >= 4);
1481 if (vertex.Mag() < 0.1) {
1491 if (pdg == 111 && vertex.Mag() < 0.1) {
1495 if (pdg == 221 && vertex.Mag() < 0.1) {
1501 for (Int_t
i = 0;
i < ngTracks;
i++) {
1503 if (NULL == gTrack)
continue;
1509 if (stsInd < 0)
continue;
1511 if (stsTrack == NULL)
continue;
1514 if (stsMatch == NULL)
continue;
1517 if (stsMcTrackId < 0)
continue;
1519 if (mcTrack1 == NULL)
continue;
1520 int pdg = TMath::Abs(mcTrack1->
GetPdgCode());
1522 double momentum = mcTrack1->
GetP();
1529 if (richInd < 0 && trdInd < 0 && tofInd < 0) {
1532 if (richInd >= 0 && trdInd >= 0) {
1535 if (richInd >= 0 && trdInd >= 0 && tofInd >= 0) {
1539 if (vertex.Mag() < 0.1) {
1541 if (richInd < 0 && trdInd < 0 && tofInd < 0) {
1544 if (richInd >= 0 && trdInd >= 0) {
1547 if (richInd >= 0 && trdInd >= 0 && tofInd >= 0) {
1561 for (Int_t iGTrack = 0; iGTrack < ngTracks; iGTrack++) {
1565 if (NULL == gTrack)
continue;
1568 if (cand.
fStsInd < 0)
continue;
1570 if (stsTrack == NULL)
continue;
1590 Bool_t isRichRT = (cand.
fRichInd < 0) ?
false :
true;
1593 if (richRing == NULL) isRichRT =
false;
1600 Bool_t isTrdRT = (cand.
fTrdInd < 0) ?
false :
true;
1601 if (isTrdRT) isTrdRT =
fUseTrd;
1604 if (trdTrack == NULL) isTrdRT =
false;
1611 Bool_t isTofRT = (cand.
fTofInd < 0) ?
false :
true;
1614 if (tofHit == NULL) isTofRT =
false;
1620 if (isRichRT || isTrdRT || isTofRT) {
1626 cout <<
"fSTCandidates.size() = " <<
fSTCandidates.size() << endl;
1627 cout <<
"fRTCandidates.size() = " <<
fRTCandidates.size() << endl;
1639 for (Int_t iGTrack = 0; iGTrack < nGTracks; iGTrack++) {
1650 if (NULL == gTrack)
continue;
1653 if (cand.
fStsInd < 0)
continue;
1655 if (stsTrack == NULL)
continue;
1664 if (stsMatch == NULL)
continue;
1669 if (mcTrack1 == NULL)
continue;
1670 Int_t pdg = TMath::Abs(mcTrack1->
GetPdgCode());
1673 const FairTrackParam* richProjection =
1674 (FairTrackParam*) (
fRichProj->At(iGTrack));
1675 if (richProjection == NULL || richProjection->GetX() == 0
1676 || richProjection->GetY() == 0)
1695 if (cand.
fTrdInd < 0)
continue;
1697 if (trdTrack == NULL)
continue;
1702 if (cand.
fTofInd < 0)
continue;
1704 if (tofHit == NULL)
continue;
1713 cout <<
"fCandidates.size() = " <<
fCandidates.size() << endl;
1714 cout <<
"fTTCandidates.size() = " <<
fTTCandidates.size() << endl;
1722 for (
int i = 0;
i < nCand;
i++) {
1730 if (stsMatch == NULL)
continue;
1736 if (mcTrack1 == NULL)
continue;
1737 int pdg = TMath::Abs(mcTrack1->
GetPdgCode());
1748 if (richMatch == NULL)
continue;
1766 if (trdMatch == NULL)
continue;
1772 if (tofInd < 0)
continue;
1774 if (tofHit == NULL)
continue;
1776 if (tofHitMatch == NULL)
continue;
1778 if (tofPointIndex < 0)
continue;
1779 FairMCPoint* tofPoint = (FairMCPoint*)
fTofPoints->At(tofPointIndex);
1780 if (tofPoint == NULL)
continue;
1786 vector<CbmLmvmCandidate>& cutCandidates) {
1787 int nCand = cutCandidates.size();
1788 for (
int i = 0;
i < nCand;
i++) {
1789 cutCandidates[
i].ResetMcParams();
1791 int stsInd = cutCandidates[
i].fStsInd;
1792 if (stsInd < 0)
continue;
1795 if (stsMatch == NULL)
continue;
1798 cutCandidates[
i].fStsMcTrackId = stsMcTrackId;
1799 if (stsMcTrackId < 0)
continue;
1801 if (mcTrack1 == NULL)
continue;
1803 int pdg = TMath::Abs(mcTrack1->
GetPdgCode());
1805 cutCandidates[
i].fMcMotherId = motherId;
1806 cutCandidates[
i].fMcPdg = pdg;
1808 cutCandidates[
i].fIsMcSignalElectron =
1810 cutCandidates[
i].fIsMcPi0Electron =
1812 cutCandidates[
i].fIsMcGammaElectron =
1814 cutCandidates[
i].fIsMcEtaElectron =
1831 (!isEta) && (!isGamma) && (!isPi0)
1841 int binNum = (double) step + 0.5;
1848 hsp->Fill(0.5, 0.5);
1852 hsp->Fill(1.5, 0.5);
1856 hsp->Fill(2.5, 0.5);
1862 hsp->Fill(0.5, 1.5);
1867 hsp->Fill(1.5, 1.5);
1871 hsp->Fill(2.5, 1.5);
1877 hsp->Fill(0.5, 2.5);
1881 hsp->Fill(1.5, 2.5);
1885 hsp->Fill(2.5, 2.5);
1896 int binNum = (double) step + 0.5;
1923 if (NULL != mctrack) {
1930 sqrt(
v.X() *
v.X() +
v.Y() *
v.Y()));
1936 }
else if (pdg == 211 || pdg == -211) {
1938 }
else if (pdg == 2212) {
1940 }
else if (pdg == 321 || pdg == -321) {
1964 (!isEta) && (!isGamma) && (!isPi0)
1982 if (isBG && !isMismatch) {
2020 for (Int_t
i = 0;
i < ncand;
i++) {
2025 if (NULL != mcTrack) pdg = mcTrack->
GetPdgCode();
2037 if (!
fUseMvd) isMvd1Cut =
true;
2038 if (!
fUseMvd) isMvd2Cut =
true;
2043 if (isChi2Prim && isEl && isGammaCut)
2045 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut)
2047 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut && isMvd2Cut)
2049 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut)
2051 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut
2054 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut
2055 && isRtCut && isTtCut)
2057 if (isChi2Prim && isEl && isGammaCut && isMvd1Cut && isMvd2Cut && isStCut
2058 && isRtCut && isTtCut && isPtCut)
2062 for (Int_t iP = 0; iP < ncand; iP++) {
2067 for (Int_t iM = 0; iM < ncand; iM++) {
2072 if (iM == iP)
continue;
2075 if (mctrackP != NULL && mctrackM != NULL)
2098 Bool_t isMvd1Cut = (
fCandidates[iP].fIsMvd1CutElectron
2100 Bool_t isMvd2Cut = (
fCandidates[iP].fIsMvd2CutElectron
2102 if (!
fUseMvd) isMvd1Cut =
true;
2103 if (!
fUseMvd) isMvd2Cut =
true;
2110 if (isChiPrimary && isEl) {
2113 if (isChiPrimary && isEl && isGammaCut) {
2117 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut) {
2121 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut) {
2125 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut
2129 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut
2130 && isStCut && isRtCut) {
2133 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut
2134 && isStCut && isRtCut && isTtCut) {
2137 if (isChiPrimary && isEl && isGammaCut && isMvd1Cut && isMvd2Cut
2138 && isStCut && isRtCut && isTtCut && isPtCut) {
2147 for (Int_t iP = 0; iP < ncand; iP++) {
2149 for (Int_t iM = 0; iM < ncand; iM++) {
2151 if (iM == iP)
continue;
2169 const string& cutName,
2170 const vector<CbmLmvmCandidate>& cutCandidates,
2171 const vector<TH2D*>& hcut,
2172 const vector<TH2D*>& hcutPion,
2173 const vector<TH2D*>& hcutTruepair,
2176 vector<Double_t> angles, mom, candInd;
2178 Int_t nCutCand = cutCandidates.size();
2179 for (Int_t iP = 0; iP < nCand; iP++) {
2186 for (Int_t iM = 0; iM < nCutCand; iM++) {
2188 if (cutCandidates[iM].fCharge !=
fCandidates[iP].fCharge) {
2192 angles.push_back(pRec.
fAngle);
2193 mom.push_back(cutCandidates[iM].fMomentum.Mag());
2194 candInd.push_back(iM);
2198 Double_t minAng = 360.;
2200 for (UInt_t
i = 0;
i < angles.size();
i++) {
2201 if (minAng > angles[
i]) {
2207 if (cutName ==
"TT")
fCandidates[iP].fIsTtCutElectron =
true;
2208 if (cutName ==
"ST")
fCandidates[iP].fIsStCutElectron =
true;
2209 if (cutName ==
"RT")
fCandidates[iP].fIsRtCutElectron =
true;
2213 TMath::Sqrt(
fCandidates[iP].fMomentum.Mag() * mom[minInd]);
2214 Double_t val = -1. * (angleCut / ppCut) * sqrt_mom + angleCut;
2215 if (!(sqrt_mom < ppCut && val > minAng)) {
2216 if (cutName ==
"TT")
fCandidates[iP].fIsTtCutElectron =
true;
2217 if (cutName ==
"ST")
fCandidates[iP].fIsStCutElectron =
true;
2218 if (cutName ==
"RT")
fCandidates[iP].fIsRtCutElectron =
true;
2221 int stsInd = cutCandidates[candInd[minInd]].fStsInd;
2222 if (stsInd < 0)
continue;
2223 int pdg = TMath::Abs(cutCandidates[candInd[minInd]].fMcPdg);
2224 int motherId = cutCandidates[candInd[minInd]].fMcMotherId;
2228 if (pdg == 211) hcutPion[
kSignal]->Fill(sqrt_mom, minAng,
fWeight);
2232 hcut[
kBg]->Fill(sqrt_mom, minAng);
2233 if (pdg == 211) hcutPion[
kBg]->Fill(sqrt_mom, minAng);
2235 if (motherId != -1 && motherId ==
fCandidates[iP].fMcMotherId)
2236 hcutTruepair[
kBg]->Fill(sqrt_mom, minAng);
2240 hcut[
kPi0]->Fill(sqrt_mom, minAng);
2241 if (pdg == 211) hcutPion[
kPi0]->Fill(sqrt_mom, minAng);
2243 if (motherId != -1 && motherId ==
fCandidates[iP].fMcMotherId)
2244 hcutTruepair[
kPi0]->Fill(sqrt_mom, minAng);
2248 hcut[
kGamma]->Fill(sqrt_mom, minAng);
2249 if (pdg == 211) hcutPion[
kGamma]->Fill(sqrt_mom, minAng);
2251 if (motherId != -1 && motherId ==
fCandidates[iP].fMcMotherId)
2252 hcutTruepair[
kGamma]->Fill(sqrt_mom, minAng);
2267 const string& source) {
2269 vector<bool> isAdded;
2270 isAdded.resize(nCand);
2271 for (UInt_t iP = 0; iP < nCand; iP++) {
2272 isAdded[iP] =
false;
2274 for (UInt_t iP = 0; iP < nCand; iP++) {
2276 if (source ==
"pi0" && !
fCandidates[iP].fIsMcPi0Electron)
continue;
2277 if (source ==
"gamma" && !
fCandidates[iP].fIsMcGammaElectron)
continue;
2283 if (isAdded[iP])
continue;
2287 h_nof_pairs->Fill(4.5);
2292 if (isAdded[iP])
continue;
2296 h_nof_pairs->Fill(5.5);
2301 if (isAdded[iP])
continue;
2305 h_nof_pairs->Fill(6.5);
2310 if (isAdded[iP])
continue;
2312 for (UInt_t iM = 0; iM <
fCandidates.size(); iM++) {
2316 h_nof_pairs->Fill(7.5);
2323 if (isAdded[iP])
continue;
2324 Int_t nofPointsSts = 0;
2325 Int_t nofMcTracks =
fMCTracks->GetEntriesFast();
2326 for (Int_t iMCTrack = 0; iMCTrack < nofMcTracks; iMCTrack++) {
2330 if (iMCTrack ==
fCandidates[iP].fStsMcTrackId)
continue;
2333 FairRun::Instance()->GetEventHeader()->GetMCEntryNumber(),
2338 FairRun::Instance()->GetEventHeader()->GetMCEntryNumber(), iMCTrack);
2343 if (nofPointsSts == 0) h_nof_pairs->Fill(0.5);
2344 if (nofPointsSts >= 1 && nofPointsSts <= 3) h_nof_pairs->Fill(1.5);
2345 if (nofPointsSts >= 4 && nofPointsSts <= 5) h_nof_pairs->Fill(2.5);
2346 if (nofPointsSts >= 6) h_nof_pairs->Fill(3.5);
2372 globalTrackIndex, momentum);
2374 globalTrackIndex, momentum);
2377 globalTrackIndex, momentum);
2379 globalTrackIndex, momentum);
2382 globalTrackIndex, momentum);
2387 if (richEl && trdEl && tofEl && momCut) {
2400 if (pdg == 11 || pdg == -11) {
2418 for (Int_t
i = 0;
i < nCand;
i++) {
2462 for (Int_t
i = 0;
i < nCand;
i++) {
2490 for (
int i = 0;
i < nCand;
i++) {
2497 if (NULL == track)
continue;
2500 double mvd1x = 0., mvd1y = 0., mvd2x = 0., mvd2y = 0.;
2502 for (Int_t ith = 0; ith < nofMvdHits; ith++) {
2506 if (NULL == pmh)
continue;
2507 if (stationNum == 1) {
2508 mvd1x = pmh->
GetX();
2509 mvd1y = pmh->
GetY();
2510 }
else if (stationNum == 2) {
2511 mvd2x = pmh->
GetX();
2512 mvd2y = pmh->
GetY();
2515 double mvd1r =
sqrt(mvd1x * mvd1x + mvd1y * mvd1y);
2516 double mvd2r =
sqrt(mvd2x * mvd2x + mvd2y * mvd2y);
2554 vector<TH2D*>& hist,
2555 vector<TH1D*>& histQa) {
2559 vector<float> candX;
2560 vector<float> candY;
2561 vector<int> candInd;
2565 Int_t nMvdHit =
fMvdHits->GetEntriesFast();
2566 for (Int_t iHit = 0; iHit < nMvdHit; iHit++) {
2568 if (NULL == pmh)
continue;
2570 if (stationNum == mvdStationNum) {
2571 mvdX.push_back(pmh->
GetX());
2572 mvdY.push_back(pmh->
GetY());
2573 mvdInd.push_back(iHit);
2578 for (Int_t
i = 0;
i < nCand;
i++) {
2583 if (NULL == track)
continue;
2585 for (Int_t ith = 0; ith < nhits; ith++) {
2589 if (NULL == pmh)
continue;
2590 if (stationNum == mvdStationNum) {
2591 candX.push_back(pmh->
GetX());
2592 candY.push_back(pmh->
GetY());
2593 candInd.push_back(
i);
2599 for (UInt_t iT = 0; iT < candInd.size(); iT++) {
2600 Float_t mind = 9999999.;
2601 Int_t minMvdInd = -1;
2602 for (UInt_t iH = 0; iH < mvdX.size(); iH++) {
2603 Float_t dx = mvdX[iH] - candX[iT];
2604 Float_t dy = mvdY[iH] - candY[iT];
2605 Float_t d2 = dx * dx + dy * dy;
2606 if (d2 < 1.e-9)
continue;
2608 minMvdInd = mvdInd[iH];
2612 Double_t dmvd =
sqrt(mind);
2618 if (NULL != hitMatch) {
2621 int mcMvdHitPdg = TMath::Abs(mct1->
GetPdgCode());
2624 int stsMcTrackId =
fCandidates[candInd[iT]].fStsMcTrackId;
2625 int stsMotherId = -2;
2626 if (stsMcTrackId >= 0) {
2632 if (mvdMotherId != -1 && mvdMotherId == stsMotherId) {
2638 if (
fCandidates[candInd[iT]].fIsMcSignalElectron) {
2639 if (mvdMotherId == stsMotherId && mcMvdHitPdg == 11) {
2650 if (
fCandidates[candInd[iT]].fIsMcSignalElectron) {
2656 histQa[
kBg]->Fill(bin);
2658 if (
fCandidates[candInd[iT]].fIsMcGammaElectron) {
2660 histQa[
kGamma]->Fill(bin);
2664 histQa[
kPi0]->Fill(bin);
2668 if (mvdStationNum == 1) {
2669 Double_t mom =
fCandidates[candInd[iT]].fMomentum.Mag();
2672 if (!(dmvd < fCuts.fMvd1CutD && val > mom)) {
2673 fCandidates[candInd[iT]].fIsMvd1CutElectron =
true;
2675 fCandidates[candInd[iT]].fIsMvd1CutElectron =
false;
2678 if (mvdStationNum == 2) {
2679 Double_t mom =
fCandidates[candInd[iT]].fMomentum.Mag();
2682 if (!(dmvd < fCuts.fMvd2CutD && val > mom)) {
2683 fCandidates[candInd[iT]].fIsMvd2CutElectron =
true;
2685 fCandidates[candInd[iT]].fIsMvd2CutElectron =
false;
2694 for (Int_t
i = 0;
i < nCand;
i++) {
2699 if (NULL == track)
continue;
2702 for (Int_t ith = 0; ith < nhits; ith++) {
2705 if (NULL == pmh1)
continue;
2709 for (
int iMvd = 0; iMvd < nofMvdHits; iMvd++) {
2712 if (NULL == hitMatch)
continue;
2714 if (stsMcTrackId != mcMvdHitId)
continue;
2717 double dx = pmh1->
GetX() - pmh2->
GetX();
2718 double dy = pmh1->
GetY() - pmh2->
GetY();
2719 double d =
sqrt(dx * dx + dy * dy);
2720 if (stationNum == 1) {
2725 }
else if (stationNum == 1) {
2738 TDirectory* oldir = gDirectory;
2739 TFile* outFile = FairRootManager::Instance()->GetOutFile();
2740 if (outFile != NULL) {
2747 gDirectory->cd(oldir->GetPath());
2751 const string& particle) {
2770 if (energy ==
"8gev") {
2772 if (particle ==
"omegaepem") this->
SetWeight(2.5 * 7.28e-5);
2774 if (particle ==
"omegadalitz") this->
SetWeight(2.5 * 7.7e-4);
2776 if (particle ==
"phi") this->
SetWeight(0.365 * 2.97e-4);
2778 if (particle ==
"inmed") this->
SetWeight(0.5 * 4.45e-2);
2780 if (particle ==
"qgp") this->
SetWeight(0.5 * 1.15e-2);
2781 }
else if (energy ==
"25gev") {
2783 if (particle ==
"rho0") this->
SetWeight(23 * 4.7e-5);
2785 if (particle ==
"omegaepem") this->
SetWeight(38 * 7.28e-5);
2787 if (particle ==
"omegadalitz") this->
SetWeight(38 * 7.7e-4);
2789 if (particle ==
"phi") this->
SetWeight(1.28 * 2.97e-4);
2791 if (particle ==
"inmed") this->
SetWeight(4.45e-2);
2793 if (particle ==
"qgp") this->
SetWeight(1.15e-2);
2794 }
else if (energy ==
"3.5gev") {
2796 if (particle ==
"rho0") this->
SetWeight(1.0 * 4.7e-5);
2798 if (particle ==
"omegaepem") this->
SetWeight(1.2 * 7.28e-5);
2800 if (particle ==
"omegadalitz") this->
SetWeight(1.2 * 7.7e-5);
2802 if (particle ==
"phi") this->
SetWeight(0.1 * 2.97e-4);
2816 }
else if (energy ==
"4.5gev") {
2818 if (particle ==
"omegadalitz") this->
SetWeight(1.2 * 7.7e-5);
2820 if (particle ==
"omegaepem") this->
SetWeight(1.2 * 7.28e-6);
2822 if (particle ==
"phi") this->
SetWeight(1.2 * 2.97e-6);
2824 if (particle ==
"inmed") this->
SetWeight(2.4 * 10e-3);
2826 cout <<
"-ERROR- CbmAnaDielectronTask::SetEnergyAndParticle energy or "
2827 "particle is not correct, energy:"
2828 << energy <<
" particle:" << particle << endl;