FairRoot/PandaRoot
PndMyAnalysisTask.cxx
Go to the documentation of this file.
1 // ******************************************************
2 // Psi' -> J/psi (-> e+ e-) pi+ pi- Analysis Example Task
3 //
4 // K.Goetzen 7/2012
5 // *******************************************************
6 
7 #include "TClonesArray.h"
8 
9 #include "FairRootManager.h"
10 #include "FairRunAna.h"
11 #include "FairRuntimeDb.h"
12 
13 #include "TVector3.h"
14 #include "TH1F.h"
15 #include "TH2F.h"
16 #include "TStopwatch.h"
17 
18 #include "FairRun.h"
19 #include "FairRuntimeDb.h"
20 #include <string>
21 #include <iostream>
22 
23 // RHO stuff
24 #include "RhoCandidate.h"
25 #include "RhoHistogram/RhoTuple.h"
26 #include "RhoFactory.h"
27 
28 // Analysis Tools
29 #include "PndMyAnalysisTask.h"
30 #include "PndAnalysis.h"
31 
32 // Fitters
33 #include "Rho4CFitter.h"
34 #include "RhoKinVtxFitter.h"
35 #include "RhoKinFitter.h"
36 
37 
38 using std::cout;
39 using std::endl;
40 
41 
42 // ----- Default constructor -------------------------------------------
44  FairTask("Panda Analysis Task")
45 {
46 }
47 // -------------------------------------------------------------------------
48 
49 // ----- Destructor ----------------------------------------------------
51 // -------------------------------------------------------------------------
52 
53 
54 
55 // ----- Public method Init --------------------------------------------
57 {
58  //cout << " Inside the Init function****" << endl;
59 
60 
61  // initialize analysis object
63 
64  // Create and register output array
65  /* hjpsimass = new TH1F("hjpsimass","jpsi cands",100,0.0,4.5);
66  hpsimass = new TH1F("hpsimass","psi cands",100,0.,4.5);*/
67 
68  hjpsim_nopid = new TH1F("hjpsim_nopid","J/#psi mass (no pid)",200,0,4);
69  hpsim_nopid = new TH1F("hpsim_nopid","#psi(2S) mass (no pid)",200,0,5);
70  hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4);
71  hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5);
72  hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (true pid)",200,0,4);
73  hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (true pid)",200,0,5);
74 
75  hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4);
76  hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass ( truth match)",200,0,5);
77  hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4);
78  hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5);
79 
80  hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2);
81  hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2);
82 
83  hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass vertex fit",200,0,4);
84  hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4);
85  hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (mass constraint fit)",200,0,4);
86 
87  hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf","J/#psi, #chi^{2} vertex fit",100,0,10);
88  hpsi_chi2_4c = new TH1F("hpsi_chi2_4c","#psi, #chi^{2} 4C fit",100,0,250);
89  hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf","J/#psi, #chi^{2} mass fit",100,0,10);
90 
91  hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2);
92 
93  ntp=new RhoTuple("ntp","the J/psi ntuple");
94  ntp2=new RhoTuple("ntp2","the psi(2S) ntuple");
95 
96  // **** mass selector and McTruthMatcher
97  //
98  jpsiMassSel = new RhoMassParticleSelector("jpsiSelector" , 3.097, 0.6);
99 
100  evcount=0;
101  epmax=0;
102  emmax=0;
103  pipmax=0;
104  pimmax=0;
105  mcmax=0;
106 
107  cout << "-I- PndMyAnalysisTask: Intialization successfull" << endl;
108 
109  timer=new TStopwatch();
110  timer->Start();
111 
112  return kSUCCESS;
113 }
114 
116 {
117 
118  // Get run and runtime database
119  FairRun* run = FairRun::Instance();
120  if ( ! run ) { Fatal("SetParContainers", "No analysis run"); }
121 
122  //FairRuntimeDb* db = run->GetRuntimeDb();
123  //if ( ! db ) Fatal("SetParContainers", "No runtime database");
124 
125 
126 }
127 
128 // -------------------------------------------------------------------------
129 
130 // ----- Public method Exec --------------------------------------------
131 void PndMyAnalysisTask::Exec(Option_t*)
132 {
133  Int_t j=0;
134 
136 
137  if (!(++evcount%100)) {
138  cout <<"evt "<<evcount;
139  timer->Stop();
140  cout <<" t="<< timer->CpuTime();
141  timer->Start();
142  cout <<" Cand Watermark = "<<RhoFactory::Instance()->GetCandidateWatermark();
143  cout <<" epmax:"<<epmax;
144  cout <<" emmax:"<<emmax;
145  cout <<" pipmax:"<<pipmax;
146  cout <<" pimmax:"<<pimmax;
147  cout <<" mcmax:"<<mcmax;
148  cout <<endl;
149  }
150 
151 
152  // **** create all the particle lists we'll need for rebuilding the decay tree
153  //
154  //RhoCandList eplus, eminus, piplus, piminus, jpsi, psi, mctrk;
155 
156  RhoCandList all, chrg, el, eplus, eminus, piplus, piminus, jpsi, jpsi2, psi2s;
157  RhoCandList epsel, emsel;
158 
159  TLorentzVector ini(0, 0, 6.231552, 7.240065);
160 
161 
162 // cout <<" #### mct="<<mctrk.GetLength()<<endl;
163 
164  // *** Select with no PID info ('All'); type and mass are set
165  theAnalysis->FillList(eplus, "ElectronAllPlus","PidAlgoEmcBayes");
166  theAnalysis->FillList(eminus, "ElectronAllMinus","PidAlgoEmcBayes");
167  theAnalysis->FillList(piplus, "PionAllPlus","PidAlgoEmcBayes");
168  theAnalysis->FillList(piminus, "PionAllMinus","PidAlgoEmcBayes");
169 
170 
171  int nep = eplus.GetLength();
172  int nem = eminus.GetLength();
173  int npip = piplus.GetLength();
174  int npim = piminus.GetLength();
175 
176  if (nep>epmax) { epmax=nep; }
177  if (nem>emmax) { emmax=nem; }
178  if (npip>pipmax) { pipmax=npip; }
179  if (npim>pimmax) { pimmax=npim; }
180 
181 // jpsi.Combine(eplus, eminus);
182 // FillMassHisto(hjpsimass, jpsi);
183 // jpsi.Select(jpsiMSel);
184 //
185 // psi.Combine(jpsi, piplus, piminus);
186 // FillMassHisto(hpsimass, psi);
187 
188 
189  jpsi.Combine(eplus, eminus);
190 
191  // *** do the truth match for jpsi
192  jpsi.SetType("J/psi");
193  for (j=0; j<jpsi.GetLength(); ++j) {
194  hjpsim_nopid->Fill( jpsi[j]->M() );
195  if (theAnalysis->McTruthMatch(jpsi[j])) { hjpsim_ftm->Fill( jpsi[j]->M() ); }
196  else { hjpsim_nm->Fill( jpsi[j]->M() ); }
197  }
198 
199  // *** do vertex fitting (J/psi)
200  for (j=0; j<jpsi.GetLength(); ++j) {
201  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
202 
203  vtxfitter.Fit();
204  const RhoCandidate* jfit = jpsi[j]->GetFit(); // access the fitted cand
205  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
206 
207  double chi2_vtx=vtxfitter.GetChi2(); // access chi2 of fit
208  hjpsi_chi2_vf->Fill(chi2_vtx);
209 
210  if (chi2_vtx<2) { // when good enough, fill some histos
211  hjpsim_vf->Fill(jfit->M());
212  hvpos->Fill(jVtx.X(),jVtx.Y());
213  }
214 
215  bool match = theAnalysis->McTruthMatch(jpsi[j]);
216 
217  //ntp->Fill(jpsi[j]->M(),jfit->M(),chi2_vtx,jVtx.X(),jVtx.Y(),jVtx.Z(),match);
218 
219  RhoCandidate* ep = jpsi[j]->Daughter(0);
220  RhoCandidate* em = jpsi[j]->Daughter(1);
221 
222  ntp->Column("jpsim", jpsi[j]->M(), -1.);
223  ntp->Column("jpsitht", jpsi[j]->P4().Theta(), -999.);
224  ntp->Column("jpsip", jpsi[j]->P(), -1.);
225  ntp->Column("jpsiphi", jpsi[j]->P4().Phi(), -999.);
226 
227  ntp->Column("eptht", ep->P4().Theta(), -999.);
228  ntp->Column("epp", ep->P(), -1.);
229  ntp->Column("epphi", ep->P4().Phi(), -999.);
230  ntp->Column("eppid0", ep->GetPidInfo(0), -1.);
231  ntp->Column("eppid1", ep->GetPidInfo(1), -1.);
232  ntp->Column("eppid2", ep->GetPidInfo(2), -1.);
233  ntp->Column("eppid3", ep->GetPidInfo(3), -1.);
234  ntp->Column("eppid4", ep->GetPidInfo(4), -1.);
235 
236  ntp->Column("emtht", em->P4().Theta(), -999.);
237  ntp->Column("emp", em->P(), -1.);
238  ntp->Column("emphi", em->P4().Phi(), -999.);
239  ntp->Column("empid0", em->GetPidInfo(0), -1.);
240  ntp->Column("empid1", em->GetPidInfo(1), -1.);
241  ntp->Column("empid2", em->GetPidInfo(2), -1.);
242  ntp->Column("empid3", em->GetPidInfo(3), -1.);
243  ntp->Column("empid4", em->GetPidInfo(4), -1.);
244 
245  ntp->Column("jpsifitm", jfit->M(), -1.);
246  ntp->Column("vchi2", chi2_vtx, 9999.);
247  ntp->Column("vx", jVtx.X(), 9999.);
248  ntp->Column("vy", jVtx.Y(), 9999.);
249  ntp->Column("vz", jVtx.Z(), 9999.);
250  ntp->Column("jmct", match, kFALSE);
251 
252  ntp->DumpData();
253  }
254 
255  // *** some rough mass selection
256  //jpsi.Select(jpsiMassSel);
257 
258  psi2s.Combine(jpsi, piplus, piminus);
259 
260  // *** do the truth match for psi2s
261  psi2s.SetType("psi(2S)0");
262 
263  for (j=0; j<psi2s.GetLength(); ++j) {
264  hpsim_nopid->Fill( psi2s[j]->M() );
265  if (theAnalysis->McTruthMatch(psi2s[j])) {
266  hpsim_ftm->Fill( psi2s[j]->M() );
267 
269  const RhoCandidate* truePsi = psi2s[j]->GetMcTruth();
270  const RhoCandidate* trueJ = psi2s[j]->Daughter(0)->GetMcTruth();
271 
272  hjpsim_diff->Fill(trueJ->M()-psi2s[j]->Daughter(0)->M());
273  hpsim_diff->Fill(truePsi->M()-psi2s[j]->M());
274  } else { hpsim_nm->Fill( psi2s[j]->M() ); }
275  }
276 
277  // *** do 4c fit (initial psi(2S) system)
278  for (j=0; j<psi2s.GetLength(); ++j) {
279  Rho4CFitter fitter(psi2s[j],ini);
280  fitter.FitConserveMasses();
281 
282  double chi2_4c=fitter.GetChi2();
283  hpsi_chi2_4c->Fill(chi2_4c);
284 
285  //const RhoCandidate* jfit = (psi2s[j]->Daughter(0))->GetFit();
286 
287  const RhoCandidate* epfit = (psi2s[j]->Daughter(0)->Daughter(0))->GetFit();
288  const RhoCandidate* emfit = (psi2s[j]->Daughter(0)->Daughter(1))->GetFit();
289 
290  TLorentzVector tlvepf = epfit->P4();
291  TLorentzVector tlvemf = emfit->P4();
292 
293  hjpsim_4cf->Fill((tlvepf+tlvemf).M());
294 
295 
296 
297  RhoCandidate* ep = psi2s[j]->Daughter(0)->Daughter(0);
298  RhoCandidate* em = psi2s[j]->Daughter(0)->Daughter(1);
299  RhoCandidate* pip = psi2s[j]->Daughter(1);
300  RhoCandidate* pim = psi2s[j]->Daughter(2);
301  RhoCandidate* jpsic = psi2s[j]->Daughter(0);
302 
303  bool match = theAnalysis->McTruthMatch(psi2s[j]);
304  bool matchj = theAnalysis->McTruthMatch(jpsic);
305 
306  ntp2->Column("psim", psi2s[j]->M(), -1.);
307  ntp2->Column("fcchi2", chi2_4c, 9999.);
308  ntp2->Column("psimct", match, kFALSE);
309 
310  ntp2->Column("jmct", matchj, kFALSE);
311  ntp2->Column("jfitm", (tlvepf+tlvemf).M(), -1.);
312  ntp2->Column("jpsim", jpsic->M(), -1.);
313  ntp->Column("jpsitht", jpsic->P4().Theta(), -999.);
314  ntp->Column("jpsip", jpsic->P(), -1.);
315  ntp->Column("jpsiphi", jpsic->P4().Phi(), -999.);
316 
317  ntp2->Column("eptht", ep->P4().Theta(), -999.);
318  ntp2->Column("epp", ep->P(), -1.);
319  ntp2->Column("epphi", ep->P4().Phi(), -999.);
320  ntp2->Column("eppid0", ep->GetPidInfo(0), -1.);
321  ntp2->Column("eppid1", ep->GetPidInfo(1), -1.);
322  ntp2->Column("eppid2", ep->GetPidInfo(2), -1.);
323  ntp2->Column("eppid3", ep->GetPidInfo(3), -1.);
324  ntp2->Column("eppid4", ep->GetPidInfo(4), -1.);
325 
326  ntp2->Column("emtht", em->P4().Theta(), -999.);
327  ntp2->Column("emp", em->P(), -1.);
328  ntp2->Column("emphi", em->P4().Phi(), -999.);
329  ntp2->Column("empid0", em->GetPidInfo(0), -1.);
330  ntp2->Column("empid1", em->GetPidInfo(1), -1.);
331  ntp2->Column("empid2", em->GetPidInfo(2), -1.);
332  ntp2->Column("empid3", em->GetPidInfo(3), -1.);
333  ntp2->Column("empid4", em->GetPidInfo(4), -1.);
334 
335  ntp2->Column("piptht", pip->P4().Theta(), -999.);
336  ntp2->Column("pipp", pip->P(), -1.);
337  ntp2->Column("pipphi", pip->P4().Phi(), -999.);
338  ntp2->Column("pippid0", pip->GetPidInfo(0), -1.);
339  ntp2->Column("pippid1", pip->GetPidInfo(1), -1.);
340  ntp2->Column("pippid2", pip->GetPidInfo(2), -1.);
341  ntp2->Column("pippid3", pip->GetPidInfo(3), -1.);
342  ntp2->Column("pippid4", pip->GetPidInfo(4), -1.);
343 
344  ntp2->Column("pimtht", pim->P4().Theta(), -999.);
345  ntp2->Column("pimp", pim->P(), -1.);
346  ntp2->Column("pimphi", pim->P4().Phi(), -999.);
347  ntp2->Column("pimpid0", pim->GetPidInfo(0), -1.);
348  ntp2->Column("pimpid1", pim->GetPidInfo(1), -1.);
349  ntp2->Column("pimpid2", pim->GetPidInfo(2), -1.);
350  ntp2->Column("pimpid3", pim->GetPidInfo(3), -1.);
351  ntp2->Column("pimpid4", pim->GetPidInfo(4), -1.);
352 
353  ntp2->DumpData();
354 
355  }
356 
357  // do mass constraint fit
358  for (j=0; j<jpsi.GetLength(); ++j) {
359  RhoKinFitter mfitter(jpsi[j]);
360  mfitter.AddMassConstraint(3.0965);
361  mfitter.Fit();
362  double chi2_m = mfitter.GetChi2();
363  hjpsi_chi2_mf->Fill(chi2_m);
364 
365  if (chi2_m<2) { hjpsim_mcf->Fill(jpsi[j]->M()); }
366  }
367 
368 
369  // *** all combinatorics again with true PID
370  jpsi.Combine(eplus, eminus);
371  for (j=0; j<jpsi.GetLength(); ++j) { hjpsim_tpid->Fill( jpsi[j]->M() ); }
372 
373  //jpsi.Select(jpsiMassSel);
374 
375  psi2s.Combine(jpsi, piplus, piminus);
376  for (j=0; j<psi2s.GetLength(); ++j) { hpsim_tpid->Fill( psi2s[j]->M() ); }
377 
378  // *** and again with PidAlgoEmcBayes and loose selection
379  theAnalysis->FillList(eplus, "ElectronLoosePlus","PidAlgoEmcBayes");
380  theAnalysis->FillList(eminus, "ElectronLooseMinus","PidAlgoEmcBayes");
381  theAnalysis->FillList(piplus, "PionLoosePlus","PidAlgoEmcBayes");
382  theAnalysis->FillList(piminus, "PionLooseMinus","PidAlgoEmcBayes");
383 
384  jpsi.Combine(eplus, eminus);
385  for (j=0; j<jpsi.GetLength(); ++j) { hjpsim_lpid->Fill( jpsi[j]->M() ); }
386 
387  //jpsi.Select(jpsiMassSel);
388  //SelectMass(jpsi,3.096,1.0);
389 
390  psi2s.Combine(jpsi, piplus, piminus);
391  for (j=0; j<psi2s.GetLength(); ++j) { hpsim_lpid->Fill( psi2s[j]->M() ); }
392 
393 }
394 // -------------------------------------------------------------------------
395 
397 {
398  Int_t i=0;
399  for (i=0; i<l.GetLength(); ++i) {
400  h->Fill(l[i]->M());
401  }
402 }
403 
404 // -------------------------------------------------------------------------
405 
406 int PndMyAnalysisTask::SelectPdgCode(RhoCandList& , RhoCandList& l) // mct //[R.K.03/2017] unused variable(s)
407 {
408  int removed = 0;
409  int pdgcode=0;
410  //int nmct = mct.GetLength();
411 
412  //PndMcTruthMatch mcm;
413 
414  if (l.GetLength()>0) { pdgcode = l[0]->PdgCode(); }
415 
416  for (int ii=0; ii<l.GetLength(); ++ii) {
417  RhoCandidate* mccnd = l[ii]->GetMcTruth();
418  if (mccnd)
419  if(mccnd->PdgCode()==pdgcode)
420  continue; // we skip if things are allright and remove if not
421  l.Remove(l[ii]);
422  removed++;
423  }
424  return removed;
425 }
426 
427 // -------------------------------------------------------------------------
428 
430 {
431  hjpsim_nopid->Write();
432  hpsim_nopid->Write();
433  hjpsim_lpid->Write();
434  hpsim_lpid->Write();
435  hjpsim_tpid->Write();
436  hpsim_tpid->Write();
437  hjpsim_ftm->Write();
438  hpsim_ftm->Write();
439  hjpsim_nm->Write();
440  hpsim_nm->Write();
441 
442  hpsim_diff->Write();
443  hjpsim_diff->Write();
444 
445  hjpsim_vf->Write();
446  hjpsim_4cf->Write();
447  hjpsim_mcf->Write();
448 
449  hjpsi_chi2_vf->Write();
450  hpsi_chi2_4c->Write();
451  hjpsi_chi2_mf->Write();
452 
453  hvpos->Write();
454 
455  ntp->GetInternalTree()->Write();
456  ntp2->GetInternalTree()->Write();
457  //ntp->WriteToFile("outtask_ntp.root");
458  //ntp2->AddToFile("outtask_ntp.root");
459 
460 // hjpsimass->Write();
461 // hpsimass->Write();
462 
463 }
464 
void AddMassConstraint(double mass)
Int_t run
Definition: autocutx.C:47
Int_t i
Definition: run_full.C:25
Bool_t FitConserveMasses()
Definition: Rho4CFitter.cxx:64
TVector3 Pos() const
Definition: RhoCandidate.h:186
virtual void SetParContainers()
Int_t GetLength() const
Definition: RhoCandList.h:46
PndAnalysis * theAnalysis
Bool_t FillList(RhoCandList &l, TString listkey="All", TString pidTcaNames="", int trackHypothesis=-1)
CandList piplus
CandList eplus
void FillMassHisto(TH1F *h, RhoCandList &l)
void Combine(RhoCandList &l1, RhoCandList &l2)
int SelectPdgCode(RhoCandList &mct, RhoCandList &l)
TLorentzVector P4() const
Definition: RhoCandidate.h:195
virtual void Exec(Option_t *opt)
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
Definition: RhoTuple.cxx:56
void SetType(const TParticlePDG *pdt, int start=0)
static RhoFactory * Instance()
Definition: RhoFactory.cxx:34
Double_t M() const
Int_t Remove(RhoCandidate *)
Double_t P() const
void DumpData()
Definition: RhoTuple.cxx:391
static Int_t GetCandidateWatermark()
Definition: RhoFactory.cxx:110
RhoMassParticleSelector * jpsiMassSel
GeV c P
CandList piminus
ClassImp(PndAnaContFact)
double GetChi2() const
Definition: RhoFitterBase.h:48
TTree * GetInternalTree()
Definition: RhoTuple.h:207
Bool_t McTruthMatch(RhoCandidate *cand, Int_t level=2, bool verbose=false)
Int_t GetEvent(Int_t n=-1)
Bool_t Fit()
virtual InitStatus Init()
double GetPidInfo(int hypo)
CandList eminus