CbmRoot
Loading...
Searching...
No Matches
CbmKresEtaMCAnalysis.cxx
Go to the documentation of this file.
1/* Copyright (C) 2018-2020 GSI Helmholtzzentrum fuer Schwerionenforschung, Darmstadt
2 SPDX-License-Identifier: GPL-3.0-only
3 Authors: Ievgenii Kres, Florian Uhlig [committer] */
4
21
22#include "CbmGlobalTrack.h"
23#include "CbmKresFunctions.h"
24#include "CbmMCTrack.h"
25#include "CbmStsTrack.h"
26#include "CbmTrackMatchNew.h"
27
28#include "FairRootManager.h"
29
30#include "TDirectory.h"
31#include "TH1.h"
32#include "TH1D.h"
33#include "TH2D.h"
34#include "TMath.h"
35
36#include <iostream>
37
38
39using namespace std;
40
42 : fMcTracks(nullptr)
43 , fGlobalTracks(nullptr)
44 , fStsTracks(nullptr)
45 , fStsTrackMatches(nullptr)
46 , El_Photon_Eta_refmomentum()
47 , El_Photon_Eta_MCtrack()
48 , El_Photon_Eta_Id()
49 , El_Photon_Pion_Eta_refmomentum()
50 , El_Photon_Pion_Eta_MCtrack()
51 , El_Photon_Pion_Eta_Id()
52 , Pion_Eta_refmomentum()
53 , Pion_Eta_MCtrack()
54 , Pion_Eta_Id()
55 , All_El_refmomentum()
56 , All_El_MCtrack()
57 , All_El_Id()
58 , All_Pion_refmomentum()
59 , All_Pion_MCtrack()
60 , All_Pion_Id()
61 , frefmomenta()
62 , fMCtracks()
63 , fMCId()
64 , EDGA_RefMom()
65 , EDGA_MC()
66 , EDGA_Id()
67 , ECPGA_leptons_RefMom()
68 , ECPGA_leptons_MC()
69 , ECPGA_leptons_Id()
70 , ECPGA_pions_RefMom()
71 , ECPGA_pions_MC()
72 , fHistoList_eta_gg()
73 , InvMass_eta_gg_mc(nullptr)
74 , InvMass_eta_gg_reffited(nullptr)
75 , InvMassPhoton_eta_gg_mc(nullptr)
76 , InvMassPhoton_eta_gg_reffited(nullptr)
77 , OpeningAnglePhoton_eta_gg_mc(nullptr)
78 , OpeningAnglePhoton_eta_gg_reffited(nullptr)
79 , OpeningAngle_eta_gg_between_gg_mc(nullptr)
80 , OpeningAngle_eta_gg_between_gg_reffited(nullptr)
81 , InvMass_eta_gg_allcombinations_mc(nullptr)
82 , InvMass_eta_gg_allcombinations_reffited(nullptr)
83 , EMT_eta_gg(nullptr)
84 , InvMass_eta_gg_reco_aftercuts(nullptr)
85 , rap_vs_pt_eta_gg_reco_aftercuts(nullptr)
86 , rap_vs_pt_NOTeta_gg_reco_aftercuts(nullptr)
87 , fHistoList_eta_ppg()
88 , InvMass_eta_ppg_mc(nullptr)
89 , InvMass_eta_ppg_reffited(nullptr)
90 , InvMassPhoton_eta_ppg_mc(nullptr)
91 , InvMassPhoton_eta_ppg_reffited(nullptr)
92 , OpeningAnglePhoton_eta_ppg_mc(nullptr)
93 , OpeningAnglePhoton_eta_ppg_reffited(nullptr)
94 , InvMass_eta_ppg_allcombinations_mc(nullptr)
95 , InvMass_eta_ppg_allcombinations_reffited(nullptr)
96 , Pion_P_fromEta_reco(nullptr)
97 , Pion_P_elsewhere_reco(nullptr)
98 , Pion_Pt_fromEta_reco(nullptr)
99 , Pion_Pt_elsewhere_reco(nullptr)
100 , OA_betweenPions_fromEta_mc(nullptr)
101 , OA_betweenPions_fromEta_reco(nullptr)
102 , OA_betweenPions_fromEta_reco_wrongcombinations(nullptr)
103 , EMT_eta_ppg(nullptr)
104 , EMT_eta_three_body(nullptr)
105 , InvMass_eta_ppg_reco_aftercuts(nullptr)
106 , rap_vs_pt_eta_ppg_reco_aftercuts(nullptr)
107 , rap_vs_pt_NOTeta_ppg_reco_aftercuts(nullptr)
108 , fHistoList_eta_ppp()
109 , InvMass_eta_ppp_mc(nullptr)
110 , InvMass_eta_ppp_reffited(nullptr)
111 , InvMass_eta_Npion_mc(nullptr)
112 , InvMass_eta_Npion_reffited(nullptr)
113 , InvMass_eta_ppp_allcombinations_mc(nullptr)
114 , InvMass_eta_ppp_allcombinations_reffited(nullptr)
115 , EMT_gg_Event()
116 , EMT_gg_pair_momenta()
117 , EMT_ppg_ee_Event()
118 , EMT_ppg_ee_pair_momenta()
119 , EMT_ppg_pp_Event()
120 , EMT_ppg_pp_pair_momenta()
121 , EMT_ppg_positive_pion_Event()
122 , EMT_ppg_positive_pion_momenta()
123 , EMT_ppg_negative_pion_Event()
124 , EMT_ppg_negative_pion_momenta()
125{
126}
127
129
131{
133
134 FairRootManager* ioman = FairRootManager::Instance();
135 if (nullptr == ioman) { Fatal("CbmKresEtaMCAnalysis::Init", "RootManager not instantised!"); }
136
137 fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
138 if (nullptr == fMcTracks) { Fatal("CbmKresEtaMCAnalysis::Init", "No MCTrack array!"); }
139
140 fGlobalTracks = (TClonesArray*) ioman->GetObject("GlobalTrack");
141 if (nullptr == fGlobalTracks) { Fatal("CbmKresEtaMCAnalysis::Init", "No GlobalTrack array!"); }
142
143 fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
144 if (nullptr == fStsTracks) { Fatal("CbmKresEtaMCAnalysis::Init", "No StsTrack array!"); }
145
146 fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
147 if (nullptr == fStsTrackMatches) { Fatal("CbmKresEtaMCAnalysis::Init", "No StsTrackMatch array!"); }
148}
149
150
151void CbmKresEtaMCAnalysis::Exec(int Event, double OpeningAngleCut, double GammaInvMassCut)
152{
153 // cout << "CbmKresEtaMCAnalysis, event No. " << Event << endl;
154
155 int counter = 0;
156 Int_t nofMcTracks = fMcTracks->GetEntriesFast();
157 for (int i = 0; i < nofMcTracks; i++) {
158 CbmMCTrack* mctrack = (CbmMCTrack*) fMcTracks->At(i);
159 if (mctrack == nullptr) continue;
160 if (mctrack->GetPdgCode() == 221) counter++;
161
162 // if (mctrack->GetMotherId() == -1) continue;
163 // CbmMCTrack* mcMotherTrack_1 = (CbmMCTrack*) fMcTracks->At(mctrack->GetMotherId());
164 // if (mcMotherTrack_1 == nullptr) continue;
165 // if (TMath::Abs(mctrack->GetPdgCode()) == 11 && mcMotherTrack_1->GetPdgCode() == 22){
166 // if(mcMotherTrack_1->GetMotherId() == -1) continue;
167 // CbmMCTrack* mcGrTrack = (CbmMCTrack*) fMcTracks->At(mcMotherTrack_1->GetMotherId());
168 // if (mcGrTrack == nullptr) continue;
169 // if (mcGrTrack->GetPdgCode() == 221){
170 // cout<< "Electron: id = " << i << "; motherid = " << mctrack->GetMotherId() << "; GrId = " << mcMotherTrack_1->GetMotherId() << endl;
171 // }
172 // }
173 }
174 // cout << "number of etas in event = " << counter << endl;
175
177 El_Photon_Eta_MCtrack.clear();
178 El_Photon_Eta_Id.clear();
179
182 El_Photon_Pion_Eta_Id.clear();
183
184 Pion_Eta_refmomentum.clear();
185 Pion_Eta_MCtrack.clear();
186 Pion_Eta_Id.clear();
187
188 All_El_refmomentum.clear();
189 All_El_MCtrack.clear();
190 All_El_Id.clear();
191
192 All_Pion_refmomentum.clear();
193 All_Pion_MCtrack.clear();
194 All_Pion_Id.clear();
195
196
198 // ========================================================================================
199 Int_t ngTracks = fGlobalTracks->GetEntriesFast();
200 for (Int_t iTr = 0; iTr < ngTracks; iTr++) {
201 CbmGlobalTrack* gTrack = (CbmGlobalTrack*) fGlobalTracks->At(iTr);
202 if (nullptr == gTrack) continue;
203 int stsInd = gTrack->GetStsTrackIndex();
204
205 // ========================================================================================
207 if (stsInd < 0) continue;
208 CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(stsInd);
209 if (stsTrack == nullptr) continue;
210 CbmTrackMatchNew* stsMatch = (CbmTrackMatchNew*) fStsTrackMatches->At(stsInd);
211 if (stsMatch == nullptr) continue;
212 if (stsMatch->GetNofLinks() <= 0) continue;
213 int McTrackId = stsMatch->GetMatchedLink().GetIndex();
214 if (McTrackId < 0) continue;
215 CbmMCTrack* mcTrack = (CbmMCTrack*) fMcTracks->At(McTrackId);
216 if (mcTrack == nullptr) continue;
217 int pdg = mcTrack->GetPdgCode();
218 int motherId = mcTrack->GetMotherId();
219 if (TMath::Abs(pdg) != 11 && TMath::Abs(pdg) != 211) continue;
221 // ========================================================================================
222
223 double chi2 = 0;
224 TVector3 refmomentum = CbmKresFunctions::FitToVertexAndGetChi(stsTrack, mcTrack->GetStartX(), mcTrack->GetStartY(),
225 mcTrack->GetStartZ(), chi2);
226
227 // cout << "" << endl;
228 // cout << " pdg = " << mcTrack->GetPdgCode() << endl;
229 // cout << "mc: P = " << mcTrack->GetP() << ";\t Pt = " << mcTrack->GetPt() << "; \t E = " << mcTrack->GetEnergy() << endl;
230 // cout << "reco: P = " << refmomentum.Mag() << ";\t Pt = " << refmomentum.Perp() << "; \t E = " << TMath::Sqrt(refmomentum.Mag2() + mcTrack->GetMass()*mcTrack->GetMass()) << endl;
231 // cout << "chi2 = " << chi2 << endl;
232 // cout << "=========================" << endl;
233
234 // if (chi2 > 3) continue;
235
236 if (TMath::Abs(pdg) == 11) {
237 All_El_refmomentum.push_back(refmomentum);
238 All_El_MCtrack.push_back(mcTrack);
239 All_El_Id.push_back(McTrackId);
240 }
241 if (TMath::Abs(pdg) == 211) {
242 All_Pion_refmomentum.push_back(refmomentum);
243 All_Pion_MCtrack.push_back(mcTrack);
244 All_Pion_Id.push_back(McTrackId);
245 Pion_P_elsewhere_reco->Fill(refmomentum.Mag());
246 Pion_Pt_elsewhere_reco->Fill(refmomentum.Perp());
247 }
248
249 if (motherId == -1) continue;
250 CbmMCTrack* mcMotherTrack = (CbmMCTrack*) fMcTracks->At(motherId);
251 if (mcMotherTrack == nullptr) continue;
252
253 if (TMath::Abs(pdg) == 11 && mcMotherTrack->GetPdgCode() == 22) {
254 if (mcMotherTrack->GetMotherId() == -1) continue;
255 CbmMCTrack* mcGrTrack = (CbmMCTrack*) fMcTracks->At(mcMotherTrack->GetMotherId());
256 if (mcGrTrack == nullptr) continue;
257 if (mcGrTrack->GetPdgCode() == 221) {
258 El_Photon_Eta_refmomentum.push_back(refmomentum);
259 El_Photon_Eta_MCtrack.push_back(mcTrack);
260 El_Photon_Eta_Id.push_back(McTrackId);
261 }
262 if (mcGrTrack->GetPdgCode() == 111) {
263 if (mcGrTrack->GetMotherId() == -1) continue;
264 CbmMCTrack* mcTopTrack = (CbmMCTrack*) fMcTracks->At(mcGrTrack->GetMotherId());
265 if (mcTopTrack->GetPdgCode() == 221) {
266 El_Photon_Pion_Eta_refmomentum.push_back(refmomentum);
267 El_Photon_Pion_Eta_MCtrack.push_back(mcTrack);
268 El_Photon_Pion_Eta_Id.push_back(McTrackId);
269 }
270 }
271 }
272
273 if (TMath::Abs(pdg) == 211 && mcMotherTrack->GetPdgCode() == 221) {
274 Pion_Eta_refmomentum.push_back(refmomentum);
275 Pion_Eta_MCtrack.push_back(mcTrack);
276 Pion_Eta_Id.push_back(McTrackId);
277 Pion_P_fromEta_reco->Fill(refmomentum.Mag());
278 Pion_Pt_fromEta_reco->Fill(refmomentum.Perp());
279 }
280 }
282 // ========================================================================================
283
284
290
291
292 EtaDoubleGammaAnalysis_plusBG(OpeningAngleCut, GammaInvMassCut, Event, All_El_refmomentum, All_El_MCtrack, All_El_Id,
294 EtaChargedPionsGammaAnalysis_plusBG(OpeningAngleCut, GammaInvMassCut, Event, All_Pion_refmomentum, All_Pion_MCtrack,
297
298 if (Event % 1000 == 0) {
299 // Mixing_gg();
300 EMT_gg_Event.clear();
301 EMT_gg_pair_momenta.clear();
302 }
303
304 if (Event % 10 == 0) {
305 // Mixing_three_body();
310
311 // Mixing_ppg();
312 EMT_ppg_ee_Event.clear();
314 EMT_ppg_pp_Event.clear();
316 }
317}
318
320 CbmMCTrack* mctrack3, CbmMCTrack* mctrack4)
321{
322 Double_t openingAngle;
323 TLorentzVector gamma1;
324 TLorentzVector gamma2;
325
326 if (mctrack1->GetMotherId() == mctrack2->GetMotherId() && mctrack3->GetMotherId() == mctrack4->GetMotherId()) {
327 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
328 mother1->Get4Momentum(gamma1);
329 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack3->GetMotherId());
330 mother2->Get4Momentum(gamma2);
331 }
332 if (mctrack1->GetMotherId() == mctrack3->GetMotherId() && mctrack2->GetMotherId() == mctrack4->GetMotherId()) {
333 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
334 mother1->Get4Momentum(gamma1);
335 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack2->GetMotherId());
336 mother2->Get4Momentum(gamma2);
337 }
338 if (mctrack1->GetMotherId() == mctrack4->GetMotherId() && mctrack2->GetMotherId() == mctrack3->GetMotherId()) {
339 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(mctrack1->GetMotherId());
340 mother1->Get4Momentum(gamma1);
341 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(mctrack2->GetMotherId());
342 mother2->Get4Momentum(gamma2);
343 }
344
345 Double_t angle = gamma1.Angle(gamma2.Vect());
346 openingAngle = 180. * angle / TMath::Pi();
347
348 return openingAngle;
349}
350
351Double_t CbmKresEtaMCAnalysis::CalculateOpeningAngleBetweenGammas_Reco(TVector3 electron1, TVector3 electron2,
352 TVector3 electron3, TVector3 electron4)
353{
354 double M2El = 2.6112004954086e-7;
355 Double_t energy1 = TMath::Sqrt(electron1.Mag2() + M2El);
356 TLorentzVector lorVec1(electron1, energy1);
357
358 Double_t energy2 = TMath::Sqrt(electron2.Mag2() + M2El);
359 TLorentzVector lorVec2(electron2, energy2);
360
361 Double_t energy3 = TMath::Sqrt(electron3.Mag2() + M2El);
362 TLorentzVector lorVec3(electron3, energy3);
363
364 Double_t energy4 = TMath::Sqrt(electron4.Mag2() + M2El);
365 TLorentzVector lorVec4(electron4, energy4);
366
367 TLorentzVector lorPhoton1 = lorVec1 + lorVec2;
368 TLorentzVector lorPhoton2 = lorVec3 + lorVec4;
369
370 Double_t angleBetweenPhotons = lorPhoton1.Angle(lorPhoton2.Vect());
371 Double_t theta = 180. * angleBetweenPhotons / TMath::Pi();
372
373 return theta;
374}
375
376
377void CbmKresEtaMCAnalysis::EtaDoubleGammaAnalysis(vector<TVector3> RefMom, vector<CbmMCTrack*> MC, vector<Int_t> Id,
378 vector<TH1*> gg)
379{
380 //================================== decay eta -> gamma gamma -> e+e- e+e-
381 // cout << "MC.size() = " << MC.size() << endl;
382 if (MC.size() < 4) return;
383
384 for (size_t i = 0; i < MC.size(); i++) {
385 for (size_t j = i + 1; j < MC.size(); j++) {
386 for (size_t k = j + 1; k < MC.size(); k++) {
387 for (size_t l = k + 1; l < MC.size(); l++) {
388
389 int pdg1 = MC.at(i)->GetPdgCode();
390 int pdg2 = MC.at(j)->GetPdgCode();
391 int pdg3 = MC.at(k)->GetPdgCode();
392 int pdg4 = MC.at(l)->GetPdgCode();
393
394 if (pdg1 + pdg2 + pdg3 + pdg4 != 0) continue;
395 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11 || TMath::Abs(pdg3) != 11 || TMath::Abs(pdg4) != 11)
396 continue;
397
398 int motherId1 = MC.at(i)->GetMotherId();
399 int motherId2 = MC.at(j)->GetMotherId();
400 int motherId3 = MC.at(k)->GetMotherId();
401 int motherId4 = MC.at(l)->GetMotherId();
402
403 int mcId1 = Id.at(i);
404 int mcId2 = Id.at(j);
405 int mcId3 = Id.at(k);
406 int mcId4 = Id.at(l);
407
408 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
409 continue; // particle is not used twice
410 if (motherId1 == -1 || motherId2 == -1 || motherId3 == -1 || motherId4 == -1) continue;
411
412 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
413 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
414 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
415 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
416
417 int mcMotherPdg1 = mother1->GetPdgCode();
418 int mcMotherPdg2 = mother2->GetPdgCode();
419 int mcMotherPdg3 = mother3->GetPdgCode();
420 int mcMotherPdg4 = mother4->GetPdgCode();
421
422 if (mcMotherPdg1 != 22 || mcMotherPdg2 != 22 || mcMotherPdg3 != 22 || mcMotherPdg4 != 22) continue;
423
424 int grandmotherId1 = mother1->GetMotherId();
425 int grandmotherId2 = mother2->GetMotherId();
426 int grandmotherId3 = mother3->GetMotherId();
427 int grandmotherId4 = mother4->GetMotherId();
428
429 if (grandmotherId1 == -1) continue;
430
431
432 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
433 && grandmotherId1 == grandmotherId4) {
434
435 CbmMCTrack* GrMother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
436 int mcGrandmotherPdg1 = GrMother1->GetPdgCode();
437 if (mcGrandmotherPdg1 != 221) continue;
438
439 Double_t InvMass_true = CbmKresFunctions::Invmass_4particles_MC(MC.at(i), MC.at(j), MC.at(k), MC.at(l));
440 Double_t InvMass_reco =
441 CbmKresFunctions::Invmass_4particles_RECO(RefMom.at(i), RefMom.at(j), RefMom.at(k), RefMom.at(l));
442 cout << "Decay eta -> gamma gamma -> e+e- e+e- detected!\t\t mc mass: " << InvMass_true
443 << "\t, reco mass: " << InvMass_reco << endl;
444 cout << "motherids: " << motherId1 << "/" << motherId2 << "/" << motherId3 << "/" << motherId4 << endl;
445 cout << "grandmotherid: " << grandmotherId1 << "/" << grandmotherId2 << "/" << grandmotherId3 << "/"
446 << grandmotherId4 << endl;
447 cout << "pdgs: " << mcGrandmotherPdg1 << "-->" << mcMotherPdg1 << "/" << mcMotherPdg2 << "/" << mcMotherPdg3
448 << "/" << mcMotherPdg4 << "-->" << pdg1 << "/" << pdg2 << "/" << pdg3 << "/" << pdg4 << endl;
449
450 gg[0]->Fill(InvMass_true);
451 gg[1]->Fill(InvMass_reco);
452
453 cout << "\t \t mc true info: " << endl;
454 cout << "particle 1: \t" << MC.at(i)->GetPdgCode() << ";\t pt = " << MC.at(i)->GetPt()
455 << ";\t X = " << MC.at(i)->GetStartX() << ";\t Y = " << MC.at(i)->GetStartY()
456 << ";\t Z = " << MC.at(i)->GetStartZ() << ";\t E = " << MC.at(i)->GetEnergy() << endl;
457 cout << "particle 2: \t" << MC.at(j)->GetPdgCode() << ";\t pt = " << MC.at(j)->GetPt()
458 << ";\t X = " << MC.at(j)->GetStartX() << ";\t Y = " << MC.at(j)->GetStartY()
459 << ";\t Z = " << MC.at(j)->GetStartZ() << ";\t E = " << MC.at(j)->GetEnergy() << endl;
460 cout << "particle 3: \t" << MC.at(k)->GetPdgCode() << ";\t pt = " << MC.at(k)->GetPt()
461 << ";\t X = " << MC.at(k)->GetStartX() << ";\t Y = " << MC.at(k)->GetStartY()
462 << ";\t Z = " << MC.at(k)->GetStartZ() << ";\t E = " << MC.at(k)->GetEnergy() << endl;
463 cout << "particle 4: \t" << MC.at(l)->GetPdgCode() << ";\t pt = " << MC.at(l)->GetPt()
464 << ";\t X = " << MC.at(l)->GetStartX() << ";\t Y = " << MC.at(l)->GetStartY()
465 << ";\t Z = " << MC.at(l)->GetStartZ() << ";\t E = " << MC.at(l)->GetEnergy() << endl;
466
467 Double_t OpeningAngle1_mc = 0;
468 Double_t OpeningAngle2_mc = 0;
469 Double_t OpeningAngle1_refitted = 0;
470 Double_t OpeningAngle2_refitted = 0;
471 Double_t InvMass_realg1_mc = 0;
472 Double_t InvMass_realg2_mc = 0;
473 Double_t InvMass_realg1_refitted = 0;
474 Double_t InvMass_realg2_refitted = 0;
475
476 Double_t openingAngleBetweenGammasReco = 0;
477
478 if (motherId1 == motherId2) {
479 OpeningAngle1_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(i), MC.at(j));
480 OpeningAngle2_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(k), MC.at(l));
481 gg[4]->Fill(OpeningAngle1_mc);
482 gg[4]->Fill(OpeningAngle2_mc);
483 OpeningAngle1_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(i), RefMom.at(j));
484 OpeningAngle2_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(k), RefMom.at(l));
485 gg[5]->Fill(OpeningAngle1_refitted);
486 gg[5]->Fill(OpeningAngle2_refitted);
487
488 InvMass_realg1_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(i), MC.at(j));
489 InvMass_realg2_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(k), MC.at(l));
490 gg[2]->Fill(InvMass_realg1_mc);
491 gg[2]->Fill(InvMass_realg2_mc);
492 InvMass_realg1_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(i), RefMom.at(j));
493 InvMass_realg2_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(k), RefMom.at(l));
494 gg[3]->Fill(InvMass_realg1_refitted);
495 gg[3]->Fill(InvMass_realg2_refitted);
496 openingAngleBetweenGammasReco =
497 CalculateOpeningAngleBetweenGammas_Reco(RefMom.at(i), RefMom.at(j), RefMom.at(k), RefMom.at(l));
498 }
499 if (motherId1 == motherId3) {
500 OpeningAngle1_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(i), MC.at(k));
501 OpeningAngle2_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(j), MC.at(l));
502 gg[4]->Fill(OpeningAngle1_mc);
503 gg[4]->Fill(OpeningAngle2_mc);
504 OpeningAngle1_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(i), RefMom.at(k));
505 OpeningAngle2_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(j), RefMom.at(l));
506 gg[5]->Fill(OpeningAngle1_refitted);
507 gg[5]->Fill(OpeningAngle2_refitted);
508
509 InvMass_realg1_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(i), MC.at(k));
510 InvMass_realg2_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(j), MC.at(l));
511 gg[2]->Fill(InvMass_realg1_mc);
512 gg[2]->Fill(InvMass_realg2_mc);
513 InvMass_realg1_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(i), RefMom.at(k));
514 InvMass_realg2_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(j), RefMom.at(l));
515 gg[3]->Fill(InvMass_realg1_refitted);
516 gg[3]->Fill(InvMass_realg2_refitted);
517 openingAngleBetweenGammasReco =
518 CalculateOpeningAngleBetweenGammas_Reco(RefMom.at(i), RefMom.at(k), RefMom.at(j), RefMom.at(l));
519 }
520 if (motherId1 == motherId4) {
521 OpeningAngle1_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(i), MC.at(l));
522 OpeningAngle2_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MC.at(j), MC.at(k));
523 gg[4]->Fill(OpeningAngle1_mc);
524 gg[4]->Fill(OpeningAngle2_mc);
525 OpeningAngle1_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(i), RefMom.at(l));
526 OpeningAngle2_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(j), RefMom.at(k));
527 gg[5]->Fill(OpeningAngle1_refitted);
528 gg[5]->Fill(OpeningAngle2_refitted);
529
530 InvMass_realg1_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(i), MC.at(l));
531 InvMass_realg2_mc = CbmKresFunctions::Invmass_2particles_MC(MC.at(j), MC.at(k));
532 gg[2]->Fill(InvMass_realg1_mc);
533 gg[2]->Fill(InvMass_realg2_mc);
534 InvMass_realg1_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(i), RefMom.at(l));
535 InvMass_realg2_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(j), RefMom.at(k));
536 gg[3]->Fill(InvMass_realg1_refitted);
537 gg[3]->Fill(InvMass_realg2_refitted);
538 openingAngleBetweenGammasReco =
539 CalculateOpeningAngleBetweenGammas_Reco(RefMom.at(i), RefMom.at(l), RefMom.at(j), RefMom.at(k));
540 }
541 Double_t openingAngleBetweenGammas =
542 CalculateOpeningAngleBetweenGammas_MC(MC.at(i), MC.at(j), MC.at(k), MC.at(l));
543 gg[6]->Fill(openingAngleBetweenGammas);
544 gg[7]->Fill(openingAngleBetweenGammasReco);
545 }
546 }
547 }
548 }
549 }
550}
551
552
553void CbmKresEtaMCAnalysis::EtaChargedPionsGammaAnalysis(vector<TVector3> RefMomPion, vector<CbmMCTrack*> MCPion,
554 vector<Int_t> /*IdPion*/, vector<TVector3> RefMomEl,
555 vector<CbmMCTrack*> MCEl, vector<Int_t> /*IdEl*/,
556 vector<TH1*> ppg)
557{
558 //================================== decay eta -> p+p- gamma -> p+p- e+e-
559 if (MCPion.size() < 2 || MCEl.size() < 2) return;
560
561 // int fDecayedParticlePdg = 221;
562
563 for (size_t i = 0; i < MCEl.size(); i++) {
564 for (size_t j = i + 1; j < MCEl.size(); j++) {
565 int pdg1 = MCEl.at(i)->GetPdgCode();
566 int pdg2 = MCEl.at(j)->GetPdgCode();
567
568 if (pdg1 + pdg2 != 0) continue;
569 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11) continue;
570
571 int motherId1 = MCEl.at(i)->GetMotherId();
572 int motherId2 = MCEl.at(j)->GetMotherId();
573
574 if (motherId1 != motherId2) continue;
575 if (motherId1 == -1 || motherId2 == -1) continue;
576
577 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
578
579 if (mother1->GetPdgCode() != 22) continue;
580
581 int grandmotherId1 = mother1->GetMotherId();
582
583 if (grandmotherId1 == -1) continue;
584
585 CbmMCTrack* GrMother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
586 int mcGrandmotherPdg1 = GrMother1->GetPdgCode();
587 if (mcGrandmotherPdg1 != 221) continue;
588
589
590 for (size_t k = 0; k < MCPion.size(); k++) {
591 for (size_t l = k + 1; l < MCPion.size(); l++) {
592 int pdg3 = MCPion.at(k)->GetPdgCode();
593 int pdg4 = MCPion.at(l)->GetPdgCode();
594
595 if (pdg3 + pdg4 != 0) continue;
596 if (TMath::Abs(pdg3) != 211 || TMath::Abs(pdg4) != 211) continue;
597
598 int motherId3 = MCPion.at(k)->GetMotherId();
599 int motherId4 = MCPion.at(l)->GetMotherId();
600
601 if (motherId3 != motherId4) continue;
602 if (motherId3 == -1 || motherId4 == -1) continue;
603 if (motherId3 != grandmotherId1) continue;
604
605 Double_t InvMass_true =
606 CbmKresFunctions::Invmass_4particles_MC(MCEl.at(i), MCEl.at(j), MCPion.at(k), MCPion.at(l));
607 Double_t InvMass_reco = CbmKresFunctions::Invmass_2el_2pions_RECO(RefMomEl.at(i), RefMomEl.at(j),
608 RefMomPion.at(k), RefMomPion.at(l));
609 cout << "Decay eta -> p+p- gamma -> p+p- e+e- detected!\t\t mc mass: " << InvMass_true
610 << "\t, reco mass: " << InvMass_reco << endl;
611 cout << "motherids: " << motherId1 << "/" << motherId2 << "/" << motherId3 << "/" << motherId4 << endl;
612 cout << "grandmotherid: " << grandmotherId1 << endl;
613 cout << "pdgs: " << mcGrandmotherPdg1 << "-->" << mother1->GetPdgCode() << "/" << pdg3 << "/" << pdg4 << "-->"
614 << pdg1 << "/" << pdg2 << "/" << pdg3 << "/" << pdg4 << endl;
615
616 ppg[0]->Fill(InvMass_true);
617 ppg[1]->Fill(InvMass_reco);
618
619 cout << "\t \t mc true info: " << endl;
620 cout << "particle 1: \t" << MCEl.at(i)->GetPdgCode() << ";\t pt = " << MCEl.at(i)->GetPt()
621 << ";\t X = " << MCEl.at(i)->GetStartX() << ";\t Y = " << MCEl.at(i)->GetStartY()
622 << ";\t Z = " << MCEl.at(i)->GetStartZ() << ";\t E = " << MCEl.at(i)->GetEnergy() << endl;
623 cout << "particle 2: \t" << MCEl.at(j)->GetPdgCode() << ";\t pt = " << MCEl.at(j)->GetPt()
624 << ";\t X = " << MCEl.at(j)->GetStartX() << ";\t Y = " << MCEl.at(j)->GetStartY()
625 << ";\t Z = " << MCEl.at(j)->GetStartZ() << ";\t E = " << MCEl.at(j)->GetEnergy() << endl;
626 cout << "particle 3: \t" << MCPion.at(k)->GetPdgCode() << ";\t pt = " << MCPion.at(k)->GetPt()
627 << ";\t X = " << MCPion.at(k)->GetStartX() << ";\t Y = " << MCPion.at(k)->GetStartY()
628 << ";\t Z = " << MCPion.at(k)->GetStartZ() << ";\t E = " << MCPion.at(k)->GetEnergy() << endl;
629 cout << "particle 4: \t" << MCPion.at(l)->GetPdgCode() << ";\t pt = " << MCPion.at(l)->GetPt()
630 << ";\t X = " << MCPion.at(l)->GetStartX() << ";\t Y = " << MCPion.at(l)->GetStartY()
631 << ";\t Z = " << MCPion.at(l)->GetStartZ() << ";\t E = " << MCPion.at(l)->GetEnergy() << endl;
632
633 double InvMass_g_mc = CbmKresFunctions::Invmass_2particles_MC(MCEl.at(i), MCEl.at(j));
634 double InvMass_g_reffited = CbmKresFunctions::Invmass_2particles_RECO(RefMomEl.at(i), RefMomEl.at(j));
635 double OpeningAngle_g_mc = CbmKresFunctions::CalculateOpeningAngle_MC(MCEl.at(i), MCEl.at(j));
636 double OpeningAngle_g_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMomEl.at(i), RefMomEl.at(j));
637
638 ppg[2]->Fill(InvMass_g_mc);
639 ppg[3]->Fill(InvMass_g_reffited);
640 ppg[4]->Fill(OpeningAngle_g_mc);
641 ppg[5]->Fill(OpeningAngle_g_refitted);
642 }
643 }
644 }
645 }
646}
647
648
649void CbmKresEtaMCAnalysis::EtaPosNegNeutralPionsAnalysis(vector<TVector3> RefMomNeutral, vector<CbmMCTrack*> MCNeutral,
650 vector<Int_t> IdNeutral, vector<TVector3> RefMomPion,
651 vector<CbmMCTrack*> MCPion, vector<Int_t> /*IdPion*/,
652 vector<TH1*> ppp)
653{
654 if (MCNeutral.size() < 4 || MCPion.size() < 2) return;
655
656 int fDecayedParticlePdg = 221;
657
658 // cout << "!!!!!!!!!! MCNeutral.size() = " << MCNeutral.size() << endl;
659
660 for (size_t i = 0; i < MCNeutral.size(); i++) {
661 for (size_t j = i + 1; j < MCNeutral.size(); j++) {
662 for (size_t k = j + 1; k < MCNeutral.size(); k++) {
663 for (size_t l = k + 1; l < MCNeutral.size(); l++) {
664
665 int pdg1 = MCNeutral.at(i)->GetPdgCode();
666 int pdg2 = MCNeutral.at(j)->GetPdgCode();
667 int pdg3 = MCNeutral.at(k)->GetPdgCode();
668 int pdg4 = MCNeutral.at(l)->GetPdgCode();
669
670 if (pdg1 + pdg2 + pdg3 + pdg4 != 0) continue;
671 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11 || TMath::Abs(pdg3) != 11 || TMath::Abs(pdg4) != 11)
672 continue;
673
674 int motherId1 = MCNeutral.at(i)->GetMotherId();
675 int motherId2 = MCNeutral.at(j)->GetMotherId();
676 int motherId3 = MCNeutral.at(k)->GetMotherId();
677 int motherId4 = MCNeutral.at(l)->GetMotherId();
678
679 int mcId1 = IdNeutral.at(i);
680 int mcId2 = IdNeutral.at(j);
681 int mcId3 = IdNeutral.at(k);
682 int mcId4 = IdNeutral.at(l);
683
684 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
685 continue; // particle is not used twice
686 if (motherId1 == -1 || motherId2 == -1 || motherId3 == -1 || motherId4 == -1) continue;
687
688 CbmMCTrack* mother1 = (CbmMCTrack*) fMcTracks->At(motherId1);
689 CbmMCTrack* mother2 = (CbmMCTrack*) fMcTracks->At(motherId2);
690 CbmMCTrack* mother3 = (CbmMCTrack*) fMcTracks->At(motherId3);
691 CbmMCTrack* mother4 = (CbmMCTrack*) fMcTracks->At(motherId4);
692
693 int mcMotherPdg1 = mother1->GetPdgCode();
694 int mcMotherPdg2 = mother2->GetPdgCode();
695 int mcMotherPdg3 = mother3->GetPdgCode();
696 int mcMotherPdg4 = mother4->GetPdgCode();
697
698 if (mcMotherPdg1 != 22 || mcMotherPdg2 != 22 || mcMotherPdg3 != 22 || mcMotherPdg4 != 22) continue;
699
700 int grandmotherId1 = mother1->GetMotherId();
701 int grandmotherId2 = mother2->GetMotherId();
702 int grandmotherId3 = mother3->GetMotherId();
703 int grandmotherId4 = mother4->GetMotherId();
704
705 if (grandmotherId1 == -1) continue;
706
707 if (grandmotherId1 == grandmotherId2 && grandmotherId1 == grandmotherId3
708 && grandmotherId1 == grandmotherId4) {
709 CbmMCTrack* GrMother1 = (CbmMCTrack*) fMcTracks->At(grandmotherId1);
710 int mcGrandmotherPdg1 = GrMother1->GetPdgCode();
711 if (mcGrandmotherPdg1 != 111) continue;
712 // cout << "***************** 111 ****************" << endl;
713 if (GrMother1->GetMotherId() == -1) continue;
714 CbmMCTrack* TopTrack = (CbmMCTrack*) fMcTracks->At(GrMother1->GetMotherId());
715 if (TopTrack->GetPdgCode() != fDecayedParticlePdg) continue;
716 // cout << "***************** from eta !!!!! ****************" << endl;
717
718 for (size_t m = 0; m < MCPion.size(); m++) {
719 for (size_t s = m + 1; s < MCPion.size(); s++) {
720 int pdg5 = MCPion.at(m)->GetPdgCode();
721 int pdg6 = MCPion.at(s)->GetPdgCode();
722
723 if (pdg5 + pdg6 != 0) continue;
724 if (TMath::Abs(pdg5) != 211 || TMath::Abs(pdg6) != 211) continue;
725
726 int motherId5 = MCPion.at(m)->GetMotherId();
727 int motherId6 = MCPion.at(s)->GetMotherId();
728
729 if (motherId5 != motherId6) continue;
730 if (motherId5 == -1 || motherId6 == -1) continue;
731 if (motherId5 != GrMother1->GetMotherId()) continue;
732
733 Double_t InvMass_six_true = CbmKresFunctions::Invmass_6particles_MC(
734 MCNeutral.at(i), MCNeutral.at(j), MCNeutral.at(k), MCNeutral.at(l), MCPion.at(m), MCPion.at(s));
735 Double_t InvMass_six_reco = CbmKresFunctions::Invmass_4el_2pions_RECO(
736 RefMomNeutral.at(i), RefMomNeutral.at(j), RefMomNeutral.at(k), RefMomNeutral.at(l), RefMomPion.at(m),
737 RefMomPion.at(s));
738 cout << "Decay eta -> p+p- p0 -> p+p- e+e-e+e- detected!\t\t "
739 "mc mass: "
740 << InvMass_six_true << "\t, reco mass: " << InvMass_six_reco << endl;
741 cout << "motherids for neutral: " << motherId1 << "/" << motherId2 << "/" << motherId3 << "/"
742 << motherId4 << endl;
743 cout << "grandmotherid: " << GrMother1->GetMotherId() << "/" << motherId5 << "/" << motherId6 << endl;
744 cout << "pdgs: " << TopTrack->GetPdgCode() << "-->" << GrMother1->GetPdgCode() << "/" << pdg5 << "/"
745 << pdg6 << endl;
746
747 ppp[0]->Fill(InvMass_six_true);
748 ppp[1]->Fill(InvMass_six_reco);
749
750 cout << "\t \t mc true info: " << endl;
751 cout << "particle 1: \t" << MCNeutral.at(i)->GetPdgCode() << ";\t pt = " << MCNeutral.at(i)->GetPt()
752 << ";\t X = " << MCNeutral.at(i)->GetStartX() << ";\t Y = " << MCNeutral.at(i)->GetStartY()
753 << ";\t Z = " << MCNeutral.at(i)->GetStartZ() << ";\t E = " << MCNeutral.at(i)->GetEnergy()
754 << endl;
755 cout << "particle 2: \t" << MCNeutral.at(j)->GetPdgCode() << ";\t pt = " << MCNeutral.at(j)->GetPt()
756 << ";\t X = " << MCNeutral.at(j)->GetStartX() << ";\t Y = " << MCNeutral.at(j)->GetStartY()
757 << ";\t Z = " << MCNeutral.at(j)->GetStartZ() << ";\t E = " << MCNeutral.at(j)->GetEnergy()
758 << endl;
759 cout << "particle 3: \t" << MCNeutral.at(k)->GetPdgCode() << ";\t pt = " << MCNeutral.at(k)->GetPt()
760 << ";\t X = " << MCNeutral.at(k)->GetStartX() << ";\t Y = " << MCNeutral.at(k)->GetStartY()
761 << ";\t Z = " << MCNeutral.at(k)->GetStartZ() << ";\t E = " << MCNeutral.at(k)->GetEnergy()
762 << endl;
763 cout << "particle 4: \t" << MCNeutral.at(l)->GetPdgCode() << ";\t pt = " << MCNeutral.at(l)->GetPt()
764 << ";\t X = " << MCNeutral.at(l)->GetStartX() << ";\t Y = " << MCNeutral.at(l)->GetStartY()
765 << ";\t Z = " << MCNeutral.at(l)->GetStartZ() << ";\t E = " << MCNeutral.at(l)->GetEnergy()
766 << endl;
767 cout << "particle 5: \t" << MCPion.at(m)->GetPdgCode() << ";\t pt = " << MCPion.at(m)->GetPt()
768 << ";\t X = " << MCPion.at(m)->GetStartX() << ";\t Y = " << MCPion.at(m)->GetStartY()
769 << ";\t Z = " << MCPion.at(m)->GetStartZ() << ";\t E = " << MCPion.at(m)->GetEnergy() << endl;
770 cout << "particle 6: \t" << MCPion.at(s)->GetPdgCode() << ";\t pt = " << MCPion.at(s)->GetPt()
771 << ";\t X = " << MCPion.at(s)->GetStartX() << ";\t Y = " << MCPion.at(s)->GetStartY()
772 << ";\t Z = " << MCPion.at(s)->GetStartZ() << ";\t E = " << MCPion.at(s)->GetEnergy() << endl;
773
774 Double_t InvMass_true = CbmKresFunctions::Invmass_4particles_MC(MCNeutral.at(i), MCNeutral.at(j),
775 MCNeutral.at(k), MCNeutral.at(l));
776 Double_t InvMass_reco = CbmKresFunctions::Invmass_4particles_RECO(
777 RefMomNeutral.at(i), RefMomNeutral.at(j), RefMomNeutral.at(k), RefMomNeutral.at(l));
778 ppp[2]->Fill(InvMass_true);
779 ppp[3]->Fill(InvMass_reco);
780 }
781 }
782 }
783 }
784 }
785 }
786 }
787}
788
789
790void CbmKresEtaMCAnalysis::EtaDoubleGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut, int Event,
791 vector<TVector3> RefMom, vector<CbmMCTrack*> MC,
792 vector<Int_t> Id, vector<TH1*> gg)
793{
794 if (MC.size() < 4) return;
795
796 EDGA_RefMom.clear();
797 EDGA_MC.clear();
798 EDGA_Id.clear();
799
800 for (size_t i = 0; i < MC.size(); i++) {
801 for (size_t j = i + 1; j < MC.size(); j++) {
802 int pdg1 = MC.at(i)->GetPdgCode();
803 int pdg2 = MC.at(j)->GetPdgCode();
804 if (pdg1 + pdg2 != 0) continue;
805 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11) continue;
806
807 double OpeningAngle_g_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMom.at(i), RefMom.at(j));
808 double InvMass_g_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMom.at(i), RefMom.at(j));
809
810 if (OpeningAngle_g_refitted > OpeningAngleCut || InvMass_g_refitted > GammaInvMassCut) continue;
811
812 frefmomenta.clear();
813 frefmomenta.push_back(RefMom.at(i));
814 frefmomenta.push_back(RefMom.at(j));
815 fMCtracks.clear();
816 fMCtracks.push_back(MC.at(i));
817 fMCtracks.push_back(MC.at(j));
818 fMCId.clear();
819 fMCId.push_back(Id.at(i));
820 fMCId.push_back(Id.at(j));
821
822 EDGA_RefMom.push_back(frefmomenta);
823 EDGA_MC.push_back(fMCtracks);
824 EDGA_Id.push_back(fMCId);
825
826 // for event mixing gg channel
827 EMT_gg_Event.push_back(Event);
829 }
830 }
831
832 if (EDGA_MC.size() < 2) return; // min 2 gammas to form eta are required
833 for (size_t gamma1 = 0; gamma1 < EDGA_RefMom.size(); gamma1++) {
834 for (size_t gamma2 = gamma1 + 1; gamma2 < EDGA_RefMom.size(); gamma2++) {
835 TVector3 e11 = EDGA_RefMom[gamma1][0];
836 TVector3 e12 = EDGA_RefMom[gamma1][1];
837 TVector3 e21 = EDGA_RefMom[gamma2][0];
838 TVector3 e22 = EDGA_RefMom[gamma2][1];
839
840 // MC true data for this particles
841 CbmMCTrack* mcTrack1 = EDGA_MC[gamma1][0];
842 CbmMCTrack* mcTrack2 = EDGA_MC[gamma1][1];
843 CbmMCTrack* mcTrack3 = EDGA_MC[gamma2][0];
844 CbmMCTrack* mcTrack4 = EDGA_MC[gamma2][1];
845
846 int mcId1 = EDGA_Id[gamma1][0];
847 int mcId2 = EDGA_Id[gamma1][1];
848 int mcId3 = EDGA_Id[gamma2][0];
849 int mcId4 = EDGA_Id[gamma2][1];
850
851 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
852 continue; // particle is not used twice
853
854 double openingAngleBetweenGammasReco = CalculateOpeningAngleBetweenGammas_Reco(e11, e12, e21, e22);
855 if (openingAngleBetweenGammasReco < 10 || openingAngleBetweenGammasReco > 40) continue;
856
857 double InvMass_true = CbmKresFunctions::Invmass_4particles_MC(mcTrack1, mcTrack2, mcTrack3, mcTrack4);
858 double InvMass_reco = CbmKresFunctions::Invmass_4particles_RECO(e11, e12, e21, e22);
860
861 gg[8]->Fill(InvMass_true);
862 gg[9]->Fill(InvMass_reco);
863
864 int CorrectEta = 0;
865 if ((mcTrack1->GetMotherId() == mcTrack2->GetMotherId() && mcTrack3->GetMotherId() == mcTrack4->GetMotherId())
866 || (mcTrack1->GetMotherId() == mcTrack3->GetMotherId() && mcTrack2->GetMotherId() == mcTrack4->GetMotherId())
867 || (mcTrack1->GetMotherId() == mcTrack4->GetMotherId()
868 && mcTrack2->GetMotherId() == mcTrack3->GetMotherId())) {
869 CbmMCTrack* MothTrack1 = (CbmMCTrack*) fMcTracks->At(mcTrack1->GetMotherId());
870 CbmMCTrack* MothTrack2 = (CbmMCTrack*) fMcTracks->At(mcTrack2->GetMotherId());
871 CbmMCTrack* MothTrack3 = (CbmMCTrack*) fMcTracks->At(mcTrack3->GetMotherId());
872 CbmMCTrack* MothTrack4 = (CbmMCTrack*) fMcTracks->At(mcTrack4->GetMotherId());
873 if (MothTrack1->GetPdgCode() == 22 && MothTrack2->GetPdgCode() == 22 && MothTrack3->GetPdgCode() == 22
874 && MothTrack4->GetPdgCode() == 22 && MothTrack1->GetMotherId() != -1
875 && MothTrack1->GetMotherId() == MothTrack2->GetMotherId()
876 && MothTrack1->GetMotherId() == MothTrack3->GetMotherId()
877 && MothTrack1->GetMotherId() == MothTrack4->GetMotherId()) {
878 CbmMCTrack* GrandMothTrack1 = (CbmMCTrack*) fMcTracks->At(MothTrack1->GetMotherId());
879 if (GrandMothTrack1->GetPdgCode() == 221) CorrectEta = 1;
880 }
881 }
882
883 if (CorrectEta == 1) {
884 InvMass_eta_gg_reco_aftercuts->Fill(InvMass_reco);
885 rap_vs_pt_eta_gg_reco_aftercuts->Fill(params.fRapidity, params.fPt);
886 }
887 else {
889 }
890 }
891 }
892}
893
894void CbmKresEtaMCAnalysis::EtaChargedPionsGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut,
895 int Event, vector<TVector3> RefMomPion,
896 vector<CbmMCTrack*> MCPion, vector<Int_t> /*IdPion*/,
897 vector<TVector3> RefMomEl, vector<CbmMCTrack*> MCEl,
898 vector<Int_t> IdEl, vector<TH1*> ppg, vector<TH1*> ppp)
899{
900 if (MCPion.size() < 2 || MCEl.size() < 2) return;
901
902 ECPGA_leptons_RefMom.clear();
903 ECPGA_leptons_MC.clear();
904 ECPGA_leptons_Id.clear();
905
906 for (size_t i = 0; i < MCEl.size(); i++) {
907 for (size_t j = i + 1; j < MCEl.size(); j++) {
908 int pdg1 = MCEl.at(i)->GetPdgCode();
909 int pdg2 = MCEl.at(j)->GetPdgCode();
910 if (pdg1 + pdg2 != 0) continue;
911 if (TMath::Abs(pdg1) != 11 || TMath::Abs(pdg2) != 11) continue;
912
913 double OpeningAngle_g_refitted = CbmKresFunctions::CalculateOpeningAngle_Reco(RefMomEl.at(i), RefMomEl.at(j));
914 double InvMass_g_refitted = CbmKresFunctions::Invmass_2particles_RECO(RefMomEl.at(i), RefMomEl.at(j));
915
916 if (OpeningAngle_g_refitted > OpeningAngleCut || InvMass_g_refitted > GammaInvMassCut) continue;
917
918 frefmomenta.clear();
919 frefmomenta.push_back(RefMomEl.at(i));
920 frefmomenta.push_back(RefMomEl.at(j));
921 fMCtracks.clear();
922 fMCtracks.push_back(MCEl.at(i));
923 fMCtracks.push_back(MCEl.at(j));
924 fMCId.clear();
925 fMCId.push_back(IdEl.at(i));
926 fMCId.push_back(IdEl.at(j));
927
929 ECPGA_leptons_MC.push_back(fMCtracks);
930 ECPGA_leptons_Id.push_back(fMCId);
931
932 // for event mixing ppg channel
933 EMT_ppg_ee_Event.push_back(Event);
935 }
936 }
937
938 ECPGA_pions_RefMom.clear();
939 ECPGA_pions_MC.clear();
940 for (size_t k = 0; k < MCPion.size(); k++) {
941 if (MCPion.at(k)->GetPdgCode() == 211) {
942 EMT_ppg_positive_pion_Event.push_back(Event);
943 EMT_ppg_positive_pion_momenta.push_back(RefMomPion.at(k));
944 }
945 else {
946 EMT_ppg_negative_pion_Event.push_back(Event);
947 EMT_ppg_negative_pion_momenta.push_back(RefMomPion.at(k));
948 }
949 for (size_t l = k + 1; l < MCPion.size(); l++) {
950 int pdg3 = MCPion.at(k)->GetPdgCode();
951 int pdg4 = MCPion.at(l)->GetPdgCode();
952 if (pdg3 + pdg4 != 0) continue;
953 if (TMath::Abs(pdg3) != 211 || TMath::Abs(pdg4) != 211) continue;
954
955 double Pions_angle_true = CbmKresFunctions::CalculateOpeningAngleBetweenPions_MC(MCPion.at(k), MCPion.at(l));
956 double Pions_angle_reco =
957 CbmKresFunctions::CalculateOpeningAngleBetweenPions_Reco(RefMomPion.at(k), RefMomPion.at(l));
958
959 if (Pions_angle_reco > 20) continue;
960
961 CbmMCTrack* mcTrack3 = MCPion.at(k);
962 CbmMCTrack* mcTrack4 = MCPion.at(l);
963
964 frefmomenta.clear();
965 frefmomenta.push_back(RefMomPion.at(k));
966 frefmomenta.push_back(RefMomPion.at(l));
967
968 fMCtracks.clear();
969 fMCtracks.push_back(MCPion.at(k));
970 fMCtracks.push_back(MCPion.at(l));
971
973 ECPGA_pions_MC.push_back(fMCtracks);
974
975 // for event mixing ppg channel
976 EMT_ppg_pp_Event.push_back(Event);
978
979 // if (pdg3 == 211){
980 // EMT_ppg_positive_pion_Event.push_back(Event);
981 // EMT_ppg_positive_pion_momenta.push_back(RefMomPion.at(k));
982 // EMT_ppg_negative_pion_Event.push_back(Event);
983 // EMT_ppg_negative_pion_momenta.push_back(RefMomPion.at(l));
984 // }
985 // if (pdg4 == 211){
986 // EMT_ppg_negative_pion_Event.push_back(Event);
987 // EMT_ppg_negative_pion_momenta.push_back(RefMomPion.at(k));
988 // EMT_ppg_positive_pion_Event.push_back(Event);
989 // EMT_ppg_positive_pion_momenta.push_back(RefMomPion.at(l));
990 // }
991
992 int correctPionspair = 0;
993 if (mcTrack3->GetMotherId() == mcTrack4->GetMotherId() && mcTrack4->GetMotherId() != -1) {
994 CbmMCTrack* EtaTrack = (CbmMCTrack*) fMcTracks->At(mcTrack4->GetMotherId());
995 if (EtaTrack->GetPdgCode() == 221) correctPionspair = 1;
996 }
997
998 if (correctPionspair == 1) {
999 OA_betweenPions_fromEta_mc->Fill(Pions_angle_true);
1000 OA_betweenPions_fromEta_reco->Fill(Pions_angle_reco);
1001 }
1002 else {
1003 OA_betweenPions_fromEta_reco_wrongcombinations->Fill(Pions_angle_reco);
1004 }
1005 }
1006 }
1007
1008
1009 if (ECPGA_leptons_MC.size() < 1) return; // min 1 gamma to form eta is required
1010 for (size_t gamma = 0; gamma < ECPGA_leptons_RefMom.size(); gamma++) {
1011 if (ECPGA_pions_MC.size() < 1) continue; // min 1 p+p- pair to form eta is required
1012 TVector3 e11 = ECPGA_leptons_RefMom[gamma][0];
1013 TVector3 e12 = ECPGA_leptons_RefMom[gamma][1];
1014
1015 CbmMCTrack* mcTrack1 = ECPGA_leptons_MC[gamma][0];
1016 CbmMCTrack* mcTrack2 = ECPGA_leptons_MC[gamma][1];
1017
1018 for (size_t pionpair = 0; pionpair < ECPGA_pions_RefMom.size(); pionpair++) {
1019
1020 CbmMCTrack* mcTrack3 = ECPGA_pions_MC[pionpair][0];
1021 CbmMCTrack* mcTrack4 = ECPGA_pions_MC[pionpair][1];
1022
1023 TVector3 e21 = ECPGA_pions_RefMom[pionpair][0];
1024 TVector3 e22 = ECPGA_pions_RefMom[pionpair][1];
1025
1026 Double_t InvMass_true = CbmKresFunctions::Invmass_4particles_MC(mcTrack1, mcTrack2, mcTrack3, mcTrack4);
1027 Double_t InvMass_reco = CbmKresFunctions::Invmass_2el_2pions_RECO(e11, e12, e21, e22);
1029
1030 ppg[6]->Fill(InvMass_true);
1031 ppg[7]->Fill(InvMass_reco);
1032
1033 int CorrectEta = 0;
1034 if (mcTrack1->GetMotherId() == mcTrack2->GetMotherId() && mcTrack1->GetMotherId() != -1) {
1035 CbmMCTrack* GammaTrack = (CbmMCTrack*) fMcTracks->At(mcTrack1->GetMotherId());
1036 if (GammaTrack->GetPdgCode() == 22 && GammaTrack->GetMotherId() != -1
1037 && mcTrack3->GetMotherId() == mcTrack4->GetMotherId()
1038 && mcTrack3->GetMotherId() == GammaTrack->GetMotherId()) {
1039 CbmMCTrack* EtaTrack = (CbmMCTrack*) fMcTracks->At(GammaTrack->GetMotherId());
1040 if (EtaTrack->GetPdgCode() == 221) CorrectEta = 1;
1041 }
1042 }
1043
1044 if (CorrectEta == 1) {
1045 InvMass_eta_ppg_reco_aftercuts->Fill(InvMass_reco);
1046 rap_vs_pt_eta_ppg_reco_aftercuts->Fill(params.fRapidity, params.fPt);
1047 }
1048 else {
1050 }
1051 }
1052 }
1053
1054
1056 if (ECPGA_leptons_MC.size() < 2)
1057 return; // min 2 gammas to form pion is required in EtaPosNegNeutralPionsAnalysis_plusBG
1058 for (size_t gamma1 = 0; gamma1 < ECPGA_leptons_RefMom.size(); gamma1++) {
1059 for (size_t gamma2 = gamma1 + 1; gamma2 < ECPGA_leptons_RefMom.size(); gamma2++) {
1060 TVector3 e11 = ECPGA_leptons_RefMom[gamma1][0];
1061 TVector3 e12 = ECPGA_leptons_RefMom[gamma1][1];
1062 TVector3 e21 = ECPGA_leptons_RefMom[gamma2][0];
1063 TVector3 e22 = ECPGA_leptons_RefMom[gamma2][1];
1064
1065 CbmMCTrack* mcTrack1 = ECPGA_leptons_MC[gamma1][0];
1066 CbmMCTrack* mcTrack2 = ECPGA_leptons_MC[gamma1][1];
1067 CbmMCTrack* mcTrack3 = ECPGA_leptons_MC[gamma2][0];
1068 CbmMCTrack* mcTrack4 = ECPGA_leptons_MC[gamma2][1];
1069
1070 int mcId1 = ECPGA_leptons_Id[gamma1][0];
1071 int mcId2 = ECPGA_leptons_Id[gamma1][1];
1072 int mcId3 = ECPGA_leptons_Id[gamma2][0];
1073 int mcId4 = ECPGA_leptons_Id[gamma2][1];
1074
1075 if (mcId1 == mcId2 || mcId1 == mcId3 || mcId1 == mcId4 || mcId2 == mcId3 || mcId2 == mcId4 || mcId3 == mcId4)
1076 continue; // particle is not used twice
1077
1078 // Double_t InvMass_pi_true = CbmKresFunctions::Invmass_4particles_MC(mcTrack1, mcTrack2, mcTrack3, mcTrack4);
1079 CbmKresFunctions::Invmass_4particles_MC(mcTrack1, mcTrack2, mcTrack3, mcTrack4);
1080 Double_t InvMass_pi_reco = CbmKresFunctions::Invmass_4particles_RECO(e11, e12, e21, e22);
1081
1082 if (InvMass_pi_reco < 0.12 || InvMass_pi_reco > 0.15) continue;
1083
1084 for (size_t k = 0; k < MCPion.size(); k++) {
1085 for (size_t l = k + 1; l < MCPion.size(); l++) {
1086 int pdg5 = MCPion.at(k)->GetPdgCode();
1087 int pdg6 = MCPion.at(l)->GetPdgCode();
1088
1089 if (pdg5 + pdg6 != 0) continue;
1090 if (TMath::Abs(pdg5) != 211 || TMath::Abs(pdg6) != 211) continue;
1091
1092 double InvMass_six_true =
1093 CbmKresFunctions::Invmass_6particles_MC(mcTrack1, mcTrack2, mcTrack3, mcTrack4, MCPion.at(k), MCPion.at(l));
1094 double InvMass_six_reco =
1095 CbmKresFunctions::Invmass_4el_2pions_RECO(e11, e12, e21, e22, RefMomPion.at(k), RefMomPion.at(l));
1096 double Pions_angle_reco =
1097 CbmKresFunctions::CalculateOpeningAngleBetweenPions_Reco(RefMomPion.at(k), RefMomPion.at(l));
1098
1099 if (Pions_angle_reco > 20) continue;
1100
1101 ppp[4]->Fill(InvMass_six_true);
1102 ppp[5]->Fill(InvMass_six_reco);
1103 }
1104 }
1105 }
1106 }
1107}
1108
1109
1111{
1112 int nof = EMT_gg_Event.size();
1113 cout << "Mixing for eta->gg channel - nof entries " << nof << endl;
1114 for (Int_t a = 0; a < nof - 1; a++) {
1115 for (Int_t b = a + 1; b < nof; b++) {
1116 if (EMT_gg_Event[a] == EMT_gg_Event[b]) continue; // to make sure that the photons are from two different events
1117 TVector3 e11 = EMT_gg_pair_momenta[a][0];
1118 TVector3 e12 = EMT_gg_pair_momenta[a][1];
1119 TVector3 e21 = EMT_gg_pair_momenta[b][0];
1120 TVector3 e22 = EMT_gg_pair_momenta[b][1];
1121
1122 double openingAngleBetweenGammasReco = CalculateOpeningAngleBetweenGammas_Reco(e11, e12, e21, e22);
1123 if (openingAngleBetweenGammasReco < 10 || openingAngleBetweenGammasReco > 40) continue;
1124
1126 EMT_eta_gg->Fill(params.fMinv);
1127 }
1128 }
1129}
1130
1132{
1133 int nof_leptons = EMT_ppg_ee_Event.size();
1134 int nof_pions = EMT_ppg_pp_Event.size();
1135 cout << "Mixing for eta->(pp)g channel - nof lepton pairs = " << nof_leptons << "; nof pion pairs = " << nof_pions
1136 << endl;
1137 for (Int_t a = 0; a < nof_leptons; a++) {
1138 for (Int_t b = 0; b < nof_pions; b++) {
1140 continue; // to make sure that the photons are from two different events
1141 TVector3 e11 = EMT_ppg_ee_pair_momenta[a][0];
1142 TVector3 e12 = EMT_ppg_ee_pair_momenta[a][1];
1143 TVector3 e21 = EMT_ppg_pp_pair_momenta[b][0];
1144 TVector3 e22 = EMT_ppg_pp_pair_momenta[b][1];
1145
1147 EMT_eta_ppg->Fill(params.fMinv);
1148 }
1149 }
1150}
1151
1152
1154{
1155 int nof_photons = EMT_ppg_ee_Event.size();
1156 int nof_positive = EMT_ppg_positive_pion_Event.size();
1157 int nof_negative = EMT_ppg_negative_pion_Event.size();
1158 cout << "Mixing 3 bodies for eta-> p p g channel - nof photons = " << nof_photons
1159 << "; nof +pions = " << nof_positive << "; nof -pions = " << nof_negative << endl;
1160
1161 for (Int_t p = 0; p < nof_positive; p++) {
1162 for (Int_t m = 0; m < nof_negative; m++) {
1164 continue; // to make sure that the photons are from two different events
1165 TVector3 e21 = EMT_ppg_positive_pion_momenta[p];
1166 TVector3 e22 = EMT_ppg_negative_pion_momenta[m];
1167
1168 double Pions_angle_reco = CbmKresFunctions::CalculateOpeningAngleBetweenPions_Reco(e21, e22);
1169 if (Pions_angle_reco > 20) continue;
1170
1171 for (Int_t a = 0; a < nof_photons; a++) {
1174 continue; // to make sure that the photons are from two different events
1175 TVector3 e11 = EMT_ppg_ee_pair_momenta[a][0];
1176 TVector3 e12 = EMT_ppg_ee_pair_momenta[a][1];
1177
1179 EMT_eta_three_body->Fill(params.fMinv);
1180 }
1181 }
1182 }
1183}
1184
1185
1187{
1188 gDirectory->mkdir("EtaMCAnalysis");
1189 gDirectory->cd("EtaMCAnalysis");
1190
1191 gDirectory->mkdir("eta_gg");
1192 gDirectory->cd("eta_gg");
1193 for (UInt_t i = 0; i < fHistoList_eta_gg.size(); i++) {
1194 fHistoList_eta_gg[i]->Write();
1195 }
1196 gDirectory->cd("..");
1197
1198 gDirectory->mkdir("eta_ppg");
1199 gDirectory->cd("eta_ppg");
1200 for (UInt_t i = 0; i < fHistoList_eta_ppg.size(); i++) {
1201 fHistoList_eta_ppg[i]->Write();
1202 }
1203 gDirectory->cd("..");
1204
1205 gDirectory->mkdir("eta_ppp");
1206 gDirectory->cd("eta_ppp");
1207 for (UInt_t i = 0; i < fHistoList_eta_ppp.size(); i++) {
1208 fHistoList_eta_ppp[i]->Write();
1209 }
1210 gDirectory->cd("..");
1211
1212
1213 gDirectory->cd("..");
1214}
1215
1217{
1219 InvMass_eta_gg_mc = new TH1D("InvMass_eta_gg_mc", "InvMass_eta_gg_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1222 new TH1D("InvMass_eta_gg_reffited", "InvMass_eta_gg_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1225 new TH1D("InvMassPhoton_eta_gg_mc", "InvMassPhoton_eta_gg_mc; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1228 "InvMassPhoton_eta_gg_reffited", "InvMassPhoton_eta_gg_reffited; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1231 new TH1D("OpeningAnglePhoton_eta_gg_mc", "OpeningAnglePhoton_eta_gg_mc (between e+e- from #gamma); angle [deg];#",
1232 100, -0.1, 9.9);
1234 OpeningAnglePhoton_eta_gg_reffited = new TH1D("OpeningAnglePhoton_eta_gg_reffited",
1235 "OpeningAnglePhoton_eta_gg_reffited (between e+e- from #gamma); "
1236 "angle [deg];#",
1237 100, -0.1, 9.9);
1239 OpeningAngle_eta_gg_between_gg_mc = new TH1D("OpeningAngle_eta_gg_between_gg_mc",
1240 "OpeningAngle_eta_gg_between_gg_mc (between #gamma#gamma from "
1241 "#eta); angle [deg];#",
1242 900, -0.1, 89.9);
1244 OpeningAngle_eta_gg_between_gg_reffited = new TH1D("OpeningAngle_eta_gg_between_gg_reffited",
1245 "OpeningAngle_eta_gg_between_gg_reffited (between #gamma#gamma "
1246 "from #eta); angle [deg];#",
1247 900, -0.1, 89.9);
1250 new TH1D("InvMass_eta_gg_allcombinations_mc", "InvMass_eta_gg_allcombinations_mc; invariant mass in GeV/c^{2};#",
1251 200, 0.0, 2.0);
1254 new TH1D("InvMass_eta_gg_allcombinations_reffited",
1255 "InvMass_eta_gg_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1257 EMT_eta_gg = new TH1D("EMT_eta_gg", "EMT_eta_gg; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1258 fHistoList_eta_gg.push_back(EMT_eta_gg);
1260 "InvMass_eta_gg_reco_aftercuts", "InvMass_eta_gg_reco_aftercuts; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1263 new TH2D("rap_vs_pt_eta_gg_reco_aftercuts", "rap_vs_pt_eta_gg_reco_aftercuts; rapidity y; p_{t} in GeV/c ", 10, 0.,
1264 4., 5, 0., 2.);
1267 new TH2D("rap_vs_pt_NOTeta_gg_reco_aftercuts", "rap_vs_pt_NOTeta_gg_reco_aftercuts; rapidity y; p_{t} in GeV/c ",
1268 10, 0., 4., 5, 0., 2.);
1270
1271
1274 new TH1D("InvMass_eta_ppg_mc", "InvMass_eta_ppg_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1277 new TH1D("InvMass_eta_ppg_reffited", "InvMass_eta_ppg_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1280 new TH1D("InvMassPhoton_eta_ppg_mc", "InvMassPhoton_eta_ppg_mc; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1283 "InvMassPhoton_eta_ppg_reffited", "InvMassPhoton_eta_ppg_reffited; invariant mass in GeV/c^{2};#", 60, -0.01, 0.05);
1286 new TH1D("OpeningAnglePhoton_eta_ppg_mc", "OpeningAnglePhoton_eta_ppg_mc (between e+e- from #gamma); angle [deg];#",
1287 100, -0.1, 9.9);
1289 OpeningAnglePhoton_eta_ppg_reffited = new TH1D("OpeningAnglePhoton_eta_ppg_reffited",
1290 "OpeningAnglePhoton_eta_ppg_reffited (between e+e- from #gamma); "
1291 "angle [deg];#",
1292 100, -0.1, 9.9);
1295 new TH1D("InvMass_eta_ppg_allcombinations_mc", "InvMass_eta_ppg_allcombinations_mc; invariant mass in GeV/c^{2};#",
1296 200, 0.0, 2.0);
1299 new TH1D("InvMass_eta_ppg_allcombinations_reffited",
1300 "InvMass_eta_ppg_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1302 Pion_P_fromEta_reco = new TH1D("Pion_P_fromEta_reco", "Pion_P_fromEta_reco; pions P in GeV/c;#", 500, 0.0, 5.0);
1304 Pion_P_elsewhere_reco = new TH1D("Pion_P_elsewhere_reco", "Pion_P_elsewhere_reco; pions P in GeV/c;#", 500, 0.0, 5.0);
1307 new TH1D("Pion_Pt_fromEta_reco", "Pion_Pt_fromEta_reco; pions P_{t} in GeV/c;#", 500, 0.0, 5.0);
1310 new TH1D("Pion_Pt_elsewhere_reco", "Pion_Pt_elsewhere_reco; pions P_{t} in GeV/c;#", 500, 0.0, 5.0);
1312 OA_betweenPions_fromEta_mc = new TH1D("OA_betweenPions_fromEta_mc",
1313 "OA_betweenPions_fromEta_mc (between #p^{+}#p^{-} from #eta); "
1314 "angle [deg];#",
1315 900, -0.1, 89.9);
1317 OA_betweenPions_fromEta_reco = new TH1D("OA_betweenPions_fromEta_reco",
1318 "OA_betweenPions_fromEta_reco (between #p^{+}#p^{-} from #eta); "
1319 "angle [deg];#",
1320 900, -0.1, 89.9);
1322 OA_betweenPions_fromEta_reco_wrongcombinations = new TH1D("OA_betweenPions_fromEta_reco_wrongcombinations",
1323 "OA_betweenPions_fromEta_reco_wrongcombinations (between "
1324 "#p^{+}#p^{-} from #eta); angle [deg];#",
1325 900, -0.1, 89.9);
1327 EMT_eta_ppg = new TH1D("EMT_eta_ppg", "EMT_eta_ppg; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1330 new TH1D("EMT_eta_three_body", "EMT_eta_three_body; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1333 "InvMass_eta_ppg_reco_aftercuts", "InvMass_eta_ppg_reco_aftercuts; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1336 new TH2D("rap_vs_pt_eta_ppg_reco_aftercuts", "rap_vs_pt_eta_ppg_reco_aftercuts; rapidity y; p_{t} in GeV/c ", 10,
1337 0., 4., 5, 0., 2.);
1340 new TH2D("rap_vs_pt_NOTeta_ppg_reco_aftercuts", "rap_vs_pt_NOTeta_ppg_reco_aftercuts; rapidity y; p_{t} in GeV/c ",
1341 10, 0., 4., 5, 0., 2.);
1343
1344
1347 new TH1D("InvMass_eta_ppp_mc", "InvMass_eta_ppp_mc; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1350 new TH1D("InvMass_eta_ppp_reffited", "InvMass_eta_ppp_reffited; invariant mass in GeV/c^{2};#", 50, 0.3, 0.8);
1353 new TH1D("InvMass_eta_Npion_mc", "InvMass_eta_Npion_mc; invariant mass in GeV/c^{2};#", 30, 0.0, 0.3);
1356 new TH1D("InvMass_eta_Npion_reffited", "InvMass_eta_Npion_reffited; invariant mass in GeV/c^{2};#", 30, 0., 0.3);
1359 new TH1D("InvMass_eta_ppp_allcombinations_mc", "InvMass_eta_ppp_allcombinations_mc; invariant mass in GeV/c^{2};#",
1360 200, 0.0, 2.0);
1363 new TH1D("InvMass_eta_ppp_allcombinations_reffited",
1364 "InvMass_eta_ppp_allcombinations_reffited; invariant mass in GeV/c^{2};#", 200, 0.0, 2.0);
1366}
Data class for STS tracks.
int32_t GetStsTrackIndex() const
vector< TH1 * > fHistoList_eta_gg
vector< CbmMCTrack * > El_Photon_Pion_Eta_MCtrack
Double_t CalculateOpeningAngleBetweenGammas_Reco(TVector3 electron1, TVector3 electron2, TVector3 electron3, TVector3 electron4)
Double_t CalculateOpeningAngleBetweenGammas_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2, CbmMCTrack *mctrack3, CbmMCTrack *mctrack4)
std::vector< std::vector< TVector3 > > EMT_ppg_pp_pair_momenta
vector< TVector3 > El_Photon_Eta_refmomentum
vector< CbmMCTrack * > fMCtracks
vector< CbmMCTrack * > All_El_MCtrack
void EtaChargedPionsGammaAnalysis(vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TVector3 > RefMomEl, vector< CbmMCTrack * > MCEl, vector< Int_t > IdEl, vector< TH1 * > ppg)
TClonesArray * fStsTrackMatches
std::vector< int > EMT_ppg_pp_Event
vector< TVector3 > Pion_Eta_refmomentum
void EtaChargedPionsGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut, int Event, vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TVector3 > RefMomEl, vector< CbmMCTrack * > MCEl, vector< Int_t > IdEl, vector< TH1 * > ppg, vector< TH1 * > ppp)
void Exec(int Event, double OpeningAngleCut, double GammaInvMassCut)
std::vector< std::vector< TVector3 > > ECPGA_pions_RefMom
std::vector< std::vector< CbmMCTrack * > > ECPGA_leptons_MC
std::vector< int > EMT_ppg_positive_pion_Event
std::vector< TVector3 > EMT_ppg_positive_pion_momenta
vector< CbmMCTrack * > Pion_Eta_MCtrack
std::vector< std::vector< CbmMCTrack * > > EDGA_MC
std::vector< int > EMT_gg_Event
void EtaPosNegNeutralPionsAnalysis(vector< TVector3 > RefMomNeutral, vector< CbmMCTrack * > MCNeutral, vector< Int_t > IdNeutral, vector< TVector3 > RefMomPion, vector< CbmMCTrack * > MCPion, vector< Int_t > IdPion, vector< TH1 * > ppp)
std::vector< std::vector< int > > EDGA_Id
std::vector< std::vector< TVector3 > > EMT_ppg_ee_pair_momenta
vector< TVector3 > frefmomenta
void EtaDoubleGammaAnalysis(vector< TVector3 > RefMom, vector< CbmMCTrack * > MC, vector< Int_t > Id, vector< TH1 * > gg)
vector< TVector3 > El_Photon_Pion_Eta_refmomentum
std::vector< std::vector< TVector3 > > ECPGA_leptons_RefMom
std::vector< std::vector< TVector3 > > EMT_gg_pair_momenta
TH1D * OA_betweenPions_fromEta_reco_wrongcombinations
std::vector< std::vector< CbmMCTrack * > > ECPGA_pions_MC
vector< TVector3 > All_El_refmomentum
void EtaDoubleGammaAnalysis_plusBG(double OpeningAngleCut, double GammaInvMassCut, int Event, vector< TVector3 > RefMom, vector< CbmMCTrack * > MC, vector< Int_t > Id, vector< TH1 * > gg)
std::vector< TVector3 > EMT_ppg_negative_pion_momenta
vector< CbmMCTrack * > All_Pion_MCtrack
std::vector< std::vector< TVector3 > > EDGA_RefMom
vector< TH1 * > fHistoList_eta_ppg
std::vector< int > EMT_ppg_ee_Event
vector< TH1 * > fHistoList_eta_ppp
vector< CbmMCTrack * > El_Photon_Eta_MCtrack
std::vector< int > EMT_ppg_negative_pion_Event
vector< int > El_Photon_Pion_Eta_Id
vector< TVector3 > All_Pion_refmomentum
std::vector< std::vector< int > > ECPGA_leptons_Id
static TVector3 FitToVertexAndGetChi(CbmStsTrack *stsTrack, double x, double y, double z, double &chi)
static LmvmKinePar CalculateKinematicParams_2el_2pions(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static Double_t CalculateOpeningAngle_Reco(TVector3 electron1, TVector3 electron2)
static double Invmass_2particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2)
static Double_t CalculateOpeningAngle_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2)
static Double_t CalculateOpeningAngleBetweenPions_MC(CbmMCTrack *mctrack1, CbmMCTrack *mctrack2)
calculate opening angle between two pions using MCtrue momenta
static Double_t CalculateOpeningAngleBetweenPions_Reco(TVector3 electron1, TVector3 electron2)
calculate opening angle between two pions using reconstructed momenta
static LmvmKinePar CalculateKinematicParams_4particles(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static double Invmass_6particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4, const CbmMCTrack *mctrack5, const CbmMCTrack *mctrack6)
static double Invmass_2particles_RECO(const TVector3 part1, const TVector3 part2)
static double Invmass_4particles_MC(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
static double Invmass_4particles_RECO(const TVector3 part1, const TVector3 part2, const TVector3 part3, const TVector3 part4)
static double Invmass_4el_2pions_RECO(const TVector3 part1El, const TVector3 part2El, const TVector3 part3El, const TVector3 part4El, const TVector3 part5Pion, const TVector3 part6Pion)
static double Invmass_2el_2pions_RECO(const TVector3 part1El, const TVector3 part2El, const TVector3 part3Pion, const TVector3 part4Pion)
double GetStartZ() const
Definition CbmMCTrack.h:75
int32_t GetMotherId() const
Definition CbmMCTrack.h:69
double GetStartX() const
Definition CbmMCTrack.h:73
int32_t GetPdgCode() const
Definition CbmMCTrack.h:68
void Get4Momentum(TLorentzVector &momentum) const
Definition CbmMCTrack.h:173
double GetStartY() const
Definition CbmMCTrack.h:74
int32_t GetNofLinks() const
Definition CbmMatch.h:42
const CbmLink & GetMatchedLink() const
Definition CbmMatch.h:41
Double_t fMinv
Definition LmvmKinePar.h:22
Double_t fPt
Definition LmvmKinePar.h:20
Double_t fRapidity
Definition LmvmKinePar.h:21
Hash for CbmL1LinkKey.