FairRoot/PandaRoot
PndSimpleAnalysis.cxx
Go to the documentation of this file.
1 /******************************************************
2 Class PndSimpleAnalysis
3 
4 Task to perform analysis (combinatorics, simple filetring)
5 by setting a configuration file
6 
7 Writes out a TTree 'ntp'
8 
9 Author: K.Goetzen, GSI, 06/2008
10 
11 *******************************************************/
12 
13 //Base headers
14 #include <string>
15 #include <iostream>
16 #include <fstream>
17 #include <algorithm>
18 
19 //FAIR headers
20 #include "FairRootManager.h"
21 #include "FairRunAna.h"
22 #include "FairRuntimeDb.h"
23 #include "FairRun.h"
24 
25 //ROOT headers
26 #include "TVector3.h"
27 #include "TLorentzVector.h"
28 #include "TH1F.h"
29 #include "TClonesArray.h"
30 #include "TDatabasePDG.h"
31 #include "TParticlePDG.h"
32 #include "TTree.h"
33 
34 //Fastsim headers
35 #include "PndSimpleAnalysis.h"
36 #include "StrTok.h"
37 #include "ArgList.h"
38 #include "PndListDefiner.h"
39 
40 //RHO headers
41 #include "RhoBase/RhoCandidate.h"
42 #include "RhoBase/RhoCandList.h"
52 #include "RhoBase/RhoFactory.h"
53 #include "PndPidCandidate.h"
54 #include "PndPidProbability.h"
55 
56 
57 using std::cout;
58 using std::endl;
59 using std::ios;
60 
61 
62 // ----- Default constructor -------------------------------------------
64  FairTask("Panda Analysis Task")
65 {
66  }
67 
68 // -------------------------------------------------------------------------
69 
71  FairTask("Panda Analysis Task")
72 {
73  SetConfigFile(filename);
74 }
75 
76 // -------------------------------------------------------------------------
77 
78 // ----- Destructor ----------------------------------------------------
80 // -------------------------------------------------------------------------
81 
82 
83 
84 // ----- Public method Init --------------------------------------------
86 {
87 
88  //cout << " Inside the Init function****" << endl;
89 
90  //FairDetector::Initialize();
91  //FairRun* sim = FairRun::Instance();
92  //FairRuntimeDb* rtdb=sim->GetRuntimeDb();
93 
94  // Get RootManager
95  FairRootManager* ioman = FairRootManager::Instance();
96  if ( ! ioman ) {
97  cout << "-E- PndSimpleAnalysis::Init: "
98  << "RootManager not instantiated!" << endl;
99  return kFATAL;
100  }
101 
102  //Prepare all generic Listnames
104  InitColumnNames();
105 
106  //ntp=new TTree("ntp","PndSimpleAnalysis NTuple");
107  ntp=0;
108 
109  // set the initial 4 vector of the pbar p system
110  // can be used for missing mass and mom in CMS
111  // has to be set in config file for the time being
112  fpInit.SetXYZT(0.,0.,0.,0.);
113 
114 // Get input array
115  fChargedArray = (TClonesArray*) ioman->GetObject("PidChargedCand");
116  fNeutralArray = (TClonesArray*) ioman->GetObject("PidNeutralCand");
117 
118  fChargedProbability = (TClonesArray*) ioman->GetObject("PidChargedProbability");
119  fNeutralProbability = (TClonesArray*) ioman->GetObject("PidNeutralProbability");
120 
121  fMcArray = (TClonesArray*) ioman->GetObject("PndMcTracks");
122  fMicroArray = 0;//(TClonesArray*) ioman->GetObject("PndPidCandidates");
123 
124  if ( !fChargedArray && !fNeutralArray && !fMcArray && !fMicroArray) {
125  cout << "-W- PndSimpleAnalysis::Init: "
126  << "None of PidChargedCand, PndNeutralCand, PidMcTracks, PndPidCandidates available!" << endl;
127  return kERROR;
128  }
129 
130  if (!SetupAnalysis()) {
131  cout << "-E- PndSimpleAnalysis::Init: "
132  << "Error reading config file "<<fCfgFileName<<"." << endl;
133  return kFATAL;
134  }
135 
136  // Create and register output array
137  cout << "-I- PndSimpleAnalysis: Intialization successfull" << endl;
138 
139  // **** create and configure the selectors/filters we'd like to use later
140  //
141 
144 
145 
146  // **** pid selectors for generic list creation
147  //
151  kSel = new RhoSimpleKaonSelector();
153 
154 
155  evcount=0;
156 
157  //fMaxEntries = 200;
158 
159  return kSUCCESS;
160 
161 }
162 // -------------------------------------------------------------------------
163 
165 {
166 
167  //very Loose Lists
168  fGenericListNames.push_back("ElectronVeryLooseP");
169  fGenericListNames.push_back("ElectronVeryLooseM");
170  fGenericListNames.push_back("MuonVeryLooseP");
171  fGenericListNames.push_back("MuonVeryLooseM");
172  fGenericListNames.push_back("PionVeryLooseP");
173  fGenericListNames.push_back("PionVeryLooseM");
174  fGenericListNames.push_back("KaonVeryLooseP");
175  fGenericListNames.push_back("KaonVeryLooseM");
176  fGenericListNames.push_back("ProtonVeryLooseP");
177  fGenericListNames.push_back("ProtonVeryLooseM");
178 
179  //Loose Lists
180  fGenericListNames.push_back("ElectronLooseP");
181  fGenericListNames.push_back("ElectronLooseM");
182  fGenericListNames.push_back("MuonLooseP");
183  fGenericListNames.push_back("MuonLooseM");
184  fGenericListNames.push_back("PionLooseP");
185  fGenericListNames.push_back("PionLooseM");
186  fGenericListNames.push_back("KaonLooseP");
187  fGenericListNames.push_back("KaonLooseM");
188  fGenericListNames.push_back("ProtonLooseP");
189  fGenericListNames.push_back("ProtonLooseM");
190 
191  //tight Lists
192  fGenericListNames.push_back("ElectronTightP");
193  fGenericListNames.push_back("ElectronTightM");
194  fGenericListNames.push_back("MuonTightP");
195  fGenericListNames.push_back("MuonTightM");
196  fGenericListNames.push_back("PionTightP");
197  fGenericListNames.push_back("PionTightM");
198  fGenericListNames.push_back("KaonTightP");
199  fGenericListNames.push_back("KaonTightM");
200  fGenericListNames.push_back("ProtonTightP");
201  fGenericListNames.push_back("ProtonTightM");
202 
203  //very tight Lists
204  fGenericListNames.push_back("ElectronVeryTightP");
205  fGenericListNames.push_back("ElectronVeryTightM");
206  fGenericListNames.push_back("MuonVeryTightP");
207  fGenericListNames.push_back("MuonVeryTightM");
208  fGenericListNames.push_back("PionVeryTightP");
209  fGenericListNames.push_back("PionVeryTightM");
210  fGenericListNames.push_back("KaonVeryTightP");
211  fGenericListNames.push_back("KaonVeryTightM");
212  fGenericListNames.push_back("ProtonVeryTightP");
213  fGenericListNames.push_back("ProtonVeryTightM");
214 
215 
216  // the neutral List
217  fGenericListNames.push_back("Neutral");
218  fGenericListNames.push_back("McTruth");
219 
220  unsigned int i=0;
221 
222  for (i=0; i<fGenericListNames.size(); i++) {
224  }
225 
226  int pcodes[10] = {-11,11,-13,13,211,-211,321,-321,2212,-2212};
227 
228  for (i=0; i<40; i++) {
229  PndListDefiner* ldef=new PndListDefiner;
230  ldef->fName = fGenericListNames[i];
231  ldef->fPdgCode = pcodes[i%10];
232  ldef->fCharge = (float)(1-(i%2*2));
233  ldef->fIsAntiList = (i%2)?true:false;
234  ldef->fAntiIdx = (i%2) ? i-1 : i+1;
235  ldef->fIsGeneric = true;
236  fListDefiners.push_back(ldef);
237  }
238 
239  // the neutral list
240  PndListDefiner* ldef=new PndListDefiner;
241  ldef->fName = fGenericListNames[40];
242  ldef->fPdgCode = 22;
243  ldef->fCharge = 0.0;
244  ldef->fIsAntiList = false;
245  ldef->fAntiIdx = 40;
246  ldef->fIsGeneric = true;
247  fListDefiners.push_back(ldef);
248 
249  // mc truth list
250  ldef=new PndListDefiner;
251  ldef->fName = fGenericListNames[41];
252  ldef->fIsAntiList = false;
253  ldef->fAntiIdx = 41;
254  ldef->fIsGeneric = false;
255  fListDefiners.push_back(ldef);
256 
257  //for (i=0;i<fListDefiners.size();i++) {cout << i <<" : "; fListDefiners[i]->Print();}
258 }
259 // -------------------------------------------------------------------------
260 
262 {
263  // mapping of column name and index
264 
265  fColKeyMap.clear();
266  fColShortKeyMap.clear();
267 
268  // indices <500 are floats
269  fColKeyMap["px"] = 100; // px
270  fColKeyMap["py"] = 101; // py
271  fColKeyMap["pz"] = 102; // pz
272  fColKeyMap["en"] = 103; // energy of 4 vector
273  fColKeyMap["m"] = 104; // mass
274  fColKeyMap["ch"] = 105; // charge
275  fColKeyMap["miss"] = 113; // missing mass
276 
277  fColKeyMap["vx"] = 106; // vertex pos x
278  fColKeyMap["vy"] = 107; // vertex pos y
279  fColKeyMap["vz"] = 108; // vertex pos z
280 
281  fColKeyMap["phi"] = 109; // phi of 3 vector
282  fColKeyMap["tht"] = 110; // theta
283  fColKeyMap["cth"] = 111; // cos(theta)
284  fColKeyMap["p"] = 112; // absolute value of momentum
285 
286  fColKeyMap["pxcm"] = 120; // px in cms
287  fColKeyMap["pycm"] = 121; // py in cms
288  fColKeyMap["pzcm"] = 122; // pz in cms
289  fColKeyMap["encm"] = 123; // energy in cms
290 
291  fColKeyMap["phicm"] = 124; // phi of 3 vector in cms
292  fColKeyMap["thtcm"] = 125; // theta in cms
293  fColKeyMap["cthcm"] = 126; // cos(theta) in cms
294  fColKeyMap["pcm"] = 127; // absolute value of momentum in cms
295 
296  fColKeyMap["dm"] = 149;
297  fColKeyMap["d1m"] = 150; // directly store daughter masses
298  fColKeyMap["d2m"] = 151; // directly store daughter masses
299  fColKeyMap["d3m"] = 152; // directly store daughter masses
300  fColKeyMap["d4m"] = 153; // directly store daughter masses
301  fColKeyMap["d5m"] = 154; // directly store daughter masses
302 
303  fColKeyMap["tcb"] = 160; // PID: Barrel DIRC theta_C
304  fColKeyMap["tcd"] = 161; // PID: Disc DIRC theta_C
305  fColKeyMap["tcr"] = 162; // PID: RICH theta_C
306  fColKeyMap["tm2"] = 163; // PID: TOF m^2
307  fColKeyMap["dem"] = 164; // PID: MVD dE/dx
308  fColKeyMap["des"] = 165; // PID: STT dE/dx
309  fColKeyMap["det"] = 166; // PID: TPC dE/dx
310 
311  fColKeyMap["elh"] = 170; // PID: Electron Likelihood
312  fColKeyMap["mulh"] = 171; // PID: Muon Likelihood
313  fColKeyMap["pilh"] = 172; // PID: Pion Likelihood
314  fColKeyMap["klh"] = 173; // PID: Kaon Likelihood
315  fColKeyMap["plh"] = 174; // PID: Proton Likelihood
316 
317  // indeces >500 are ints
318  fColKeyMap["d1"] = 501; // index in daughter list 1
319  fColKeyMap["d2"] = 502; // index in daughter list 2
320  fColKeyMap["d3"] = 503; // index in daughter list 3
321  fColKeyMap["d4"] = 504; // index in daughter list 4
322  fColKeyMap["d5"] = 505; // index in daughter list 5
323 
324  fColKeyMap["mci"] = 510; // index in MC truth list
325  fColKeyMap["sel"] = 515; // PID:selector level accept (1=veryLoose, 2=loose, 3= tight, 4=veryTight)
326 
327  fColKeyMap["pdg"] = 520; // PDG code
328  fColKeyMap["mothi"] = 521; // index of mother particle MCTruth
329 
330 
331  // indeces >1000 map to lists of col names
332 
333  std::vector<std::string> shortcut;
334 
335  fColKeyMap["p4"] = 1000;
336  shortcut.clear();
337  shortcut.push_back("px");
338  shortcut.push_back("py");
339  shortcut.push_back("pz");
340  shortcut.push_back("en");
341  fColShortKeyMap[1000]=shortcut;
342 
343  fColKeyMap["p4cm"] = 1001;
344  shortcut.clear();
345  shortcut.push_back("pxcm");
346  shortcut.push_back("pycm");
347  shortcut.push_back("pzcm");
348  shortcut.push_back("encm");
349  fColShortKeyMap[1001]=shortcut;
350 
351  fColKeyMap["p3"] = 1002;
352  shortcut.clear();
353  shortcut.push_back("phi");
354  shortcut.push_back("tht");
355  shortcut.push_back("cth");
356  shortcut.push_back("p");
357  fColShortKeyMap[1002]=shortcut;
358 
359  fColKeyMap["p3cm"] = 1003;
360  shortcut.clear();
361  shortcut.push_back("phicm");
362  shortcut.push_back("thtcm");
363  shortcut.push_back("cthcm");
364  shortcut.push_back("pcm");
365  fColShortKeyMap[1003]=shortcut;
366 
367  fColKeyMap["pos"] = 1004;
368  shortcut.clear();
369  shortcut.push_back("vx");
370  shortcut.push_back("vy");
371  shortcut.push_back("vz");
372  fColShortKeyMap[1004]=shortcut;
373 
374  fColKeyMap["pid"] = 1005;
375  shortcut.clear();
376  shortcut.push_back("sel");
377  shortcut.push_back("des");
378  shortcut.push_back("dem");
379  shortcut.push_back("det");
380  shortcut.push_back("tcb");
381  shortcut.push_back("tcd");
382  shortcut.push_back("tcr");
383  shortcut.push_back("tm2");
384  fColShortKeyMap[1005]=shortcut;
385 
386  fColKeyMap["pidlh"] = 1006;
387  shortcut.clear();
388  shortcut.push_back("elh");
389  shortcut.push_back("mulh");
390  shortcut.push_back("pilh");
391  shortcut.push_back("klh");
392  shortcut.push_back("plh");
393  fColShortKeyMap[1006]=shortcut;
394 }
395 
396 // -------------------------------------------------------------------------
397 
399 {
400 
401  // Get run and runtime database
402  FairRun* run = FairRun::Instance();
403  if ( ! run ) { Fatal("SetParContainers", "No analysis run"); }
404 
405  //FairRuntimeDb* db = run->GetRuntimeDb();
406  //if ( ! db ) Fatal("SetParContainers", "No runtime database");
407 
408 
409 }
410 
411 // -------------------------------------------------------------------------
412 
413 // ----- Public method Exec --------------------------------------------
414 void PndSimpleAnalysis::Exec(Option_t*)
415 {
416  // set the boost vector
417  TVector3 pInitBoost(0,0,0);
418  if (fpInit.Z()!=0.0) { pInitBoost=-fpInit.BoostVector(); }
419 
420  unsigned int i=0;
421  int j=0,k=0,k2=0;
422 
423  if (!(++evcount%100)) {
424  cout <<"evt "<<evcount<<endl;
425  }
426 
428 
429  bool dumpAList=false;
430 
431  // cache for boosted 4-vectors
432 
433  std::vector<TLorentzVector> lvcache;
434 
435  for (i=0; i<fListDefiners.size(); i++) {
437 
438  if (!cur->fIsUsed) { continue; }
439 
440  int nd=cur->GetNDau();
441 
442  PndListDefiner* ld1,*ld2,*ld3,*ld4,*ld5;
443 
444  // *********************
445  // do combinatorics and fill list
446  // *********************
447  switch (nd) {
448  case 2:
449  ld1=fListDefiners[cur->fDauIdx[0]];
450  ld2=fListDefiners[cur->fDauIdx[1]];
451 
452  cur->fList.Combine(ld1->fList,ld2->fList);
453 
454  break; // case 2
455 
456  case 3:
457  ld1=fListDefiners[cur->fDauIdx[0]];
458  ld2=fListDefiners[cur->fDauIdx[1]];
459  ld3=fListDefiners[cur->fDauIdx[2]];
460 
461  cur->fList.Combine(ld1->fList, ld2->fList, ld3->fList);
462  break; // case 3
463 
464  case 4:
465  ld1=fListDefiners[cur->fDauIdx[0]];
466  ld2=fListDefiners[cur->fDauIdx[1]];
467  ld3=fListDefiners[cur->fDauIdx[2]];
468  ld4=fListDefiners[cur->fDauIdx[3]];
469 
470  cur->fList.Combine(ld1->fList, ld2->fList, ld3->fList, ld4->fList);
471  break; //case 4
472 
473  case 5:
474  ld1=fListDefiners[cur->fDauIdx[0]];
475  ld2=fListDefiners[cur->fDauIdx[1]];
476  ld3=fListDefiners[cur->fDauIdx[2]];
477  ld4=fListDefiners[cur->fDauIdx[3]];
478  ld5=fListDefiners[cur->fDauIdx[4]];
479 
480  cur->fList.Combine(ld1->fList, ld2->fList, ld3->fList, ld4->fList, ld5->fList);
481  break; //case 5
482 
483  } //switch
484 
485  // *********************
486  // select candidates
487  // *********************
488  for (j=0; j<cur->GetNSels(); j++) {
489  if (0!=cur->fSelector[j]) {
490  cur->fList.Select(cur->fSelector[j]);
491  }
492  }
493 
494  // *********************
495  // fill histograms
496  // *********************
497  if (cur->fHisto.size()>0) {
498  for (j=0; j<cur->GetLength(); j++) { cur->fHisto[0]->Fill(cur->fList[j]->M()); }
499  }
500 
501  // *********************
502  // fill ntuple
503  // *********************
504  if (cur->fDumpList) {
505  dumpAList|=cur->fDumpList;
506 
507  PndListDefiner* curdump=cur;
508 
509  // clear the cache of boosted lorentz vectors
510  lvcache.clear();
511 
512  // if AntiList, then append information to arrays of List
513  int off=0;
514  if (cur->fIsAntiList) {
515  // this is the pointer to the corresponding LIST!
516  curdump=fListDefiners[cur->fAntiIdx];
517  off=curdump->GetLength();
518  }
519 
520  // set the number of entries; if list has antilist, take sum of both numbers
521  curdump->fNEntries=off+cur->GetLength();
522 
523  // Dump out the float columns
524  for (j=0; j<(int)curdump->fNtpFNames.size(); j++) {
525  //hmm, unknown name; should not be!
526  if (fColKeyMap.find(curdump->fNtpFNames[j])==fColKeyMap.end()) { continue; }
527 
528  int key=fColKeyMap[curdump->fNtpFNames[j]];
529 
530  bool cachefilled=(lvcache.size()>0);
531 
532  for (k=0; k<cur->GetLength(); k++) {
533  if ( (off+k)>=fMaxEntries ) { continue; }
534  float* theArF=curdump->fNtpFArrays[j];
535 
536  PndPidCandidate* mic=(PndPidCandidate*)((cur->fList[k]->GetRecoCandidate()));
537 
538  // set default
539  theArF[off+k]=0.0;
540 
541  switch (key) {
542  case 100:
543  theArF[off+k]=cur->fList[k]->Px();
544  break;
545  case 101:
546  theArF[off+k]=cur->fList[k]->Py();
547  break;
548  case 102:
549  theArF[off+k]=cur->fList[k]->Pz();
550  break;
551  case 103:
552  theArF[off+k]=cur->fList[k]->E();
553  break;
554 
555  case 104:
556  theArF[off+k]=cur->fList[k]->M();
557  break;
558  case 105:
559  theArF[off+k]=cur->fList[k]->Charge();
560  break;
561 
562  case 106:
563  theArF[off+k]=cur->fList[k]->Pos().X();
564  break;
565  case 107:
566  theArF[off+k]=cur->fList[k]->Pos().Y();
567  break;
568  case 108:
569  theArF[off+k]=cur->fList[k]->Pos().Z();
570  break;
571 
572  case 109:
573  theArF[off+k]=cur->fList[k]->GetVect().Phi();
574  break;
575  case 110:
576  theArF[off+k]=cur->fList[k]->GetVect().Theta();
577  break;
578  case 111:
579  theArF[off+k]=cur->fList[k]->GetVect().CosTheta();
580  break;
581  case 112:
582  theArF[off+k]=cur->fList[k]->GetVect().Mag();
583  break;
584  case 113:
585  theArF[off+k]=( fpInit - (cur->fList[k]->P4()) ).M();
586  break;
587 
588  case 120: // Px in CMS
589  if (!cachefilled) {
590  TLorentzVector lv=cur->fList[k]->P4();
591  lv.Boost(pInitBoost);
592  lvcache.push_back(lv);
593  }
594  theArF[off+k]=lvcache[k].Px();
595  break;
596 
597  case 121: // Py in CMS
598  if (!cachefilled) {
599  TLorentzVector lv=cur->fList[k]->P4();
600  lv.Boost(pInitBoost);
601  lvcache.push_back(lv);
602  }
603  theArF[off+k]=lvcache[k].Py();
604  break;
605 
606  case 122: // Pz in CMS
607  if (!cachefilled) {
608  TLorentzVector lv=cur->fList[k]->P4();
609  lv.Boost(pInitBoost);
610  lvcache.push_back(lv);
611  }
612  theArF[off+k]=lvcache[k].Pz();
613  break;
614 
615  case 123: // E in CMS
616  if (!cachefilled) {
617  TLorentzVector lv=cur->fList[k]->P4();
618  lv.Boost(pInitBoost);
619  lvcache.push_back(lv);
620  }
621  theArF[off+k]=lvcache[k].E();
622  break;
623 
624  case 124: // phi in CMS
625  if (!cachefilled) {
626  TLorentzVector lv=cur->fList[k]->P4();
627  lv.Boost(pInitBoost);
628  lvcache.push_back(lv);
629  }
630  theArF[off+k]=lvcache[k].Vect().Phi();
631  break;
632 
633  case 125: // theta in CMS
634  if (!cachefilled) {
635  TLorentzVector lv=cur->fList[k]->P4();
636  lv.Boost(pInitBoost);
637  lvcache.push_back(lv);
638  }
639  theArF[off+k]=lvcache[k].Vect().Theta();
640  break;
641 
642  case 126: // cos theta in CMS
643  if (!cachefilled) {
644  TLorentzVector lv=cur->fList[k]->P4();
645  lv.Boost(pInitBoost);
646  lvcache.push_back(lv);
647  }
648  theArF[off+k]=lvcache[k].Vect().CosTheta();
649  break;
650 
651  case 127: // |p| in CMS
652  if (!cachefilled) {
653  TLorentzVector lv=cur->fList[k]->P4();
654  lv.Boost(pInitBoost);
655  lvcache.push_back(lv);
656  }
657  theArF[off+k]=lvcache[k].Vect().Mag();
658  break;
659 
660  case 150:
661  theArF[off+k]=cur->fList[k]->Daughter(0)->M();
662  break;
663  case 151:
664  theArF[off+k]=cur->fList[k]->Daughter(1)->M();
665  break;
666  case 152:
667  theArF[off+k]=cur->fList[k]->Daughter(2)->M();
668  break;
669  case 153:
670  theArF[off+k]=cur->fList[k]->Daughter(3)->M();
671  break;
672  case 154:
673  theArF[off+k]=cur->fList[k]->Daughter(4)->M();
674  break;
675 
676  // PID Info based on PndPidCandidate entries
677  case 160:
678  if (mic) { theArF[off+k]=mic->GetDrcThetaC(); }
679  break; // Barrel DIRC tht_c
680  case 161:
681  if (mic) { theArF[off+k]=mic->GetDiscThetaC(); }
682  break; // Disc DIRC tht_c
683  case 162:
684  if (mic) { theArF[off+k]=mic->GetRichThetaC(); }
685  break; // RICH tht_c
686  case 163:
687  if (mic) { theArF[off+k]=mic->GetTofM2(); }
688  break; // TOF m^2
689  case 164:
690  if (mic) { theArF[off+k]=mic->GetMvdDEDX(); }
691  break; // MVD dEdx
692  case 165:
693  if (mic) { theArF[off+k]=mic->GetSttMeanDEDX(); }
694  break; // STT dEdx
695  // fetch likelihoods from proper place, e.g. the rhocand. now with meaning
696  case 170:
697  if (mic) { theArF[off+k]=cur->fList[k]->GetPidInfo(0); }
698  break; // Electron LH
699  case 171:
700  if (mic) { theArF[off+k]=cur->fList[k]->GetPidInfo(1); }
701  break; // Muon LH
702  case 172:
703  if (mic) { theArF[off+k]=cur->fList[k]->GetPidInfo(2); }
704  break; // Pion LH
705  case 173:
706  if (mic) { theArF[off+k]=cur->fList[k]->GetPidInfo(3); }
707  break; // Kaon LH
708  case 174:
709  if (mic) { theArF[off+k]=cur->fList[k]->GetPidInfo(4); }
710  break; // Proton LH
711 
712  // PID Info based in RhoCandidate entries... obsolete
713  /*
714  case 160: theArF[off+k]=cur->fList[k]->GetPidInfo(5); break; // Barrel DIRC tht_c
715  case 161: theArF[off+k]=cur->fList[k]->GetPidInfo(6); break; // Disc DIRC tht_c
716  case 162: theArF[off+k]=cur->fList[k]->GetPidInfo(7); break; // RICH tht_c
717  case 163: theArF[off+k]=cur->fList[k]->GetPidInfo(8); break; // TOF mÏ€
718  case 164: theArF[off+k]=cur->fList[k]->GetPidInfo(9); break; // MVD dEdx
719  case 165: theArF[off+k]=cur->fList[k]->GetPidInfo(10); break; // STT dEdx
720  case 166: theArF[off+k]=cur->fList[k]->GetPidInfo(11); break; // TPC dEdx
721 
722  case 170: theArF[off+k]=cur->fList[k]->GetPidInfo(0); break; // Electron LH
723  case 171: theArF[off+k]=cur->fList[k]->GetPidInfo(1); break; // Muon LH
724  case 172: theArF[off+k]=cur->fList[k]->GetPidInfo(2); break; // Pion LH
725  case 173: theArF[off+k]=cur->fList[k]->GetPidInfo(3); break; // Kaon LH
726  case 174: theArF[off+k]=cur->fList[k]->GetPidInfo(4); break; // Proton LH
727  */
728  }
729 
730  }
731  }
732 
733 
734  // Dump out the int columns
735  for (j=0; j<(int)curdump->fNtpINames.size(); j++) {
736  //hmm, unknown name; should not be!
737  if (fColKeyMap.find(curdump->fNtpINames[j])==fColKeyMap.end()) { continue; }
738 
739  int key=fColKeyMap[curdump->fNtpINames[j]];
740 
741  for (k=0; k<cur->GetLength(); k++) {
742  if ( (off+k)>=fMaxEntries ) { continue; }
743  int* theArI=curdump->fNtpIArrays[j];
744 
745  //set default
746  theArI[k+off]=0;
747 
748  switch (key) {
749  // make the daughter index match
750  case 501:
751  case 502:
752  case 503:
753  case 504:
754  case 505: {
755  int dauindex=-1;
756  int dauoff=0;
757  int dnum=key-501; //number of the daughter to find index of
758 
759  // the dnum-th daughter of our particle
760  RhoCandidate* dauC=cur->fList[k]->Daughter(dnum);
761  if (0==dauC) { break; }
762 
763  int lidx=cur->fDauIdx[dnum];
764  PndListDefiner* dld=fListDefiners[lidx];
765  if (dld->fIsAntiList) {
766  dauoff=fListDefiners[dld->fAntiIdx]->fList.GetLength();
767  }
768 
769  for (k2=0; k2<dld->fList.GetLength(); k2++)
770  if (dauC->IsCloneOf(*(dld->fList[k2]))) {
771  dauindex=k2+dauoff;
772  k2=1000;
773  }
774  theArI[k+off]=dauindex;
775  }
776  break;
777 
778  case 510:
779  Error("PndSimpleAnalysis::Exec","case 510 requested and MC truth handling changed! Please fix it.");
780  //theArI[k+off]=cur->fList[k]->GetMcIdx();
781  break;
782 
783  case 515: // determine the best pid level
784  theArI[k+off]=0;
785  if (cur->fIsGeneric) {
786  int pidl=0;
787  RhoParticleSelectorBase* theSel=0;
788 
789  switch (abs(cur->fPdgCode)) {
790  case 11:
791  theSel=eSel;
792  break;
793  case 13:
794  theSel=muSel;
795  break;
796  case 211:
797  theSel=piSel;
798  break;
799  case 321:
800  theSel=kSel;
801  break;
802  case 2212:
803  theSel=pSel;
804  break;
805  }
806  if (theSel) {
807  pidl=1; //list at least veryLoose
808  if (i>=30) { pidl=4; } //list is already veryTight
809  else if (i>=20) { pidl=3; } //list is already tight
810  else if (i>=10) { pidl=2; } //list is already loose
811 
812  if (pidl<4) {
813  theSel->SetCriterion(veryTight);
814  if (theSel->Accept(cur->fList[k])) { pidl=4; }
815  }
816  if (pidl<3) {
817  theSel->SetCriterion(tight);
818  if (theSel->Accept(cur->fList[k])) { pidl=3; }
819  }
820  if (pidl<2) {
821  theSel->SetCriterion(loose);
822  if (theSel->Accept(cur->fList[k])) { pidl=2; }
823  }
824  }
825  theArI[k+off]=pidl;
826  }
827  break;
828 
829  case 520:
830  theArI[k+off]=cur->fList[k]->PdgCode();
831  break;
832 
833  //case 521: theArI[k+off]=cur->fList[k]->MotherIdx(); break;
834 
835  }
836 
837  }
838  }
839  }
840 
841  }
842 
843  if (dumpAList) { ntp->Fill(); }
844 
845 // for (i=0;i<fListDefiners.size();i++)
846 // if (fListDefiners[i]->fIsUsed) {cout <<i<<" : "; fListDefiners[i]->Print(); }
847 
848 }
849 
850 // -------------------------------------------------------------------------
851 
853 {
854  int i=0;
855  int nd=tc->NDaughters();
856  if (nd==0) { return; }
857  cout <<tc->Uid()<<"("<<tc->GetPidInfo(29)<<") -> ";
858  for (i=0; i<nd; i++) { cout<<tc->Daughter(i)->Uid()<<"("<<tc->Daughter(i)->GetPidInfo(29)<<") "; }
859  cout <<endl;
860  for (i=0; i<nd; i++) { PrintTree(tc->Daughter(i),level+1); }
861  if (level==0) { cout <<endl; }
862 }
863 
864 // -------------------------------------------------------------------------
865 
867 {
869 }
870 
871 // -------------------------------------------------------------------------
872 
874 {
875  if (n=="Charged" || n=="Neutral") { return true; }
876 
877  if (n=="ElectronVeryLoose" || n=="MuonVeryLoose" || n=="PionVeryLoose"
878  || n=="KaonVeryLoose" || n=="ProtonVeryLoose") { return true; }
879 
880  if (n=="ElectronLoose" || n=="MuonLoose" || n=="PionLoose"
881  || n=="KaonLoose" || n=="ProtonLoose") { return true; }
882 
883  if (n=="ElectronTight" || n=="MuonTight" || n=="PionTight"
884  || n=="KaonTight" || n=="ProtonTight") { return true; }
885 
886  if (n=="ElectronVeryTight" || n=="MuonVeryTight" || n=="PionVeryTight"
887  || n=="KaonVeryTight" || n=="ProtonVeryTight") { return true; }
888  else { return false; }
889 }
890 
891 // -------------------------------------------------------------------------
892 
894 {
896 
897  int i, uid=1;
898 
899  // **** loop over all Candidates and add them to the list allCands
900  //
903  mcCands.Cleanup();
904 
905  // when we have a PndPidCandidates Array, take that one
906  // ********** DEPRECATED ***************
907  /*
908  if (fMicroArray)
909  {
910  for (i=0; i<fMicroArray->GetEntriesFast(); i++)
911  {
912  PndPidCandidate *mic = (PndPidCandidate *)fMicroArray->At(i);
913  RhoCandidate tc(*mic,i);
914 
915  if (fabs(tc.Charge())>0.01)
916  chargedCands.Add(tc);
917  else
918  neutralCands.Add(tc);
919  }
920 
921  }
922  else //otherwise look for the older RhoCandidate arrays
923  */
924  // read the charged candidates
925  if (fChargedArray)
926  for (i=0; i<fChargedArray->GetEntriesFast(); i++) {
928  RhoCandidate tcc(*mic,uid++);
929 
930  // are pid data available?
932  if (i<fChargedProbability->GetEntriesFast()) {
934  // numbering see PndPidListMaker
935  tcc.SetPidInfo(0,chProb->GetElectronPidProb());
936  tcc.SetPidInfo(1,chProb->GetMuonPidProb());
937  tcc.SetPidInfo(2,chProb->GetPionPidProb());
938  tcc.SetPidInfo(3,chProb->GetKaonPidProb());
939  tcc.SetPidInfo(4,chProb->GetProtonPidProb());
940  }
941  chargedCands.Add(&tcc);
942  }
943 
944  // read the neutral candidates
945  if (fNeutralArray)
946  for (i=0; i<fNeutralArray->GetEntriesFast(); i++) {
948  RhoCandidate tcn(*mic,uid++);
949 
950  // are pid data available?
952  if (i<fNeutralProbability->GetEntriesFast()) {
954  // numbering see PndPidListMaker
955  tcn.SetPidInfo(0,neuProb->GetElectronPidProb());
956  tcn.SetPidInfo(1,neuProb->GetMuonPidProb());
957  tcn.SetPidInfo(2,neuProb->GetPionPidProb());
958  tcn.SetPidInfo(3,neuProb->GetKaonPidProb());
959  tcn.SetPidInfo(4,neuProb->GetProtonPidProb());
960  }
961  neutralCands.Add(&tcn);
962  }
963 
964  // read the mc truth list
965  if (fMcArray)
966  for (i=0; i<fMcArray->GetEntriesFast(); i++) {
967  RhoCandidate* tc = (RhoCandidate*)fMcArray->At(i);
968  mcCands.Add(tc);
969  }
970 
971  for (i=0; i<40; i++) {
972  if (!fListDefiners[i]->fIsUsed) { continue; }
973  int crithint=i/10;
974  int pidhint=(i%10)/2;
975  //VAbsPidSelector *chSel;
976  RhoParticleSelectorBase* pidSel=0;
977 
978  if (i%2) {
979  fListDefiners[i]->fList.Select(chargedCands,minusSel);
980  } else {
981  fListDefiners[i]->fList.Select(chargedCands,plusSel);
982  }
983 
984  switch (pidhint) {
985  case 0 :
986  pidSel = eSel;
987  break;
988  case 1 :
989  pidSel = muSel;
990  break;
991  case 2 :
992  pidSel = piSel;
993  break;
994  case 3 :
995  pidSel = kSel;
996  break;
997  case 4 :
998  pidSel = pSel;
999  break;
1000  }
1001 
1002  switch (crithint) {
1003  case 0 :
1004  pidSel->SetCriterion(veryLoose);
1005  break;
1006  case 1 :
1007  pidSel->SetCriterion(loose);
1008  break;
1009  case 2 :
1010  pidSel->SetCriterion(tight);
1011  break;
1012  case 3 :
1013  pidSel->SetCriterion(veryTight);
1014  break;
1015  }
1016 
1017  fListDefiners[i]->fList.Select(pidSel);
1018  }
1019 
1020  fListDefiners[40]->fList=neutralCands;
1021  fListDefiners[41]->fList=mcCands;
1022 
1023 
1024  return;
1025 }
1026 
1027 // -------------------------------------------------------------------------
1028 
1030 {
1031  cout<<"setupAnalysis"<<endl;
1032  std::ifstream cfgFile(fCfgFileName.c_str(),ios::in);
1033 
1034  int nLists=fListDefiners.size();
1035 
1036  unsigned int i=0,j=0;
1037 
1038  std::vector<int> daupdgs;
1039  unsigned int daucnt=0;
1040  unsigned int daulistcnt=0;
1041 
1042  char line[200];
1043  unsigned int linecnt=0;
1044 
1045  //bool inList=false;
1046  bool hasAnti=false;
1047 
1048  bool defineOpen=false;
1049  bool decaySet=false;
1050  //bool dauListsSet=false; //[R.K.03/2017] unused variable
1051 
1052  PndListDefiner* currentList;
1053  PndListDefiner* currentAntiList;
1054 
1055  // **************************************************
1056  // the big loop through the .cfg file
1057  // **************************************************
1058  while (!cfgFile.eof()) {
1059  cfgFile.getline(line,200);
1060  linecnt++;
1061 
1062  ArgVector tokenVec;
1063  CStrTok tokenizer;
1064 
1065  char* token = tokenizer.GetFirst(line, " \t");
1066 
1067  while(token) {
1068  tokenVec.push_back(token);
1069  token=tokenizer.GetNext(" \t");
1070  }
1071 
1072  if (tokenVec.size()==0) { continue; }
1073 
1074 
1075  // **************************************************
1076  // ***** BEGIN LIST DEFINITION **** Token = DefineList
1077  // **************************************************
1078 
1079  if (tokenVec[0]=="DefineList") {
1080  // are we already in a DefineList?
1081  if (defineOpen) { return ErrorMessage(100,linecnt); }
1082 
1083  // list already defined?
1084  if (fListMap.find(tokenVec[1])!=fListMap.end()) { return ErrorMessage(101,linecnt,tokenVec[1]); }
1085 
1086  currentList=new PndListDefiner(tokenVec[1]);
1087  currentList->fIsUsed=true;
1088 
1089  fListMap[tokenVec[1]]=nLists++;
1090 
1091  defineOpen=true;
1092  decaySet=false;
1093  //dauListsSet=false; //[R.K.03/2017] unused variable
1094  daupdgs.clear();
1095  daucnt=0;
1096  daulistcnt=0;
1097 
1098  } // ******** DefineList *********** END
1099 
1100 
1101 
1102  // **************************************************
1103  // ***** READ IN THE DECAY MODE **** Token = decayMode
1104  // **************************************************
1105 
1106  if (tokenVec[0]=="decayMode") {
1107  // are we in a DefineList statement?
1108  if (!defineOpen) { return ErrorMessage(200,linecnt); }
1109 
1110  // decaymode already defined?
1111  if (decaySet) { return ErrorMessage(201,linecnt); }
1112 
1113  currentList->fPdgCode=GetPdgCode(tokenVec[1]);
1114 
1115  for (i=3; i<tokenVec.size(); i++) {
1116  daupdgs.push_back(GetPdgCode(tokenVec[i]));
1117  }
1118 
1119  daucnt=daupdgs.size();
1120  decaySet=true;
1121 
1122  if (daucnt>5) { return ErrorMessage(203,linecnt); }
1123 
1124  } // ******* decayMode *********** END
1125 
1126 
1127 
1128  // **************************************************
1129  // ***** SET DAUGHTER LISTS **** Token = daughterList
1130  // **************************************************
1131 
1132  if (tokenVec[0]=="daughterList") {
1133  // are we in a DefineList statement?
1134  if (!defineOpen) { return ErrorMessage(200,linecnt); }
1135 
1136  // decaymode already defined?
1137  if (!decaySet) { return ErrorMessage(301,linecnt); }
1138 
1139  // more daughterLists than daughters in decaymode definition?
1140  if (daulistcnt>daucnt) { return ErrorMessage(302,linecnt); }
1141 
1142  // cannot use McTruth list for combinatorics
1143  if (tokenVec[1]=="McTruth") { return ErrorMessage(305,linecnt); }
1144 
1145  std::string dName=tokenVec[1];
1146  bool isgeneric=false;
1147 
1148  // the generic charged lists have a 'P'=plus or 'M'=minus as last letter
1149  // we have to modifiy accordingly
1150  if (IsGenericListName(dName)) {
1151  isgeneric=true;
1152 
1153  // replace 'Charged' with the corresponding particle list
1154  if (dName=="Charged")
1155  switch (abs(daupdgs[daulistcnt])) {
1156  case 11:
1157  dName="ElectronVeryLoose";
1158  break;
1159  case 13:
1160  dName="MuonVeryLoose";
1161  break;
1162  case 211:
1163  dName="PionVeryLoose";
1164  break;
1165  case 321:
1166  dName="KaonVeryLoose";
1167  break;
1168  case 2212:
1169  dName="ProtonVeryLoose";
1170  break;
1171  }
1172 
1173  float charge=TDatabasePDG::Instance()->GetParticle(daupdgs[daulistcnt])->Charge();
1174  if (charge<-0.1) { dName+="M"; }
1175  else if (charge>0.1) { dName+="P"; }
1176  currentList->fDauIdx.push_back(fListMap[dName]);
1177  fListDefiners[fListMap[dName]]->fIsUsed=true;
1178  }
1179 
1180  // here we have to check whether the list or corresponding anti list is meant
1181  if (!isgeneric) {
1182  // does list exist?
1183  if (fListMap.find(dName)==fListMap.end()) { return ErrorMessage(303,linecnt,dName); }
1184 
1185  int idx=fListMap[dName];
1186 
1187  int decpdg=daupdgs[daulistcnt];
1188  int listpdg=fListDefiners[idx]->fPdgCode;
1189  int alistpdg=fListDefiners[fListDefiners[idx]->fAntiIdx]->fPdgCode;
1190 
1191  // do particle types match between decayMode and daughterList definition?
1192  if (decpdg!=listpdg && decpdg!=alistpdg) { return ErrorMessage(304,linecnt); }
1193 
1194  // decide whether the list or the charged conjugate has to be used
1195  // be aware, that for pdg==0 (unspecified particles) always the list is used;
1196  // I don't know how to decide better at the moment
1197  if (decpdg==listpdg) {
1198  currentList->fDauIdx.push_back(idx);
1199  } else if (decpdg==alistpdg) {
1200  currentList->fDauIdx.push_back(fListDefiners[idx]->fAntiIdx);
1201  }
1202  }
1203 
1204  daulistcnt++;
1205  } // ******* daughterList *********** END
1206 
1207 
1208 
1209  // **************************************************
1210  // ***** END LIST DEFINITION **** Token = End
1211  // **************************************************
1212 
1213  if (tokenVec[0]=="End" || tokenVec[0]=="EndDefine") {
1214  // are we in a DefineList statement?
1215  if (!defineOpen) { return ErrorMessage(200,linecnt); }
1216  // same number of daughterLists as daughters in decaymode definition?
1217  if (!decaySet || daucnt<2 || daulistcnt!=daucnt) { return ErrorMessage(400,linecnt); }
1218 
1219  // everything seems to be ok, take care whether antilist is needed or not
1220  // therefore check whether the final state can be distinguished from it's c.c.
1221  // i.e. the set of list indices from the list of the anti-indices
1222 
1223  //sort the daughter index list
1224  //std::sort(currentList->fDauIdx.begin(),currentList->fDauIdx.end());
1225 
1226  fListDefiners.push_back(currentList);
1227 
1228  std::vector<int> finstate;
1229  std::vector<int> ccstate;
1230 
1231  double charge=0.0;
1232 
1233  for (i=0; i<daulistcnt; i++) {
1234  int idx = currentList->fDauIdx[i];
1235  int ccidx = fListDefiners[idx]->fAntiIdx;
1236 
1237  finstate.push_back(idx);
1238  ccstate.push_back(ccidx);
1239  charge+=fListDefiners[idx]->fCharge;
1240  }
1241 
1242  currentList->fCharge=charge;
1243 
1244  // to find out whether the final state and the cc final state are
1245  // different sort the lists of indices and compare pair-wise
1246  std::sort(finstate.begin(),finstate.end());
1247  std::sort(ccstate.begin(),ccstate.end());
1248 
1249  hasAnti=false;
1250 
1251  cout <<currentList->fName<<" : list indices of list and cc : ";
1252  for (i=0; i<finstate.size(); i++) {
1253  cout <<"( "<<finstate[i]<<" , "<<ccstate[i]<<" ) ";
1254  if (finstate[i]!=ccstate[i]) { hasAnti=true; }
1255  }
1256  cout <<endl;
1257 
1258  if (!hasAnti) { // set the idx of the antilist to own idx
1259  currentList->fAntiIdx=fListMap[currentList->fName];
1260  } else {
1261  // create the c.c. list
1262  std::string aName=currentList->fName+"_A";
1263  currentAntiList=new PndListDefiner(aName);
1264  fListMap[aName]=nLists++;
1265 
1266  // set the idx for both lists to each other
1267  currentList->fAntiIdx=fListMap[aName];
1268  currentAntiList->fAntiIdx=fListMap[currentList->fName];
1269 
1270  charge=0.;
1271  // set the daughterlist idx
1272  for (i=0; i<ccstate.size(); i++) {
1273  int idx = currentList->fDauIdx[i];
1274  int ccidx = fListDefiners[idx]->fAntiIdx;
1275 
1276  currentAntiList->fDauIdx.push_back(ccidx);
1277 
1278  //also sum charges of daughters of anti list
1279  charge+=fListDefiners[ccstate[i]]->fCharge;
1280 
1281  // all the lists used as input for antilist have to be filled as well
1282  fListDefiners[ccstate[i]]->fIsUsed=true;
1283  }
1284 
1285  // set the properties of Anti-List according to list properties
1286  currentAntiList->fPdgCode=GetAntiPdgCode(currentList->fPdgCode);
1287  currentAntiList->fCharge=charge;
1288  currentAntiList->fIsAntiList=true;
1289  currentAntiList->fIsUsed=true;
1290  currentAntiList->fHisto=currentList->fHisto;
1291  currentAntiList->fSelector=currentList->fSelector;
1292  currentAntiList->fDumpList=currentList->fDumpList;
1293 
1294  fListDefiners.push_back(currentAntiList);
1295  }
1296 
1297 
1298  defineOpen=false;
1299 
1300  }// ******* End *********** END
1301 
1302 
1303  // **************************************************
1304  // ***** Misc Command Definitions
1305  // **************************************************
1306 
1307  // set the initial 4 vector i.e. pbar momentum or E_cms
1308  if (tokenVec[0]=="pbarMom") {
1309  if (fpInit.E()!=0) { ErrorMessage(503,linecnt); }
1310 
1311  double p = atof(tokenVec[1].c_str());
1312  double mp = 0.938272;
1313  double E = sqrt(mp*mp+p*p)+mp;
1314 
1315  fpInit.SetXYZT(0.0,0.0,p,E);
1316  }
1317 
1318  if (tokenVec[0]=="Ecms") {
1319  if (fpInit.E()!=0) { ErrorMessage(503,linecnt); }
1320 
1321  double M = atof(tokenVec[1].c_str());
1322  double mp = 0.938272;
1323 
1324  double X = (M*M-2*mp*mp)/(2*mp);
1325  double p = sqrt(X*X-mp*mp);
1326  double E = sqrt(mp*mp+p*p)+mp;
1327 
1328  fpInit.SetXYZT(0.0,0.0,p,E);
1329  }
1330 
1331  // set maximum for entries in column
1332 
1333  /* if (tokenVec[0]=="MaxEntries")
1334  {
1335  int num=atoi(tokenVec[1].c_str());
1336  if (num<=0 || num>32768)
1337  {
1338  ErrorMessage(510,linecnt);
1339  }
1340  else
1341  fMaxEntries=num;
1342  }*/
1343 
1344 
1345 
1346  // **************************************************
1347  // dump the list into nTuple
1348  // **************************************************
1349 
1350  if (tokenVec[0]=="dumpList") {
1351  // are we in a DefineList statement? Then we need at least table prefix
1352  if (defineOpen && tokenVec.size()<2) { return ErrorMessage(501,linecnt); }
1353 
1354  // are we outside of DefineList statement? We need in addition list name!
1355  if (!defineOpen && tokenVec.size()<3) { return ErrorMessage(502,linecnt); }
1356 
1357  // tokenVec[startcols] is the first column name to dump
1358  // in DefineList this is position 2; outside it is pos 3
1359  unsigned int startcols=2;
1360 
1361  std::string dName=tokenVec[1];
1362 
1363  if (!defineOpen) {
1364  // List has to be defined before
1365  //for Generic Names like PionVeryLoose add the "P" at the end
1366  if (dName=="Charged") { dName="PionVeryLoose"; }
1367 
1368  if (IsGenericListName(dName) && dName!="Neutral") { dName+="P"; }
1369 
1370  if (fListMap.find(dName)==fListMap.end()) { return ErrorMessage(303,linecnt,dName); }
1371 
1372  // which list is it?
1373  int idx=fListMap[dName];
1374 
1375  currentList=fListDefiners[idx];
1376 
1377  // we also want to dump Anti list if existing!
1378  fListDefiners[currentList->fAntiIdx]->fDumpList=true;
1379  fListDefiners[currentList->fAntiIdx]->fIsUsed=true;
1380 
1381  startcols=3;
1382  }
1383 
1384  currentList->fDumpList=true;
1385  currentList->fIsUsed=true;
1386  currentList->fColName=tokenVec[startcols-1];
1387 
1388  if (tokenVec.size()==startcols) {
1389  tokenVec.push_back("p4");
1390  tokenVec.push_back("ch");
1391  tokenVec.push_back("m");
1392  }
1393 
1394  for (i=startcols; i<tokenVec.size(); i++) {
1395  std::string col=tokenVec[i];
1396 
1397  if (fColKeyMap.find(col)==fColKeyMap.end()) {
1398  cout << "-W- PndSimpleAnalysis::SetupAnalysis: In '"<< fCfgFileName
1399  << "', line "<<linecnt<<": ";
1400  cout << "Unknown column name '"<<col<<"'."<<endl;
1401 
1402  continue;
1403  }
1404 
1405  int key=fColKeyMap[col];
1406 
1407  // 'mothi' and 'pdg' cannot be added explicit
1408  if (key==520 || key==521) { continue; }
1409 
1410  // daughter masses cannot be added explicit
1411  if (key>=150 && key<=154) { continue; }
1412 
1413  // daughter list indices cannot be added explicit
1414  if (key>=501 && key<=505) { continue; }
1415 
1416  // short cut!
1417  if (key>=1000) {
1418  std::vector<std::string> keylist=fColShortKeyMap[key];
1419  for (j=0; j<keylist.size(); j++) { tokenVec.push_back(keylist[j]); }
1420 
1421  } else if (key<500) { currentList->fNtpFNames.push_back(col); }
1422  else { currentList->fNtpINames.push_back(col); }
1423 
1424  }
1425 
1426  if (dName=="McTruth") {
1427  currentList->fNtpINames.push_back("pdg");
1428  currentList->fNtpINames.push_back("mothi");
1429  }
1430 
1431  } // dumpList
1432 
1433 
1434  // ***** WRITE OUT MASS HISTOGRAM **** Token = histogram
1435 
1436  if (tokenVec[0]=="histogram") {
1437  // are we in a DefineList statement?
1438  if (!defineOpen) { return ErrorMessage(200,linecnt); }
1439 
1440  double width=0.1;
1441  double mass=1.0;
1442 
1443  if (currentList->fPdgCode) {
1444  mass=TDatabasePDG::Instance()->GetParticle(currentList->fPdgCode)->Mass();
1445  }
1446 
1447  if (tokenVec.size()==2) {
1448  width=atof(tokenVec[1].c_str());
1449  }
1450 
1451  double min=mass - width;
1452  double max=mass + width;
1453 
1454  if (tokenVec.size()==3) {
1455  double num1=atof(tokenVec[1].c_str());
1456  double num2=atof(tokenVec[2].c_str());
1457  if (num1<num2) {min=num1; max=num2;}
1458  else {min=num2; max=num1;}
1459  }
1460 
1461  std::string hname=currentList->fName+"_M";
1462  std::string htitle=currentList->fName+" mass";
1463 
1464  TH1F* h=new TH1F(hname.c_str(),htitle.c_str(),100,min,max);
1465  currentList->fHisto.push_back(h);
1466 
1467  } // histogram
1468 
1469  // ***** SELECTOR DEFINITIONS **** Token = select
1470 
1471  if (tokenVec[0]=="select") {
1472  // add a mass selection
1473 
1474  if (tokenVec[1]=="Mass") {
1475  // particle unknown (pdg==0) and only one argument given?
1476  if (currentList->fPdgCode==0 && tokenVec.size()<4) {
1477  return ErrorMessage(500,linecnt, currentList->fName);
1478  }
1479 
1480  double mean=TDatabasePDG::Instance()->GetParticle(currentList->fPdgCode)->Mass();
1481  double width=0.1;
1482 
1483  // parameter specifies width around nominal mass
1484  if (tokenVec.size()==3) {
1485  width=atof(tokenVec[2].c_str());
1486  }
1487 
1488  if (tokenVec.size()>3) {
1489  mean=atof(tokenVec[2].c_str());
1490  width=atof(tokenVec[3].c_str());
1491  }
1492  RhoMassParticleSelector* sel=new RhoMassParticleSelector((currentList->fName+"sel").c_str(),mean,width);
1493  currentList->fSelector.push_back(sel);
1494  } //Mass
1495 
1496  } //selector
1497 
1498 
1499 
1500  } // while cfg file
1501 
1502  cfgFile.close();
1503 
1504  std::map<std::string, int> colMap;
1505 
1506  for (j=0; j<fListDefiners.size(); j++) {
1507  currentList=fListDefiners[j];
1508 
1509  if (!currentList->fIsAntiList && currentList->fDumpList) {
1510  if (0==ntp) { ntp=new TTree("ntp","PndSimpleAnalysis NTuple"); }
1511  std::string pre=currentList->fColName;
1512 
1513  int nd=currentList->GetNDau();
1514 
1515  // add the daughter idx branches
1516  if (nd>0) { currentList->fNtpINames.push_back("d1"); }
1517  if (nd>1) { currentList->fNtpINames.push_back("d2"); }
1518  if (nd>2) { currentList->fNtpINames.push_back("d3"); }
1519  if (nd>3) { currentList->fNtpINames.push_back("d4"); }
1520  if (nd>4) { currentList->fNtpINames.push_back("d5"); }
1521 
1522  for (i=0; i<currentList->fNtpFNames.size(); i++) {
1523  if (currentList->fNtpFNames[i]=="dm") {
1524  currentList->fNtpFNames.erase(currentList->fNtpFNames.begin()+i);
1525 
1526  if (nd>0) { currentList->fNtpFNames.push_back("d1m"); }
1527  if (nd>1) { currentList->fNtpFNames.push_back("d2m"); }
1528  if (nd>2) { currentList->fNtpFNames.push_back("d3m"); }
1529  if (nd>3) { currentList->fNtpFNames.push_back("d4m"); }
1530  if (nd>4) { currentList->fNtpFNames.push_back("d5m"); }
1531 
1532  i=1000;
1533  }
1534  }
1535 
1536  std::string brname="n"+pre;
1537  if (colMap.find(brname)!=colMap.end()) { return ErrorMessage(504,0,brname); }
1538  colMap[brname]=1;
1539 
1540  //create the branch with the number of entries
1541  ntp->Branch(brname.c_str(),&(currentList->fNEntries),(brname+"/I").c_str());
1542 
1543  //create the branches holding the values (float)
1544  for (i=0; i<currentList->fNtpFNames.size(); i++) {
1545  brname=pre+currentList->fNtpFNames[i];
1546  cout <<brname<<" ";
1547  if (colMap.find(brname)!=colMap.end()) { return ErrorMessage(504,0,brname); }
1548  colMap[brname]=1;
1549 
1550  float* theAr=new float[fMaxEntries];
1551 
1552  currentList->fNtpFArrays.push_back(theAr);
1553  ntp->Branch(brname.c_str(),currentList->fNtpFArrays[i],(brname+"[n"+pre+"]/F").c_str());
1554  }
1555 
1556  //create the branches holding the values (int, at the moment only the daughter indices)
1557  for (i=0; i<currentList->fNtpINames.size(); i++) {
1558  brname=pre+currentList->fNtpINames[i];
1559  cout <<brname<<" ";
1560  if (colMap.find(brname)!=colMap.end()) { return ErrorMessage(504,0,brname); }
1561  colMap[brname]=1;
1562 
1563  int* theAr=new int[fMaxEntries];
1564  currentList->fNtpIArrays.push_back(theAr);
1565  ntp->Branch(brname.c_str(), currentList->fNtpIArrays[i], (brname+"[n"+pre+"]/I").c_str());
1566  }
1567  cout <<endl;
1568 
1569 
1570  } // dump List
1571 
1572  }
1573 
1574 
1575 
1576  for (i=0; i<fListDefiners.size(); i++)
1577  if (fListDefiners[i]->fIsUsed) {cout <<i<<" : "; fListDefiners[i]->Print(); }
1578 
1579  //if (ntp) ntp->Print();
1580 
1581  return true;
1582 }
1583 
1584 
1585 // -------------------------------------------------------------------------
1586 
1588 {
1589  for (unsigned int i=0; i<fListDefiners.size(); i++) {
1591  if (cur->fIsAntiList) { continue; }
1592  for (int j=0; j<cur->GetNHistos(); j++) {
1593  cur->fHisto[j]->Write();
1594  }
1595  }
1596 
1597  if (ntp) {
1598  ntp->Write();
1599  }
1600 }
1601 
1602 // -------------------------------------------------------------------------
1603 
1605 {
1606  int code=0;
1607  if (TDatabasePDG::Instance()->GetParticle(name.c_str())) {
1608  code=TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
1609  }
1610 
1611  //cout <<name<<" "<<code<<endl;
1612 
1613  return code;
1614 }
1615 
1616 // -------------------------------------------------------------------------
1617 
1619 {
1620  int code=0, pdgcode=0;
1621  if (TDatabasePDG::Instance()->GetParticle(name.c_str())) {
1622  pdgcode=TDatabasePDG::Instance()->GetParticle(name.c_str())->PdgCode();
1623  if (TDatabasePDG::Instance()->GetParticle(pdgcode)->AntiParticle()) {
1624  code=-pdgcode;
1625  } else {
1626  code=pdgcode;
1627  }
1628  }
1629 
1630  //cout <<name<<"("<<pdgcode<<") A="<<code<<endl;
1631 
1632  return code;
1633 }
1634 
1635 // -------------------------------------------------------------------------
1636 
1638 {
1639  int code=0;
1640  if (TDatabasePDG::Instance()->GetParticle(pdgcode)) {
1641  if (TDatabasePDG::Instance()->GetParticle(pdgcode)->AntiParticle()) {
1642  code=-pdgcode;
1643  } else {
1644  code=pdgcode;
1645  }
1646  }
1647 
1648  //cout <<pdgcode<<" A="<<code<<endl;
1649 
1650  return code;
1651 }
1652 
1653 // -------------------------------------------------------------------------
1654 
1655 bool PndSimpleAnalysis::ErrorMessage(int mid, int line, std::string arg)
1656 {
1657  cout << "-E- PndSimpleAnalysis::SetupAnalysis: In '"<< fCfgFileName << "', line "<<line<<": ";
1658 
1659  switch (mid) {
1660  case 100:
1661  cout << "Missing 'End' statement!";
1662  break;
1663 
1664  case 101:
1665  cout << "List with name '"<< arg<< "' already defined";
1666  break;
1667 
1668  case 200:
1669  cout << "Missing 'DefineList' statement!";
1670  break;
1671 
1672  case 201:
1673  cout << "DecayMode already defined!";
1674  break;
1675 
1676  case 203:
1677  cout << "More than 5 daughters not allowed!";
1678  break;
1679 
1680  case 301:
1681  cout << "Missing 'decayMode' statement!";
1682  break;
1683 
1684  case 302:
1685  cout << "Too many daughterLists defined!";
1686  break;
1687 
1688  case 303:
1689  cout << "Unknown List '"<<arg<<"'!";
1690  break;
1691 
1692  case 304:
1693  cout << "Particle type mismatch between daughterList and decayMode!";
1694  break;
1695 
1696  case 305:
1697  cout << "McTruth cannot be used as daughterList!";
1698  break;
1699 
1700  case 400:
1701  cout << "Wrong number of daughters or daughterLists!";
1702  break;
1703 
1704  case 500:
1705  cout << "Mass of particle '"<<arg <<"' unknown; must specify selector mean value!";
1706  break;
1707 
1708  case 501:
1709  cout << "Must specify name of tuple column!";
1710  break;
1711 
1712  case 502:
1713  cout << "Must specify name of list and tuple column!";
1714  break;
1715 
1716  case 503:
1717  cout << "Warning: Initial 4 vector already set before!";
1718  break;
1719 
1720  case 504:
1721  cout << "Branch '"<<arg<<"' already exists in NTuple!";
1722  break;
1723 
1724  case 510:
1725  cout << "Invalid number for parameter 'MaxEntries'.";
1726  break;
1727 
1728  default:
1729  cout << "Unexpected error!";
1730  break;
1731 
1732  }
1733 
1734  cout <<endl;
1735 
1736  return false;
1737 }
1738 
1739 
1740 // -------------------------------------------------------------------------
1741 
1743 
1744 
1745 /*
1746 // shortcut 4 vector
1747 if (col=="p4")
1748 {
1749  currentList->fNtpFNames.push_back("px");
1750  currentList->fNtpFNames.push_back("py");
1751  currentList->fNtpFNames.push_back("pz");
1752  currentList->fNtpFNames.push_back("en");
1753 }
1754 else if (col=="p4cm")
1755 {
1756  currentList->fNtpFNames.push_back("pxcm");
1757  currentList->fNtpFNames.push_back("pycm");
1758  currentList->fNtpFNames.push_back("Cmpzcm");
1759  currentList->fNtpFNames.push_back("encm");
1760 }
1761 // shortcut position
1762 else if (col=="pos")
1763 {
1764  currentList->fNtpFNames.push_back("vx");
1765  currentList->fNtpFNames.push_back("vy");
1766  currentList->fNtpFNames.push_back("vz");
1767 }
1768 else if (col=="p3")
1769 {
1770  currentList->fNtpFNames.push_back("phi");
1771  currentList->fNtpFNames.push_back("tht");
1772  currentList->fNtpFNames.push_back("cth");
1773  currentList->fNtpFNames.push_back("p");
1774 }
1775 else if (col=="p3cm")
1776 {
1777  currentList->fNtpFNames.push_back("phicm");
1778  currentList->fNtpFNames.push_back("thtcm");
1779  currentList->fNtpFNames.push_back("cthcm");
1780  currentList->fNtpFNames.push_back("pcm");
1781 }
1782 // add all known columns
1783 else
1784 {
1785  // unknown name; skip entry
1786  if (fColKeyMap.find(col)==fColKeyMap.end())
1787  {
1788  cout << "-W- PndSimpleAnalysis::SetupAnalysis: In '"<< fCfgFileName
1789  << "', line "<<linecnt<<": ";
1790  cout << "Unknown column name '"<<col<<"'."<<endl;
1791 
1792  continue;
1793  }
1794 
1795  int key=fColKeyMap[col];
1796 
1797  if (key<200) currentList->fNtpFNames.push_back(col);
1798  else currentList->fNtpINames.push_back(col);
1799 }
1800 */
1801 
int GetPdgCode(std::string name)
RhoCandList neutralCands
void Add(const RhoCandidate *c)
Definition: RhoCandList.h:49
Double_t p
Definition: anasim.C:58
RhoMinusParticleSelector * minusSel
Double_t GetProtonPidProb(PndPidProbability *flux=NULL) const
Int_t run
Definition: autocutx.C:47
void Cleanup()
Definition: RhoCandList.cxx:62
void PrintTree(RhoCandidate *tc, int level=0)
int GetAntiPdgCode(std::string name)
Double_t GetKaonPidProb(PndPidProbability *flux=NULL) const
std::vector< TH1F * > fHisto
Int_t i
Definition: run_full.C:25
std::vector< std::string > fNtpINames
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
Float_t GetSttMeanDEDX() const
std::vector< std::string > fGenericListNames
Int_t GetLength() const
Definition: RhoCandList.h:46
void Boost(const TVector3 &)
int col
Definition: anaLmdDigi.C:67
TClonesArray * fNeutralProbability
static void Reset()
Definition: RhoFactory.cxx:28
#define fMaxEntries
int n
std::string fColName
RhoSimpleElectronSelector * eSel
TClonesArray * fChargedArray
virtual Bool_t Accept(RhoCandidate *)=0
RhoCandidate * Daughter(Int_t n)
RhoSimplePionSelector * piSel
cout<< "POINTs for new FwEndCap == "<< tsim-> GetEntriesFast()
friend F32vec4 max(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:26
TClonesArray * fMicroArray
Float_t GetRichThetaC() const
Int_t Uid() const
Definition: RhoCandidate.h:419
TLorentzVector fpInit
bool IsGenericListName(std::string n)
std::vector< float * > fNtpFArrays
static const double mp
Definition: mzparameters.h:11
Float_t GetDrcThetaC() const
std::string fCfgFileName
Float_t GetDiscThetaC() const
int idx[MAX]
Definition: autocutx.C:38
int uid(int lev, int lrun, int lmode)
Definition: autocutx.C:122
void Combine(RhoCandList &l1, RhoCandList &l2)
void SetPidInfo(double *pidinfo=0)
std::string fName
RhoCandList fList
Double_t GetMuonPidProb(PndPidProbability *flux=NULL) const
TClonesArray * fChargedProbability
std::map< int, std::vector< std::string > > fColShortKeyMap
char * GetFirst(char *lpsz, const char *lpcszDelimiters)
Definition: StrTok.cxx:29
void Select(RhoParticleSelectorBase *pidmgr)
Float_t GetMvdDEDX() const
std::vector< int > fDauIdx
friend F32vec4 min(const F32vec4 &a, const F32vec4 &b)
Definition: P4_F32vec4.h:25
static RhoFactory * Instance()
Definition: RhoFactory.cxx:34
void SetConfigFile(std::string filename="analysis.cfg")
bool ErrorMessage(int mid, int line=0, std::string arg="")
TString name
std::map< std::string, int > fColKeyMap
RhoSimpleMuonSelector * muSel
double X
Definition: anaLmdDigi.C:68
virtual void SetCriterion(const char *crit)
TClonesArray * fMcArray
std::vector< int * > fNtpIArrays
Definition: StrTok.h:11
Int_t NDaughters() const
ClassImp(PndAnaContFact)
RhoPlusParticleSelector * plusSel
RhoCandList chargedCands
RhoSimpleKaonSelector * kSel
RhoSimpleProtonSelector * pSel
Double_t mean[nsteps]
Definition: dedx_bands.C:65
Bool_t IsCloneOf(const RhoCandidate &, Bool_t checkType=kFALSE) const
virtual InitStatus Init()
char * GetNext(const char *lpcszDelimiters)
Definition: StrTok.cxx:37
std::vector< std::string > ArgVector
Definition: ArgList.h:8
virtual void Exec(Option_t *opt)
std::vector< std::string > fNtpFNames
virtual void SetParContainers()
std::vector< RhoParticleSelectorBase * > fSelector
Double_t GetElectronPidProb(PndPidProbability *flux=NULL) const
Double_t GetPionPidProb(PndPidProbability *flux=NULL) const
Float_t GetTofM2() const
double GetPidInfo(int hypo)
TClonesArray * fNeutralArray
std::map< std::string, int > fListMap
const string filename
std::vector< PndListDefiner * > fListDefiners