CbmRoot
CbmAnaConversionKF.cxx
Go to the documentation of this file.
1 
11 #include "CbmAnaConversionKF.h"
12 
13 
14 // included from CbmRoot
15 #include "CbmKFParticleFinder.h"
16 #include "CbmKFParticleFinderQA.h"
17 #include "CbmMCTrack.h"
18 #include "CbmStsTrack.h"
19 #include "CbmTrackMatchNew.h"
20 #include "FairRootManager.h"
21 #include "KFParticleTopoReconstructor.h"
22 #include "KFTopoPerformance.h"
23 #include "TDatabasePDG.h"
24 
25 
26 #define M2E 2.6112004954086e-7
27 
28 using namespace std;
29 
30 
32  : fKFMcParticles(NULL)
33  , fMcTracks(NULL)
34  , fStsTracks(NULL)
35  , fStsTrackMatches(NULL)
36  , fKFparticle(NULL)
37  , fKFparticleFinderQA(NULL)
38  , fKFtopo(NULL)
39  , fKFtopoPerf(0)
40  , trackindexarray()
41  , particlecounter(0)
42  , particlecounter_2daughters(0)
43  , particlecounter_3daughters(0)
44  , particlecounter_4daughters(0)
45  , particlecounter_all(0)
46  , fhPi0_NDaughters(NULL)
47  , fNofGeneratedPi0_allEvents(0)
48  , fNofPi0_kfparticle_allEvents(0)
49  , fNofGeneratedPi0(0)
50  , fNofPi0_kfparticle(0)
51  , fhPi0Ratio(NULL)
52  , fhPi0_mass(NULL)
53  , fSignalIds()
54  , fGhostIds()
55  , fHistoList_kfparticle()
56  , particlevector()
57  , particleMatch()
58  , electronIDs()
59  , gammaIDs()
60  , fhInvMassPi0WithFullReco(NULL)
61  , fhInvMass2Gammas(NULL)
62  , fhInvMass2Gammas_cut(NULL)
63  , fhKF_particlevector(NULL)
64  , fhKF_trackvector(NULL)
65  , fhKF_NofPi0(NULL)
66  , fhKF_NofPi0_signal(NULL)
67  , fhKF_NofPi0_trackvector(NULL)
68  , fKF_photon_pairs()
69  , fhKF_invmass_fullReco(NULL)
70  , timer()
71  , fTime(0.) {}
72 
74 
75 
77  FairRootManager* ioman = FairRootManager::Instance();
78  if (NULL == ioman) {
79  Fatal("CbmAnaConversionKF::Init", "RootManager not instantised!");
80  }
81 
82  fKFMcParticles = (TClonesArray*) ioman->GetObject("KFMCParticles");
83  if (NULL == fKFMcParticles) {
84  Fatal("CbmAnaConversionKF::Init", "No KFMCParticles array!");
85  }
86 
87  fMcTracks = (TClonesArray*) ioman->GetObject("MCTrack");
88  if (NULL == fMcTracks) {
89  Fatal("CbmAnaConversionKF::Init", "No MCTrack array!");
90  }
91 
92  fStsTracks = (TClonesArray*) ioman->GetObject("StsTrack");
93  if (NULL == fStsTracks) {
94  Fatal("CbmAnaConversionKF::Init", "No StsTrack array!");
95  }
96 
97  fStsTrackMatches = (TClonesArray*) ioman->GetObject("StsTrackMatch");
98  if (NULL == fStsTrackMatches) {
99  Fatal("CbmAnaConversionKF::Init", "No StsTrackMatch array!");
100  }
101 
102 
104  fKFtopoPerf = new KFTopoPerformance;
105  fKFtopoPerf->SetTopoReconstructor(fKFtopo);
106 
107  InitHistos();
108 }
109 
110 
112  fHistoList_kfparticle.clear();
113 
114  // #############################################
115  // Histograms related to KFParticle results
116  fhPi0_NDaughters = new TH1D(
117  "fhPi0_NDaughters", "fhPi0_NDaughters;number of daughers;#", 4, 0.5, 4.5);
118  fhPi0Ratio = new TH1D("fhPi0Ratio", "fhPi0Ratio; ratio;#", 1000, 0., 0.02);
119  fhPi0_mass = new TH1D("fhPi0_mass", "fhPi0_mass;mass;#", 500, 0., 0.5);
123 
124  fhInvMassPi0WithFullReco = new TH1D("fhInvMassPi0WithFullReco",
125  "fhInvMassPi0WithFullReco;invmass;#",
126  400,
127  0.,
128  2.);
130 
132  new TH1D("fhInvMass2Gammas", "fhInvMass2Gammas;invmass;#", 400, 0., 2.);
134  fhInvMass2Gammas_cut = new TH1D(
135  "fhInvMass2Gammas_cut", "fhInvMass2Gammas_cut;invmass;#", 400, 0., 2.);
137 
139  new TH1D("fhKF_particlevector", "fhKF_particlevector;;#", 10, 0., 10.);
141  fhKF_particlevector->GetXaxis()->SetBinLabel(1, "electrons");
142  fhKF_particlevector->GetXaxis()->SetBinLabel(2, "gamma");
143  fhKF_particlevector->GetXaxis()->SetBinLabel(3, "pi0");
145  new TH1D("fhKF_trackvector", "fhKF_trackvector;;#", 10, 0., 10.);
147  fhKF_trackvector->GetXaxis()->SetBinLabel(1, "electrons");
148  fhKF_trackvector->GetXaxis()->SetBinLabel(2, "gamma");
149  fhKF_trackvector->GetXaxis()->SetBinLabel(3, "pi0");
150 
151  fhKF_NofPi0 = new TH1D("fhKF_NofPi0", "fhKF_NofPi0;nof;#", 10, -0.5, 9.5);
154  new TH1D("fhKF_NofPi0_signal", "fhKF_NofPi0_signal;nof;#", 10, -0.5, 9.5);
156  fhKF_NofPi0_trackvector = new TH1D(
157  "fhKF_NofPi0_trackvector", "fhKF_NofPi0_trackvector;nof;#", 10, -0.5, 9.5);
159 
160 
161  fhKF_invmass_fullReco = new TH1D(
162  "fhKF_invmass_fullReco", "fhKF_invmass_fullReco;invmass;#", 400, 0., 2.);
164 }
165 
166 
168  //gDirectory->cd("analysis-conversion");
169  gDirectory->mkdir("KF");
170  gDirectory->cd("KF");
171  for (UInt_t i = 0; i < fHistoList_kfparticle.size(); i++) {
172  fHistoList_kfparticle[i]->Write();
173  }
174  gDirectory->cd("..");
175 
176  cout << "CbmAnaConversionKF: Realtime - " << fTime << endl;
177 }
178 
179 
181  timer.Start();
182 
183  electronIDs.clear();
184  gammaIDs.clear();
185  fKF_photon_pairs.clear();
186 
187  //KFParticle_Analysis();
188  test2();
189 
191  CombinePhotons();
192 
193  Reconstruct();
195 
196  timer.Stop();
197  fTime += timer.RealTime();
198 }
199 
200 
202  CbmKFParticleFinderQA* kfparticleQA) {
203  fKFparticle = kfparticle;
204  fKFparticleFinderQA = kfparticleQA;
205  if (fKFparticle) {
206  cout << "CbmAnaConversionKF: kf works" << endl;
207  } else {
208  cout << "CbmAnaConversionKF: kf works not" << endl;
209  }
210 }
211 
212 
213 void CbmAnaConversionKF::SetSignalIds(std::vector<int>* signalids) {
214  fSignalIds.clear();
215  fSignalIds = *signalids;
216 }
217 
218 void CbmAnaConversionKF::SetGhostIds(std::vector<int>* ghostids) {
219  fGhostIds.clear();
220  fGhostIds = *ghostids;
221 }
222 
223 
224 /*
225 void CbmAnaConversionKF::KFParticle_Analysis()
226 {
227  timer.Start();
228 
229  int testkf = fKFtopo->NPrimaryVertices();
230  cout << "KFParticle_Analysis - test kf NPrimaryVertices: " << testkf << endl;
231 
232  const KFPTrackVector* kftrackvector;
233  kftrackvector = fKFtopo->GetTracks();
234  cout << "KFParticle_Analysis - trackvector: " << kftrackvector->Size() << endl;
235  KFParticle testparticle;
236  KFPTrack testtrack;
237  const int bla = 2;
238  kftrackvector->GetTrack(testtrack, bla);
239  cout << "KFParticle_Analysis - test p: " << testtrack.GetP() << endl;
240  kfvector_int pdgvector;
241  pdgvector = kftrackvector->PDG();
242  //cout << "pdg tests: " << pdgvector[1] << "\t" << pdgvector[2] << "\t" << pdgvector[3] << "\t" << pdgvector[4] << endl;
243 
244 
245  int nofpions = 0;
246  //nofpions = kftrackvector->NPions();
247  cout << "KFParticle_Analysis - number of pions: " << nofpions << endl;
248 
249  int testpi = 0;
250  for(int i=0; i<kftrackvector->Size(); i++) {
251  if(TMath::Abs(pdgvector[i]) == 211) testpi++;
252  }
253  cout << "KFParticle_Analysis - testpi: " << testpi << endl;
254 
255  //kftrackvector->PrintTracks();
256 
257 
258  vector<KFParticle> particlevector;
259  particlevector = fKFtopo->GetParticles();
260 
261 
262  cout << "KFParticle_Analysis - particlevector size: " << particlevector.size() << "\t";
263  for(int i=0; i<particlevector.size(); i++) {
264  if(particlevector[i].NDaughters() > 3) {
265  cout << particlevector[i].NDaughters() << "\t";
266  }
267  if(particlevector[i].GetPDG() == 111) {
268  cout << "blubb" << particlevector[i].NDaughters() << "\t";
269  if(particlevector[i].NDaughters() > 2) {
270  particlecounter++;
271  }
272 
273  std::vector<int> daughterid;
274  daughterid = particlevector[i].DaughterIds();
275  cout << "daughterids: ";
276  for(int j=0; j<daughterid.size(); j++) {
277  cout << daughterid[j] << "/";
278  cout << "pdg-" << particlevector[daughterid[j]].GetPDG() << "-";
279  }
280  particlecounter_all++;
281  }
282  }
283  cout << endl;
284 
285 
286 
287  cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
288  cout << "KFParticle_Analysis - SignalIds size: " << fSignalIds.size() << "\t";
289  for(int i=0; i<fSignalIds.size(); i++) {
290  Int_t pdg = (fKFtopo->GetParticles()[fSignalIds[i]]).GetPDG();
291  Int_t daughters = (fKFtopo->GetParticles()[fSignalIds[i]]).NDaughters();
292  vector<int> daughterIds = (fKFtopo->GetParticles()[fSignalIds[i]]).DaughterIds();
293  float mass = 0;
294  float mass_sigma = 0;
295  (fKFtopo->GetParticles()[fSignalIds[i]]).GetMass(mass, mass_sigma);
296  if(pdg == 111) {
297  cout << "CbmAnaConversionKF: pi0 found, signal id: " << (fKFtopo->GetParticles()[fSignalIds[i]]).Id() << "\t daughters: " << daughters << "\t";
298  cout << "daughter ids: ";
299  Int_t electronids[4] = {0};
300  for(int j=0; j<daughterIds.size(); j++) {
301  vector<int> granddaughterIds = (fKFtopo->GetParticles()[daughterIds[j]]).DaughterIds();
302  cout << daughterIds[j] << " (" << granddaughterIds[0] << ",pdg" << (fKFtopo->GetParticles()[granddaughterIds[0]]).GetPDG() << "/" << granddaughterIds[1] << ",pdg" << (fKFtopo->GetParticles()[granddaughterIds[1]]).GetPDG() << ")" << " / ";
303  electronids[j*2] = granddaughterIds[0];
304  electronids[j*2+1] = granddaughterIds[1];
305  }
306  cout << endl;
307 
308  cout << "the 4 electrons and their grandmotherids: ";
309  CbmMCTrack * electrontracks[4];
310  for(int k=0; k<4; k++) {
311  electrontracks[k] = (CbmMCTrack*) fMcTracks->At(electronids[k]);
312  CbmMCTrack* mother = (CbmMCTrack*) fMcTracks->At(electrontracks[k]->GetMotherId());
313  if(mother == NULL) continue;
314  CbmMCTrack* grandmother = (CbmMCTrack*) fMcTracks->At(mother->GetMotherId());
315  if(grandmother == NULL) continue;
316  cout << "gmid" << mother->GetMotherId() << "pdg" << grandmother->GetPdgCode() << " / ";
317 
318  }
319  cout << endl;
320 
321  particlecounter++;
322  if(daughters == 2) { particlecounter_2daughters++; }
323  if(daughters == 3) { particlecounter_3daughters++; }
324  if(daughters == 4) { particlecounter_4daughters++; }
325  fhPi0_NDaughters->Fill(daughters);
326  fhPi0_mass->Fill(mass);
327  fNofPi0_kfparticle++;
328  }
329  }
330 
331  cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
332  cout << "KFParticle_Analysis - GhostIds size: " << fGhostIds.size() << "\t";
333  for(int i=0; i<fGhostIds.size(); i++) {
334  Int_t pdg = (fKFtopo->GetParticles()[fGhostIds[i]]).GetPDG();
335  Int_t daughters = (fKFtopo->GetParticles()[fGhostIds[i]]).NDaughters();
336  if(pdg == 111) {
337  cout << "pi0 found!\t";
338  particlecounter++;
339  if(daughters == 2) { particlecounter_2daughters++; }
340  if(daughters == 3) { particlecounter_3daughters++; }
341  if(daughters == 4) { particlecounter_4daughters++; }
342  //fhPi0_NDaughters->Fill(daughters);
343  }
344  }
345 
346  cout << endl;
347  cout << "%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%" << endl;
348 
349 
350 
351 
352 // trackindexarray.clear();
353 // trackindexarray = fKFtopo->GetPVTrackIndexArray();
354 // cout << "trackindexarray: " << trackindexarray.size() << endl;
355 
356 // KFVertex vertex;
357 // vertex = fKFtopo->GetPrimKFVertex();
358 // int testnof = 0;
359 // testnof = vertex.GetNContributors();
360 // cout << "nof contributors: " << testnof << endl;
361 
362 // KFParticle particle;
363 // particle = fKFtopo->GetPrimVertex();
364 
365 
366 // for(int i=0; i < kftrackvector->Size(); i++) {
367  //const kfvector_int& tempPDG;
368  //tempPDG = kftrackvector->PDG;
369  //fKFVertex->At(trackindexarray[i]);
370 
371 // }
372 
373 
374 
375  test();
376 
377  timer.Stop();
378  fTime += timer.RealTime();
379 }
380 */
381 
382 
384  int nofparticles = fKFMcParticles->GetEntriesFast();
385  cout << "CbmAnaConversionKF: test nof " << nofparticles << endl;
386  for (int i = 0; i < nofparticles; i++) {
387  CbmMCTrack* mcTrack1 = (CbmMCTrack*) fKFMcParticles->At(i);
388  int pdg = TMath::Abs(mcTrack1->GetPdgCode());
389  if (pdg == 111) { cout << "CbmAnaConversionKF: test successful!" << endl; }
390  }
391 
392 
393  vector<vector<int>> ids;
394  const vector<KFParticle>& particles =
395  fKFparticle->GetTopoReconstructor()->GetParticles();
396  for (unsigned int iPart = 0; iPart < fSignalIds.size(); iPart++) {
397  if (particles[fSignalIds[iPart]].GetPDG() != 111) continue;
398  //some cuts on pi0 if needed
399 
400  const KFParticle& pi0 = particles[fSignalIds[iPart]];
401  vector<int> electrons;
402  for (int iGamma = 0; iGamma < pi0.NDaughters(); iGamma++) {
403  const int GammaID = pi0.DaughterIds()[iGamma];
404  const KFParticle& Gamma = particles[GammaID];
405  for (int iElectron = 0; iElectron < Gamma.NDaughters(); iElectron++) {
406  int ElectronID = Gamma.DaughterIds()[iElectron];
407  const KFParticle& Electron = particles[ElectronID];
408  int STStrackID = Electron.DaughterIds()[0];
409  electrons.push_back(STStrackID);
410  }
411  }
412  ids.push_back(electrons);
413  }
414 
415  if (ids.size() > 0) {
416  cout << "NEW TEST: (sts ids) ";
417  for (unsigned int i = 0; i < ids.size(); i++) {
418  for (int j = 0; j < 4; j++) {
419  cout << " " << ids[i][j];
420  }
421  cout << " | ";
422  }
423  cout << endl;
424 
425 
426  cout << "MC-pdgs: ";
427  CbmMCTrack* mcTracks[4];
428  for (unsigned int i = 0; i < ids.size(); i++) {
429  for (int j = 0; j < 4; j++) {
430  CbmStsTrack* stsTrack = (CbmStsTrack*) fStsTracks->At(ids[i][j]);
431  if (stsTrack == NULL) return;
432  CbmTrackMatchNew* stsMatch =
433  (CbmTrackMatchNew*) fStsTrackMatches->At(ids[i][j]);
434  if (stsMatch == NULL) return;
435  int stsMcTrackId = stsMatch->GetMatchedLink().GetIndex();
436  if (stsMcTrackId < 0) return;
437  mcTracks[j] = (CbmMCTrack*) fMcTracks->At(stsMcTrackId);
438  if (mcTracks[j] == NULL) return;
439 
440  CbmMCTrack* mother =
441  (CbmMCTrack*) fMcTracks->At(mcTracks[j]->GetMotherId());
442  CbmMCTrack* grandmother =
443  (CbmMCTrack*) fMcTracks->At(mother->GetMotherId());
444 
445 
446  cout << " " << stsMcTrackId << "/" << mcTracks[j]->GetPdgCode()
447  << "(motherid" << mcTracks[j]->GetMotherId() << ",motherpdg"
448  << mother->GetPdgCode() << ",grandmotherpdg"
449  << grandmother->GetPdgCode() << ",grandmotherid"
450  << mother->GetMotherId() << ")";
451  }
452  }
453  cout << endl;
454  Double_t mass =
456  cout << "mass: " << mass << endl;
457  }
458 }
459 
460 
462  const CbmMCTrack* mctrack2,
463  const CbmMCTrack* mctrack3,
464  const CbmMCTrack* mctrack4)
465 // calculation of invariant mass from four electrons/positrons
466 {
467  TLorentzVector lorVec1;
468  mctrack1->Get4Momentum(lorVec1);
469 
470  TLorentzVector lorVec2;
471  mctrack2->Get4Momentum(lorVec2);
472 
473  TLorentzVector lorVec3;
474  mctrack3->Get4Momentum(lorVec3);
475 
476  TLorentzVector lorVec4;
477  mctrack4->Get4Momentum(lorVec4);
478 
479 
480  TLorentzVector sum;
481  sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
482  cout << "mc: \t" << sum.Px() << " / " << sum.Py() << " / " << sum.Pz()
483  << " / " << sum.E() << "\t => mag = " << sum.Mag() << endl;
484 
485  return sum.Mag();
486 }
487 
488 
490  Int_t pi0counter_trackvector = 0;
491  const KFPTrackVector* kftrackvector;
492  kftrackvector = fKFtopo->GetTracks();
493  Int_t kftv_size = kftrackvector->Size();
494  cout << "CbmAnaConversionKF: size of kftrackvector: " << kftv_size << endl;
495 
496  kfvector_int vector_pdgs = kftrackvector->PDG();
497 
498  cout << "CbmAnaConversionKF: pdgs in trackvector: ";
499  for (int i = 0; i < kftv_size; i++) {
500  cout << vector_pdgs[i] << " / ";
501  if (TMath::Abs(vector_pdgs[i]) == 111) {
502  pi0counter_trackvector++;
503  fhKF_trackvector->Fill(2);
504  }
505  if (TMath::Abs(vector_pdgs[i]) == 11) { fhKF_trackvector->Fill(0); }
506  if (TMath::Abs(vector_pdgs[i]) == 22) { fhKF_trackvector->Fill(1); }
507  }
508  cout << endl;
509 
510  fhKF_NofPi0_trackvector->Fill(pi0counter_trackvector);
511 
512 
513  // particlevector to get all pi0 detected by KFParticlePackage (including signal, ghost, background); trackvector DOES NOT WORK
514  particlevector.clear();
515  particleMatch.clear();
516  //vector<KFParticle> particlevector;
517 
518 
519  particlevector = fKFtopo->GetParticles();
520 
521 
522  particleMatch = fKFtopoPerf->ParticlesMatch();
523  Int_t pv_size = particlevector.size();
524  cout << "CbmAnaConversionKF: size of particlevector: " << pv_size << endl;
525  cout << "CbmAnaConversionKF: size of particleMatch: " << particleMatch.size()
526  << endl;
527 
528  Int_t electroncounter = 0;
529  Int_t pi0counter = 0;
530  Int_t pi0counter_signal = 0;
531  Int_t gammacounter = 0;
532  for (int i = 0; i < pv_size; i++) {
533  if (TMath::Abs(particlevector[i].GetPDG()) == 11) {
534  electroncounter++;
535  electronIDs.push_back(i);
536  fhKF_particlevector->Fill(0);
537  }
538  if (TMath::Abs(particlevector[i].GetPDG()) == 111) {
539  pi0counter++;
540  fhKF_particlevector->Fill(2);
541  //if(particleMatch[i].IsMatchedWithPdg()) {
542  // pi0counter_signal++;
543  //}
544  }
545  if (TMath::Abs(particlevector[i].GetPDG()) == 22) {
546  gammacounter++;
547  gammaIDs.push_back(i);
548  fhKF_particlevector->Fill(1);
549  }
550  }
551 
552 
553  fhKF_NofPi0->Fill(pi0counter);
554  fhKF_NofPi0_signal->Fill(pi0counter_signal);
555 
556  cout << "CbmAnaConversionKF: nof electrons in particlevector: "
557  << electroncounter << endl;
558  cout << "CbmAnaConversionKF: nof pi0 in particlevector: " << pi0counter
559  << endl;
560  cout << "CbmAnaConversionKF: nof gamma in particlevector: " << gammacounter
561  << endl;
562 }
563 
564 
566  Int_t nof = electronIDs.size();
567  if (nof >= 2) {
568  for (int a = 0; a < nof - 1; a++) {
569  for (int b = a + 1; b < nof; b++) {
570  KFParticle electron1 = particlevector[electronIDs[a]];
571  KFParticle electron2 = particlevector[electronIDs[b]];
572  Double_t openingAngle = electron1.GetAngle(electron2);
573 
574  TLorentzVector lorentzE1;
575  lorentzE1.SetPxPyPzE(electron1.GetPx(),
576  electron1.GetPy(),
577  electron1.GetPz(),
578  electron1.GetE());
579  TLorentzVector lorentzE2;
580  lorentzE2.SetPxPyPzE(electron2.GetPx(),
581  electron2.GetPy(),
582  electron2.GetPz(),
583  electron2.GetE());
584  TLorentzVector sum = lorentzE1 + lorentzE2;
585  Double_t invmass = sum.Mag();
586 
587  if (openingAngle < 1 && invmass < 0.03) {
588  vector<int> pair; // = {a, b};
589  pair.push_back(a);
590  pair.push_back(b);
591  fKF_photon_pairs.push_back(pair);
592  }
593  }
594  }
595  }
596 }
597 
598 
600  Int_t nof = fKF_photon_pairs.size();
601  if (nof >= 2) {
602  for (int a = 0; a < nof - 1; a++) {
603  for (int b = a + 1; b < nof; b++) {
604  Int_t electron11 = fKF_photon_pairs[a][0];
605  Int_t electron12 = fKF_photon_pairs[a][1];
606  Int_t electron21 = fKF_photon_pairs[b][0];
607  Int_t electron22 = fKF_photon_pairs[b][1];
608 
609  Double_t invmass = Invmass_4particlesRECO(particlevector[electron11],
610  particlevector[electron12],
611  particlevector[electron21],
612  particlevector[electron22]);
613  fhKF_invmass_fullReco->Fill(invmass);
614  }
615  }
616  }
617 }
618 
619 
621  Int_t nof = electronIDs.size();
622  if (nof >= 4) {
623  for (int a = 0; a < nof - 3; a++) {
624  for (int b = a + 1; b < nof - 2; b++) {
625  for (int c = b + 1; c < nof - 1; c++) {
626  for (int d = c + 1; d < nof; d++) {
627  Int_t check1 = (particlevector[electronIDs[a]].GetPDG() > 0);
628  Int_t check2 = (particlevector[electronIDs[b]].GetPDG() > 0);
629  Int_t check3 = (particlevector[electronIDs[c]].GetPDG() > 0);
630  Int_t check4 = (particlevector[electronIDs[d]].GetPDG() > 0);
631  Int_t test1234 = check1 + check2 + check3 + check4;
632  if (test1234 != 2)
633  continue; // need two electrons and two positrons
634 
635 
636  KFParticle particle1 = particlevector[electronIDs[a]];
637  KFParticle particle2 = particlevector[electronIDs[b]];
638  KFParticle particle3 = particlevector[electronIDs[c]];
639  KFParticle particle4 = particlevector[electronIDs[d]];
640  Double_t invmass = Invmass_4particlesRECO(
641  particle1, particle2, particle3, particle4);
642 
643 
644  Bool_t fill = false;
645  Double_t invmass_cut = 0.02;
646 
647  if (particle1.GetZ() == particle2.GetZ()
648  && particle3.GetZ() == particle4.GetZ()) {
649  Double_t invmass_12 = Invmass_2electrons(particle1, particle2);
650  Double_t invmass_34 = Invmass_2electrons(particle3, particle4);
651  if (invmass_12 < invmass_cut && invmass_34 < invmass_cut) {
652  fill = true;
653  }
654  }
655  if (particle1.GetZ() == particle3.GetZ()
656  && particle2.GetZ() == particle4.GetZ()) {
657  Double_t invmass_13 = Invmass_2electrons(particle1, particle3);
658  Double_t invmass_24 = Invmass_2electrons(particle2, particle4);
659  if (invmass_13 < invmass_cut && invmass_24 < invmass_cut) {
660  fill = true;
661  }
662  }
663  if (particle1.GetZ() == particle4.GetZ()
664  && particle2.GetZ() == particle3.GetZ()) {
665  Double_t invmass_14 = Invmass_2electrons(particle1, particle4);
666  Double_t invmass_23 = Invmass_2electrons(particle2, particle3);
667  if (invmass_14 < invmass_cut && invmass_23 < invmass_cut) {
668  fill = true;
669  }
670  }
671 
672 
673  if (fill) fhInvMassPi0WithFullReco->Fill(invmass);
674 
675  /*
676  CbmLmvmKinematicParams params1 = CalculateKinematicParamsReco(fElectrons_momenta[a], fElectrons_momenta[b]);
677  CbmLmvmKinematicParams params2 = CalculateKinematicParamsReco(fElectrons_momenta[a], fElectrons_momenta[c]);
678  CbmLmvmKinematicParams params3 = CalculateKinematicParamsReco(fElectrons_momenta[a], fElectrons_momenta[d]);
679  CbmLmvmKinematicParams params4 = CalculateKinematicParamsReco(fElectrons_momenta[b], fElectrons_momenta[c]);
680  CbmLmvmKinematicParams params5 = CalculateKinematicParamsReco(fElectrons_momenta[b], fElectrons_momenta[d]);
681  CbmLmvmKinematicParams params6 = CalculateKinematicParamsReco(fElectrons_momenta[c], fElectrons_momenta[d]);
682 
683  Double_t openingAngleCut = 1;
684  Int_t IsPhoton_openingAngle1 = (params1.fAngle < openingAngleCut);
685  Int_t IsPhoton_openingAngle2 = (params2.fAngle < openingAngleCut);
686  Int_t IsPhoton_openingAngle3 = (params3.fAngle < openingAngleCut);
687  Int_t IsPhoton_openingAngle4 = (params4.fAngle < openingAngleCut);
688  Int_t IsPhoton_openingAngle5 = (params5.fAngle < openingAngleCut);
689  Int_t IsPhoton_openingAngle6 = (params6.fAngle < openingAngleCut);
690 
691  Double_t invMassCut = 0.03;
692  Int_t IsPhoton_invMass1 = (params1.fMinv < invMassCut);
693  Int_t IsPhoton_invMass2 = (params2.fMinv < invMassCut);
694  Int_t IsPhoton_invMass3 = (params3.fMinv < invMassCut);
695  Int_t IsPhoton_invMass4 = (params4.fMinv < invMassCut);
696  Int_t IsPhoton_invMass5 = (params5.fMinv < invMassCut);
697  Int_t IsPhoton_invMass6 = (params6.fMinv < invMassCut);
698 
699  if(IsPhoton_openingAngle1 && IsPhoton_openingAngle6 && IsPhoton_invMass1 && IsPhoton_invMass6 && (check1 + check2 == 1) && (check3 + check4 == 1)) {
700  fhElectrons_invmass_cut->Fill(invmass);
701  }
702  if(IsPhoton_openingAngle2 && IsPhoton_openingAngle5 && IsPhoton_invMass2 && IsPhoton_invMass5 && (check1 + check3 == 1) && (check2 + check4 == 1)) {
703  fhElectrons_invmass_cut->Fill(invmass);
704  }
705  if(IsPhoton_openingAngle3 && IsPhoton_openingAngle4 && IsPhoton_invMass3 && IsPhoton_invMass4 && (check1 + check4 == 1) && (check2 + check3 == 1)) {
706  fhElectrons_invmass_cut->Fill(invmass);
707  }
708  */
709  }
710  }
711  }
712  }
713  }
714 }
715 
716 
718  Int_t nof = gammaIDs.size();
719  if (nof >= 2) {
720  for (int a = 0; a < nof - 1; a++) {
721  for (int b = a + 1; b < nof; b++) {
722  KFParticle particle1 = particlevector[gammaIDs[a]];
723  KFParticle particle2 = particlevector[gammaIDs[b]];
724  Double_t invmass = Invmass_2gamma(particle1, particle2);
725  Double_t openingAngle =
726  OpeningAngleBetweenPhotons(particle1, particle2);
727 
728  fhInvMass2Gammas->Fill(invmass);
729 
730  if ((TMath::Abs(particle1.GetZ() - particle2.GetZ()) < 0.005)
731  && (openingAngle < 5)) {
732  fhInvMass2Gammas_cut->Fill(invmass);
733  }
734  }
735  }
736  }
737 }
738 
739 
741  KFParticle part2,
742  KFParticle part3,
743  KFParticle part4)
744 // calculation of invariant mass from four electrons/positrons
745 {
746  TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
747  Double_t energy1 = TMath::Sqrt(momentum1.Mag2() + M2E);
748  TLorentzVector lorVec1(momentum1, energy1);
749 
750  TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
751  Double_t energy2 = TMath::Sqrt(momentum2.Mag2() + M2E);
752  TLorentzVector lorVec2(momentum2, energy2);
753 
754  TVector3 momentum3(part3.GetPx(), part3.GetPy(), part3.GetPz());
755  Double_t energy3 = TMath::Sqrt(momentum3.Mag2() + M2E);
756  TLorentzVector lorVec3(momentum3, energy3);
757 
758  TVector3 momentum4(part4.GetPx(), part4.GetPy(), part4.GetPz());
759  Double_t energy4 = TMath::Sqrt(momentum4.Mag2() + M2E);
760  TLorentzVector lorVec4(momentum4, energy4);
761 
762 
763  TLorentzVector sum;
764  sum = lorVec1 + lorVec2 + lorVec3 + lorVec4;
765 
766  return sum.Mag();
767 }
768 
769 
770 Double_t CbmAnaConversionKF::Invmass_2gamma(KFParticle part1,
771  KFParticle part2) {
772  TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
773  Double_t energy1 = TMath::Sqrt(momentum1.Mag2());
774  TLorentzVector lorVec1(momentum1, energy1);
775 
776  TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
777  Double_t energy2 = TMath::Sqrt(momentum2.Mag2());
778  TLorentzVector lorVec2(momentum2, energy2);
779 
780 
781  TLorentzVector sum;
782  sum = lorVec1 + lorVec2;
783 
784  return sum.Mag();
785 }
786 
787 
788 Double_t CbmAnaConversionKF::Invmass_2electrons(KFParticle part1,
789  KFParticle part2) {
790  TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
791  Double_t energy1 = TMath::Sqrt(momentum1.Mag2() + M2E);
792  TLorentzVector lorVec1(momentum1, energy1);
793 
794  TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
795  Double_t energy2 = TMath::Sqrt(momentum2.Mag2() + M2E);
796  TLorentzVector lorVec2(momentum2, energy2);
797 
798 
799  TLorentzVector sum;
800  sum = lorVec1 + lorVec2;
801 
802  return sum.Mag();
803 }
804 
805 
807  KFParticle part2) {
808  TVector3 momentum1(part1.GetPx(), part1.GetPy(), part1.GetPz());
809  Double_t energy1 = TMath::Sqrt(momentum1.Mag2());
810  TLorentzVector lorVec1(momentum1, energy1);
811 
812  TVector3 momentum2(part2.GetPx(), part2.GetPy(), part2.GetPz());
813  Double_t energy2 = TMath::Sqrt(momentum2.Mag2());
814  TLorentzVector lorVec2(momentum2, energy2);
815 
816 
817  Double_t angleBetweenPhotons = lorVec1.Angle(lorVec2.Vect());
818  Double_t theta = 180. * angleBetweenPhotons / TMath::Pi();
819 
820  return theta;
821 }
CbmAnaConversionKF::OpeningAngleBetweenPhotons
Double_t OpeningAngleBetweenPhotons(KFParticle part1, KFParticle part2)
Definition: CbmAnaConversionKF.cxx:806
CbmMCTrack::GetMotherId
Int_t GetMotherId() const
Definition: CbmMCTrack.h:71
CbmMatch::GetMatchedLink
const CbmLink & GetMatchedLink() const
Definition: CbmMatch.h:37
CbmAnaConversionKF::fSignalIds
std::vector< int > fSignalIds
Definition: CbmAnaConversionKF.h:101
CbmAnaConversionKF::fhKF_NofPi0
TH1D * fhKF_NofPi0
Definition: CbmAnaConversionKF.h:121
CbmAnaConversionKF::test
void test()
Definition: CbmAnaConversionKF.cxx:383
CbmKFParticleFinder.h
CbmKFParticleFinder
Definition: CbmKFParticleFinder.h:26
CbmAnaConversionKF::particlevector
std::vector< KFParticle > particlevector
Definition: CbmAnaConversionKF.h:109
CbmKFParticleFinderQA
Definition: CbmKFParticleFinderQA.h:21
CbmMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: CbmMCTrack.h:70
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmAnaConversionKF::Invmass_4particles
Double_t Invmass_4particles(const CbmMCTrack *mctrack1, const CbmMCTrack *mctrack2, const CbmMCTrack *mctrack3, const CbmMCTrack *mctrack4)
Definition: CbmAnaConversionKF.cxx:461
CbmAnaConversionKF::fhPi0Ratio
TH1D * fhPi0Ratio
Definition: CbmAnaConversionKF.h:97
CbmAnaConversionKF::fKF_photon_pairs
std::vector< std::vector< int > > fKF_photon_pairs
Definition: CbmAnaConversionKF.h:126
CbmAnaConversionKF::fKFparticle
CbmKFParticleFinder * fKFparticle
Definition: CbmAnaConversionKF.h:77
CbmAnaConversionKF::fTime
Double_t fTime
Definition: CbmAnaConversionKF.h:132
CbmAnaConversionKF::Invmass_4particlesRECO
Double_t Invmass_4particlesRECO(KFParticle part1, KFParticle part2, KFParticle part3, KFParticle part4)
Definition: CbmAnaConversionKF.cxx:740
CbmAnaConversionKF::InitHistos
void InitHistos()
Definition: CbmAnaConversionKF.cxx:111
CbmAnaConversionKF::Finish
void Finish()
Definition: CbmAnaConversionKF.cxx:167
CbmAnaConversionKF.h
CbmAnaConversionKF::fhInvMassPi0WithFullReco
TH1D * fhInvMassPi0WithFullReco
Definition: CbmAnaConversionKF.h:113
CbmAnaConversionKF::fhKF_NofPi0_trackvector
TH1D * fhKF_NofPi0_trackvector
Definition: CbmAnaConversionKF.h:123
CbmAnaConversionKF::Init
void Init()
Definition: CbmAnaConversionKF.cxx:76
CbmAnaConversionKF::fStsTracks
TClonesArray * fStsTracks
Definition: CbmAnaConversionKF.h:74
CbmStsTrack.h
Data class for STS tracks.
CbmAnaConversionKF::Exec
void Exec()
Definition: CbmAnaConversionKF.cxx:180
CbmAnaConversionKF::gammaIDs
std::vector< int > gammaIDs
Definition: CbmAnaConversionKF.h:112
CbmAnaConversionKF::fhKF_trackvector
TH1D * fhKF_trackvector
Definition: CbmAnaConversionKF.h:119
CbmAnaConversionKF::fKFtopo
const KFParticleTopoReconstructor * fKFtopo
Definition: CbmAnaConversionKF.h:80
CbmAnaConversionKF::CombinePhotons
void CombinePhotons()
Definition: CbmAnaConversionKF.cxx:599
d
double d
Definition: P4_F64vec2.h:24
mcTracks
static vector< vector< QAMCTrack > > mcTracks
Definition: CbmTofHitFinderTBQA.cxx:112
CbmAnaConversionKF::SetSignalIds
void SetSignalIds(std::vector< int > *signalids)
Definition: CbmAnaConversionKF.cxx:213
CbmTrackMatchNew.h
CbmAnaConversionKF::fhInvMass2Gammas_cut
TH1D * fhInvMass2Gammas_cut
Definition: CbmAnaConversionKF.h:116
CbmAnaConversionKF::test2
void test2()
Definition: CbmAnaConversionKF.cxx:489
CbmAnaConversionKF::fhInvMass2Gammas
TH1D * fhInvMass2Gammas
Definition: CbmAnaConversionKF.h:115
CbmAnaConversionKF::SetKF
void SetKF(CbmKFParticleFinder *kfparticle, CbmKFParticleFinderQA *kfparticleQA)
Definition: CbmAnaConversionKF.cxx:201
CbmAnaConversionKF::Invmass_2gamma
Double_t Invmass_2gamma(KFParticle part1, KFParticle part2)
Definition: CbmAnaConversionKF.cxx:770
CbmAnaConversionKF::SetGhostIds
void SetGhostIds(std::vector< int > *ghostids)
Definition: CbmAnaConversionKF.cxx:218
CbmAnaConversionKF::Invmass_2electrons
Double_t Invmass_2electrons(KFParticle part1, KFParticle part2)
Definition: CbmAnaConversionKF.cxx:788
CbmAnaConversionKF::fhPi0_mass
TH1D * fhPi0_mass
Definition: CbmAnaConversionKF.h:98
xMath::Pi
double Pi()
Definition: xMath.h:5
CbmAnaConversionKF::fMcTracks
TClonesArray * fMcTracks
Definition: CbmAnaConversionKF.h:73
CbmAnaConversionKF::fhPi0_NDaughters
TH1D * fhPi0_NDaughters
Definition: CbmAnaConversionKF.h:88
CbmAnaConversionKF::~CbmAnaConversionKF
virtual ~CbmAnaConversionKF()
Definition: CbmAnaConversionKF.cxx:73
CbmAnaConversionKF::fKFMcParticles
TClonesArray * fKFMcParticles
Definition: CbmAnaConversionKF.h:72
CbmKFParticleFinderQA.h
CbmAnaConversionKF::fStsTrackMatches
TClonesArray * fStsTrackMatches
Definition: CbmAnaConversionKF.h:75
CbmAnaConversionKF::ReconstructGammas
void ReconstructGammas()
Definition: CbmAnaConversionKF.cxx:717
CbmAnaConversionKF::timer
TStopwatch timer
Definition: CbmAnaConversionKF.h:131
CbmMCTrack.h
CbmAnaConversionKF::electronIDs
std::vector< int > electronIDs
Definition: CbmAnaConversionKF.h:111
CbmAnaConversionKF::fKFparticleFinderQA
CbmKFParticleFinderQA * fKFparticleFinderQA
Definition: CbmAnaConversionKF.h:78
CbmKFParticleFinder::GetTopoReconstructor
const KFParticleTopoReconstructor * GetTopoReconstructor() const
Definition: CbmKFParticleFinder.h:46
CbmMCTrack
Definition: CbmMCTrack.h:34
CbmAnaConversionKF::fHistoList_kfparticle
std::vector< TH1 * > fHistoList_kfparticle
Definition: CbmAnaConversionKF.h:106
CbmAnaConversionKF::fhKF_particlevector
TH1D * fhKF_particlevector
Definition: CbmAnaConversionKF.h:118
CbmAnaConversionKF::fGhostIds
std::vector< int > fGhostIds
Definition: CbmAnaConversionKF.h:102
CbmAnaConversionKF::CombineElectrons
void CombineElectrons()
Definition: CbmAnaConversionKF.cxx:565
CbmAnaConversionKF::fhKF_invmass_fullReco
TH1D * fhKF_invmass_fullReco
Definition: CbmAnaConversionKF.h:127
CbmTrackMatchNew
Definition: CbmTrackMatchNew.h:19
CbmMCTrack::Get4Momentum
void Get4Momentum(TLorentzVector &momentum) const
Definition: CbmMCTrack.h:177
CbmAnaConversionKF::fKFtopoPerf
KFTopoPerformance * fKFtopoPerf
Definition: CbmAnaConversionKF.h:81
CbmAnaConversionKF::particleMatch
std::vector< KFPartMatch > particleMatch
Definition: CbmAnaConversionKF.h:110
CbmStsTrack
Definition: CbmStsTrack.h:37
M2E
#define M2E
Definition: CbmAnaConversionKF.cxx:26
CbmAnaConversionKF::Reconstruct
void Reconstruct()
Definition: CbmAnaConversionKF.cxx:620
CbmAnaConversionKF::CbmAnaConversionKF
CbmAnaConversionKF()
Definition: CbmAnaConversionKF.cxx:31
CbmAnaConversionKF::fhKF_NofPi0_signal
TH1D * fhKF_NofPi0_signal
Definition: CbmAnaConversionKF.h:122