FairRoot/PandaRoot
softrig/PndSoftTriggerTask.cxx
Go to the documentation of this file.
1 // ************************************************************************
2 //
3 // Online SoftTrigger reco and tagging for several channels
4 //
5 // K.Goetzen 10/2013
6 // ************************************************************************
7 
8 
9 // The header file
10 #include "PndSoftTriggerTask.h"
11 
12 // Package headers
13 #include "PndOnlineFilterInfo.h"
14 #include "PndSoftTriggerLine.h"
15 
16 // C++ headers
17 #include <string>
18 #include <iostream>
19 #include <fstream>
20 #include <map>
21 #include <stdlib.h> /* getenv */
22 
23 // FAIR headers
24 #include "FairRootManager.h"
25 #include "FairRunAna.h"
26 #include "FairRuntimeDb.h"
27 #include "FairRun.h"
28 #include "FairRuntimeDb.h"
29 
30 // ROOT headers
31 #include "TClonesArray.h"
32 #include "TVector3.h"
33 #include "TH1F.h"
34 #include "TH2F.h"
35 #include "TString.h"
36 #include "TRegexp.h"
37 #include "TMVA/Reader.h"
38 
39 // RHO headers
40 #include "RhoCandidate.h"
41 #include "RhoHistogram/RhoTuple.h"
42 #include "RhoFactory.h"
46 #include "RhoTuple.h"
47 
48 // Analysis headers
49 #include "PndAnalysis.h"
50 #include "Rho4CFitter.h"
51 #include "RhoKinVtxFitter.h"
52 #include "RhoKalmanVtxFitter.h"
53 #include "RhoKinFitter.h"
54 #include "RhoVtxPoca.h"
55 #include "PndPidCandidate.h"
56 #include "PndEventShape.h"
57 #include "PndRhoTupleQA.h"
58 
59 const int STMAXCUT = 50;
60 const int STMAXEVVARS = 50;
61 
62 using std::cout;
63 using std::endl;
64 
65 // *** Holds the cut set for a certain mode
66 struct STCutSet
67 {
68  int ncut;
70  int op[STMAXCUT];
71  double cutval[STMAXCUT];
72 
73  TMVA::Reader *tmvaread;
74  int tmvanvar;
76  float tmvacut;
78 };
79 
80 // *** this is the array of input variables for the cut selection or TMVA
82 
83 // *** maps needed for selection parsing
84 std::map<TString, int> fSTVarmap;
85 std::map<int, STCutSet> fSTSelmap;
86 //std::map<int, TString> fSTOps;
87 TString fSTOps[5]={">=",">","==","<=","<"};
88 std::map<int, PndSoftTriggerLine*> fSTTriggers;
89 std::map<int, int> fSTMapListIndex;
90 
91 typedef std::map<int, PndSoftTriggerLine*>::iterator TrigIt;
92 
93 // e+ e- mu+ mu- pi+ pi- K+ K- p pb gam pi0 KS eta aux anti-aux
94 int fSTPidIndex[16] = {-11, 11, -13, 13, 211, -211, 321, -321, 2212, -2212, 22, 111, 310, 221, 0, 0};
95 
96 // *** Event Shape Vars
97 
98 // ( 0) eslnpide ( 1) eslnpidmu ( 2) eslnpidpi ( 3) eslnpidk ( 4) eslnpidp
99 // ( 5) esfw1 ( 6) esfw2 ( 7) esfw3 ( 8) esfw4 ( 9) esfw5
100 // ( 10) esapl ( 11) escir ( 12) espla ( 13) esnneut ( 14) esnpart
101 // ( 15) esnchrg ( 16) espmax ( 17) espmaxl ( 18) espmin ( 19) espminl
102 // ( 20) esprapmax ( 21) esptmax ( 22) essumen ( 23) essumenl ( 24) essumetn
103 // ( 25) essumpc ( 26) essumpc05 ( 27) essumpcl ( 28) essumpt ( 29) essumptc
104 // ( 30) esthr ( 31) essph ( 32) esptmin ( 33) esncp10l ( 34) esnne10l
105 
106 int fSTNEvVars = 35;
108  "eslnpide" , "eslnpidmu" , "eslnpidpi" , "eslnpidk" , "eslnpidp" ,
109  "esfw1" , "esfw2" , "esfw3" , "esfw4" , "esfw5" ,
110  "esapl" , "escir" , "espla" , "esnneut" , "esnpart" ,
111  "esnchrg" , "espmax" , "espmaxl" , "espmin" , "espminl" ,
112  "esprapmax" , "esptmax" , "essumen" , "essumenl" , "essumetn" ,
113  "essumpc" , "essumpc05" , "essumpcl" , "essumpt" , "essumptc" ,
114  "esthr" , "essph" , "esptmin" , "esncp10l" , "esnne10l"
115 };
116 
117 // coding for variable suffix (to translate variable names to code)
118 int fSTNQuant = 17;
119 TString fSTVarSuff[] = {"pt","tht","pcm","thtcm","pide","pidmu","pidpi","pidk","pidp", "pidmax","pocqa","pocdist","pocctau","oang","cdecang","decang","p" };
120 int fSTQuantCode[] = { 11, 12, 13, 14, 20, 21, 22, 23, 24, 25, 30, 31, 32, 40, 42, 41, 10 };
121 
122 
123 // cache for event shape variables
125 
126 // cache map for candidate variables
127 std::map<int, double> fSTVarCandArray;
128 
129 // codes for the available energies
130 std::vector<int> fSTencode;
131 int fSTModeIndex = 0;
132 
133 // ----- Default constructor -------------------------------------------
134 PndSoftTriggerTask::PndSoftTriggerTask(double pmom, int mode, int runnum, TString trigfilename) :
135  FairTask("Panda Softtrigger Task"),
136  fVerbose(0), fMode(mode), fEvtCount(0), fRunNum(runnum), fSigCount(0), fNsigTag(8.0), fNsigAux(3.0), fDstMDiffCut(10.),
137  fTriggerFileName(trigfilename), fPhotosMax(0), fPhotosThresh(0.05),
138  fIniP4(0,0,0,0), fEcm(0.), fPbarMom(pmom),
139  fQAPi0(false),fQAEta(false),fQAKs0(false),fQAEvent(false), fQAMc(false), fQAMctOnly(false),
140  fGammaMinE(0.03), fPi0MinE(0.0), fEtaMinE(0.0), fTrackMinP(0.15), fIniPidCut(0.0),
141  fPi0Sel(NULL), fEtaSel(NULL), fKs0Sel(NULL),
142  fMomentumSel(NULL), fEnergySel(NULL),
143  ntp(0), nks0(0), npi0(0), neta(0), nmc(0),
144  fEventShape(NULL),
145  fQA(NULL)
146 {
147  fPdg = TDatabasePDG::Instance();
148 
149  // *** add several pbar p/n/dd Systems for MC truth match
150  double ppwidth = 0.01;
151  fPdg->AddParticle("pbarpSystem", "pbar p", fEcm, false, ppwidth,0,"",88888);
152  fPdg->AddParticle("pbarpSystem0","pbar p", fEcm, false, ppwidth,0,"",88880);
153  fPdg->AddParticle("pbarpSystem1","pbar p", fEcm, false, ppwidth,0,"",88881);
154  fPdg->AddParticle("pbarpSystem2","pbar p", fEcm, false, ppwidth,0,"",88882);
155 
156  fPdg->AddParticle("pbarnSystem", "pbar n", fEcm, false, ppwidth,0,"",88887);
157  fPdg->AddParticle("pbardSystem", "pbar d", fEcm, false, ppwidth,0,"",88889);
158 
159  double mp = fPdg->GetParticle("proton")->Mass(); //Proton mass for computation of p4_ini
160 
161  // set 4-vector of pbar-p-system
162  fIniP4.SetPz(pmom);
163  fIniP4.SetE(sqrt(pmom*pmom+mp*mp)+mp);
164  fEcm = fIniP4.M();
165 
166  // set default algorithms for pid
167  fAlgoElectron = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts";
172 
173  fIniPidCut = 0.0;
174 
175  fCfgFileName = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/selection_10ch_tight.cfg";
176  fApplyFullSelection = false;
177 
178  // *** set default signal parameters (mean, sigma)
181 
182  // if no trigger definition file is given, use the default one
183  if (fTriggerFileName.Length()==0)
184  {
185  fTriggerFileName = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/triggerlines.cfg";
186  cout <<"Reading default trigger lines file "<<fTriggerFileName.Data()<<endl;
187  }
188  // map pdg codes to RhoCandList index
189  for (int i=0;i<14;++i) fSTMapListIndex[fSTPidIndex[i]] = i;
190 
192 }
193 // -------------------------------------------------------------------------
194 
195 
196 // ----- Destructor ----------------------------------------------------
198 {
199  delete fPocaVertexer;
200 }
201 // -------------------------------------------------------------------------
202 
204 {
205  TString selectioncfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/selection_fsim_dec2014.cfg";
206 
207  SetConfigurationFile(selectioncfg);
208  ApplyFullSelection(1); // apply selection defined in 'TString selectioncfg'
209 
210  SetPi0SignalParams(0.134, 0.0035); // set parameters for pi0
211  SetKs0SignalParams(0.497, 0.0055); // set parameters for KS
212  SetEtaSignalParams(0.549, 0.0088); // set parameters for eta
213 
214  SetGammaMinE(0.10); // global energy pre-cut for neutrals
215  SetTrackMinP(0.10); // global momentum pre-cut for charged
216  SetInitialPidCut(0.1); // global PID pre-cut for charged
217  SetDstMDiffCut(0.1); // special cut on D*-D mass difference (to reduce comb and output file size)
218 
219  SetPidAlgoAll("PidChargedProbability");
220 
221  SetTagAll(true); // tag all modes
222  SetQAAll(false); // don't write any QA tuple
223 }
224 
225 // -------------------------------------------------------------------------
227 {
228  TString selectioncfg = TString(gSystem->Getenv("VMCWORKDIR"))+"/softrig/selection_full_jan2015.cfg";
229 
230  SetConfigurationFile(selectioncfg);
231  ApplyFullSelection(1); // apply selection defined in 'TString selectioncfg'
232 
233  SetPi0SignalParams(0.134, 0.0049);
234  SetKs0SignalParams(0.497, 0.0085);
235  SetEtaSignalParams(0.549, 0.0092);
236 
237  SetGammaMinE(0.10); // global energy pre-cut for neutrals
238  SetTrackMinP(0.10); // global momentum pre-cut for charged
239  SetInitialPidCut(0.1); // global PID pre-cut for charged
240  SetDstMDiffCut(0.1); // special cut on D*-D mass difference (to reduce comb and output file size)
241 
242  SetPidAlgoAll("PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts;PidAlgoMvd");
243 
244  SetTagAll(true); // tag all modes
245  SetQAAll(false); // don't write any QA tuple
246 }
247 
248 // -------------------------------------------------------------------------
249 // ----- Public method Init --------------------------------------------
251 {
252  fRootManager = FairRootManager::Instance();
253 
254  // Register TCA for tagging info
255  fTcaOnlineFilterInfo = new TClonesArray ( "PndOnlineFilterInfo" );
256  if (fRootManager)
257  fRootManager->Register ( "OnlineFilterInfo","PndOnlineFolder", fTcaOnlineFilterInfo, kTRUE );
258 
259  // *** initialize analysis object
260  fAnalysis = new PndAnalysis();
262 
263  // *** RhoTuple QA helper
265 
266  // *** Save current gDirectory
267  TDirectory *dir = gDirectory;
268  fRootManager->GetOutFile()->cd();
269 
270  // *** create ntuple
271  if (fQAEvent) { ntp = new RhoTuple("ntpev","Soft Trigger Common"); ntp->GetInternalTree()->SetDirectory(gDirectory);}
272  if (fQAKs0) { nks0 = new RhoTuple("nks0","K_S -> pi+ pi-"); nks0->GetInternalTree()->SetDirectory(gDirectory);}
273  if (fQAPi0) { npi0 = new RhoTuple("npi0","pi0 -> gamma gamma"); npi0->GetInternalTree()->SetDirectory(gDirectory);}
274  if (fQAEta) { neta = new RhoTuple("neta","eta -> gamma gamma"); neta->GetInternalTree()->SetDirectory(gDirectory);}
275  if (fQAMc) { nmc = new RhoTuple("nmc", "MC info"); nmc->GetInternalTree()->SetDirectory(gDirectory);}
276 
277 
278  // *** create mass pre selectors for QA (formular takes into account RhoSelector definition mean +- win/2
279  fPi0PreSel = new RhoMassParticleSelector("pi0PreSel", (fPi0QaMax + fPi0QaMin)/2.0, fPi0QaMax - fPi0QaMin );
280  fEtaPreSel = new RhoMassParticleSelector("etaPreSel", (fEtaQaMax + fEtaQaMin)/2.0, fEtaQaMax - fEtaQaMin );
282 
283  // *** number of sigmas deviation for tag
284  //double tagNumSig = fNsigTag;
285 
286  // *** create final selectors for pi0, eta, KS
288  fEtaSel = new RhoMassParticleSelector("etaSel", fEtaMean, fEtaSigma*2.0*fNsigAux);
289  fKs0Sel = new RhoMassParticleSelector("Ks0Sel", fKs0Mean, fKs0Sigma*2.0*fNsigAux);
290 
291  // *** basic selectors for preselection
292  fMomentumSel = new RhoMomentumParticleSelector("PSel",50.+fTrackMinP,100.);
293  fEnergySel = new RhoEnergyParticleSelector("ESel",50.+fGammaMinE,100.);
294 
295  // *** the poca vertexer
296  fPocaVertexer = new RhoVtxPoca();
297 
298  // *** read selection from configuration file
299  fSTencode.clear();
301 
302  // *** set mode index for current beam momentum
303  fSTModeIndex = 0;
304  double diff = 1000.;
305 
306  for (size_t i=0; i<fSTencode.size(); ++i)
307  {
308  double en = (double)fSTencode[i]/100.;
309  if (fabs(fEcm-en)<diff)
310  {
311  diff=fabs(fEcm-en);
312  fSTModeIndex = i;
313  }
314  }
315 
316  if (fApplyFullSelection && diff>1.0) cout <<"[PndSoftTriggerTask::Init] **** WARNING: Energy of best matching selection differs "<<diff<<" GeV from current energy!"<<endl;
317 
318  // *** initialize triggers
319  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
320  {
321  PndSoftTriggerLine *tl = it->second;
322 
323  if (tl->GetThreshold()<=fEcm)
324  {
325  tl->Init();
326  if (tl->GetRhoTuple()) tl->GetRhoTuple()->GetInternalTree()->SetDirectory(gDirectory);
327 
328  if (fVerbose>0)
329  {
330  cout <<"*** MODE: "<<it->first<<endl;
331  tl->Print();
332  cout <<endl<<endl;
333  }
334  }
335  }
336 
337  // *** Resore original gDirectory
338  dir->cd();
339 
340  if (fVerbose>0)
341  {
342  // print selection setup
343  cout <<"[PndSoftTriggerTask] **** Selection setup:"<<endl;
344  for (std::map<int,STCutSet>::iterator it=fSTSelmap.begin(); it!=fSTSelmap.end(); ++it)
345  {
346  std::cout << it->first << " => ";
347  STCutSet cs = it->second;
348  for (int i=0;i<cs.ncut;++i) cout <<"v["<<cs.varid[i]<<"]"<<fSTOps[cs.op[i]]<<cs.cutval[i]<<" ";
349  cout <<endl;
350  if (cs.tmvaread)
351  {
352  cout<<" TMVA cut:"<<cs.tmvacut<<" vars:";
353  for (int i=0;i<cs.tmvanvar;++i) cout <<cs.tmvavarid[i]<<" ";
354  cout<< endl;
355  }
356  }
357  }
358 
359  return kSUCCESS;
360 }
361 
362 // ----Defaul parameters for QA--------------------------------------------------------------
364 {
365  double Pi0Mass = DbMass("pi0");
366  double EtaMass = DbMass("eta");
367  double Ks0Mass = DbMass("K_S0");
368 
369  // default windows are roughly +-15*sigma, and 15*sigma_max, if different channels are reconstructed
370  // should not depend explicitly from set mean and sigma values, therefore fixed values are chosen
371 
372  fPi0QaMin = Pi0Mass - 0.06;
373  fPi0QaMax = Pi0Mass + 0.06;
374 
375  fEtaQaMin = EtaMass - 0.10;
376  fEtaQaMax = EtaMass + 0.10;
377 
378  fKs0QaMin = Ks0Mass - 0.10;
379  fKs0QaMax = Ks0Mass + 0.10;
380 
381 }
382 
383 // ----Defaul signal parameters --------------------------------------------------------------
385 {
386  // *** fitted peak values from previous studies
387 
388  fPi0Mean = 0.136; // mean value for pi0 signal
389  fPi0Sigma = 0.0045; // sigma value for pi0 signal
390 
391  fEtaMean = 0.552; // mean value for eta(gg) signal
392  fEtaSigma = 0.009; // sigma value for eta(gg) signal
393 
394  fKs0Mean = 0.497; // mean value for Ks signal
395  fKs0Sigma = 0.008; // sigma value for Ks signal
396 
397 }
398 
399 // ----Set all PID algos at once --------------------------------------------------------------
401 {
402  SetPidAlgoElectron(algo);
403  SetPidAlgoMuon(algo);
404  SetPidAlgoPion(algo);
405  SetPidAlgoKaon(algo);
406  SetPidAlgoProton(algo);
407 }
408 
409 
410 // ----Method to enable/disable QA for single mode --------------------------------------------------
412 {
413  int divi = 1;
414  if (mode<10) divi=100;
415  else if (mode<100) divi=10;
416 
417  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
418  if (it->first/divi == mode) it->second->SetWriteQA(qa);
419 
420 // if (fSTTriggers.find(mode) == fSTTriggers.end()) return;
421 // fSTTriggers[mode]->SetWriteQA(qa);
422 }
423 
424 
425 // ----Method to enable/disable full QA--------------------------------------------------------------
427 {
428  SetQAPi0(qa);
429  SetQAEta(qa);
430  SetQAKs0(qa);
431  SetQAEvent(qa);
432  SetQAMc(qa);
433 
434  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
435  SetQAMode(it->first, qa);
436 }
437 
438 // ----Method to enable/disable tagging for single mode --------------------------------------------------
440 {
441  int divi = 1;
442  if (mode<10) divi=100;
443  else if (mode<100) divi=10;
444 
445  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
446  if (it->first/divi == mode) it->second->SetTagActive(tag);
447 
448 // if (fSTTriggers.find(mode) == fSTTriggers.end()) return;
449 // fSTTriggers[mode]->SetTagActive(tag);
450 }
451 
452 // ----Method to enable/disable full Tagging--------------------------------------------------------------
454 {
455  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
456  it->second->SetTagActive(tag);
457 }
458 
459 
460 // ----Method to set mass selection n_sigmas for 'mode' --------------------------------------------------
462 {
463  if (fSTTriggers.find(mode) == fSTTriggers.end()) return;
464 
465  fSTTriggers[mode]->SetTagNSig(nsig);
466 }
467 
468 
469 // ----Method to set mass selection n_sigmas for all modes -----------------------------------------------
471 {
472  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
473  it->second->SetTagNSig(nsig);
474 }
475 
476 
477 // ----- Method to select true PID candidates
479 {
480  int removed = 0;
481 
482  for (int ii=l.GetLength()-1;ii>=0;--ii)
483  {
484  if ( !(fAnalysis->McTruthMatch(l[ii])) )
485  {
486  l.Remove(l[ii]);
487  removed++;
488  }
489  }
490 
491  return removed;
492 }
493 // -------------------------------------------------------------------------
494 
495 // ----- Method to select candidate with certain PID probability
497 {
498  int removed = 0;
499 
500  if (pididx>=0 && pididx<5)
501  for (int ii=l.GetLength()-1;ii>=0;--ii)
502  if (l[ii]->GetPidInfo(pididx)<cut )
503  {
504  l.Remove(l[ii]);
505  removed++;
506  }
507 
508  return removed;
509 }
510 // -------------------------------------------------------------------------
511 
512 // ----- Method to count candidates with certain PID probability
513 int PndSoftTriggerTask::MultPidProb(RhoCandList &l, int pididx, double prob)
514 {
515  int mult = 0;
516 
517  if (pididx>=0 && pididx<5)
518  for (int ii=0;ii<l.GetLength();++ii)
519  if (l[ii]->GetPidInfo(pididx)>=prob) ++mult;
520 
521  return mult;
522 }
523 // -------------------------------------------------------------------------
524 
525 
526 
527 // -------------------------------------------------------------------------
528 // reads in trigger line definitions
529 // defaults are in $VMCWORKDIR/softrig/triggerlines.cfg
531 {
532  std::ifstream file(fTriggerFileName.Data());
533  if (!file.is_open())
534  {
535  cout <<"[PndSoftTriggerTask] **** Unable to open trigger file "<<fTriggerFileName.Data()<<endl;
536  fApplyFullSelection=false;
537  }
538  else cout <<"[PndSoftTriggerTask] **** Reading trigger lines from "<<fTriggerFileName.Data()<<endl;
539 
540  char line[500];
541  TString toks[13];
542 
543  // loop through file line by line
544  while (!file.eof())
545  {
546  file.getline(line,499);
547  TString sline(line);
548 
549  // empty lines or comments are skipped (a comment has to have a # as _first_ char)
550  if (sline.BeginsWith("#")||sline=="") continue;
551 
552  // split the line into tokens, has to be
553  // 0-mode, 1-name, 2-decay pattern, 3-prefix, 4-ntpname, 5-qamin, 6-qamax, 7-mean, 8-sigma, 9-threshold, 10-nsig, 11-fqa, 12-ftag
554  int N = SplitString(sline, ":", toks,13);
555  if (N!=13) {cout <<"[PndSoftTriggerTask] **** Invalid trigger line: "<<sline.Data()<<endl; continue;}
556  //if (fVerbose>1) for (int i=0;i<N;++i) cout <<toks[i]<<":"; cout <<endl;
557 
558  int mode = toks[0].Atoi();
559  PndSoftTriggerLine *tl = new PndSoftTriggerLine(mode, toks[1], toks[2], toks[3], toks[4]);
560 
561  double mean=0, qamin = toks[5].Atof(), qamax = toks[6].Atof();
562 
563  // check whether mean = Ecm
564  if (toks[7]=="Ecm")
565  mean = fEcm;
566  else
567  mean = toks[7].Atof();
568 
569  // check whether qa mass window values are relative to mean
570  if (toks[5].BeginsWith("-")) qamin = mean + toks[5].Atof(); // sum, since toks[5].Atof()<0 due to '-' sign
571  if (toks[6].BeginsWith("+")) qamax = mean + toks[6].Atof();
572 
573  tl->SetQAMassWindow(qamin, qamax);
574  tl->SetMeanSigma(mean, toks[8].Atof());
575  tl->SetThreshold(toks[9].Atof());
576  tl->SetTagNSig(toks[10].Atof());
577  tl->SetWriteQA(toks[11].Atoi());
578  tl->SetTagActive(toks[12].Atoi());
579 
580  fSTTriggers[mode] = tl;
581  }
582 
583  cout << "[PndSoftTriggerTask] **** Found "<<fSTTriggers.size()<<" trigger definitions: ";
584  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
585  cout << it->second->GetName()<<" ";
586  cout <<endl;
587 
588  return true;
589 }
590 
591 
592 // -------------------------------------------------------------------------
593 // reads in specific trigger selections for certain modes and energies
594 // example selection (from ST report - full sim) can be found in the 2 files
595 // $VMCWORKDIR/softrig/selection_10ch_tight.cfg (suppression of 1/1000)
596 // $VMCWORKDIR/softrig/selection_10ch_loose.cfg (eff of 95%)
598 {
599  for (int i=0;i<fSTNEvVars;++i) fSTVarmap[fSTenames[i]] = i;
600 
601  std::ifstream file(fCfgFileName.Data());
602 
603  if (!file.is_open())
604  {
605  cout <<"[PndSoftTriggerTask] **** Unable to open selection file "<<fCfgFileName.Data()<<endl;
606  fApplyFullSelection=false;
607  }
608  else cout <<"[PndSoftTriggerTask] **** Reading selection from "<<fCfgFileName.Data()<<endl;
609 
610  char line[1500];
611 
612  TString toks[200];
613  TString cuts[200];
614  TString toks2[200];
615 
616  // loop through file line by line
617  while (!file.eof())
618  {
619  file.getline(line,1499);
620  TString sline(line);
621 
622  // remove tabs everywhere
623  sline.ReplaceAll("\t","");
624 
625  // if lines contains a '#' somewhere, cut string from that position (=comment)
626  if (sline.Contains("#")) sline = sline(0,sline.Index("#")-1);
627 
628  // remove whitespace at begin and end
629  sline = sline.Strip(TString::kBoth);
630 
631  if (sline=="") continue;
632 
633  // split the line into tokens; token 0 ist the mode code, token 1 is the complete cut string
634  // e.g. '38400 : eslnpidp>0&&abs(pcm-2.105)<0.695&&esthr>0.9&&pt>0.8'
635  int N1 = SplitString(sline, ":", toks,200);
636 
637  // if N==2, only simple cut; N==3 includes TMVA selector
638  if (N1<2 || N1>3) {cout <<"invalid line: "<<sline.Data()<<endl; continue;}
639 
640  // extract the mode code by converting to integer
641  int mcode = toks[0].Atoi();
642  int en = mcode / 1000; // energy prefix (e.g. 550 for sqs=5.5GeV)
643 
644  // if new energy, store the energy prefix in fSTencode
645  if ( std::find(fSTencode.begin(), fSTencode.end(), en) == fSTencode.end() ) fSTencode.push_back(en);
646 
647  // now split the cut string into single cuts
648  int N2 = SplitString(toks[1],"&&",cuts,200);
649 
650  // replace windows cuts 'abs(<name>-x)<y' by 2 single cuts
651  for (int i=0;i<N2;++i)
652  {
653  cuts[i].ReplaceAll(" ","");
654  if (cuts[i]=="tag") cuts[i]="tag>0";
655 
656  // is cut a window cut?
657  if (cuts[i].BeginsWith("abs"))
658  {
659  int pos1 = cuts[i].First('-'); // position of '-'
660  int pos2 = cuts[i].Last(')'); // position of closing abs bracket
661  int pos3 = cuts[i].Last('<'); // position of operator '<'
662 
663  TString name = TString(cuts[i](4,pos1-4)); // name is from first char after 'abs(' to '-' sign
664  double val = TString(cuts[i](pos1+1,(pos2-pos1)-1)).Atof(); // val to cut on is from '-' to closing bracket ')'
665  double win = TString(cuts[i](pos3+1,1000)).Atof(); // window size is from '<' to end
666  cuts[i]= TString::Format("%s>%f",name.Data(), val-win); // construct the 1st and
667  cuts[N2++] = TString::Format("%s<%f",name.Data(), val+win); // 2nd single cut
668  }
669  }
670 
671  // *** now encode the string cuts into sets of (variable index, operator, cutvalue)
672  // *** and put the cut sets into a map with key being the mode code
673 
674  STCutSet cs;
675  cs.ncut = N2;
676 
677  cs.tmvaread = 0;
678  cs.tmvanvar = 0;
679 
680  bool ok = true; // if a variable is not know switch to false
681 
682  for (int i=0;i<N2;++i)
683  {
684  int op=-1, j=0;
685  while (j<5 && op==-1)
686  {
687  if (cuts[i].Contains(fSTOps[j]))
688  {
689  SplitString(cuts[i],fSTOps[j],toks2,5);
690  op = j;
691  }
692  j++;
693  }
694 
695  if (op<0) // fail
696  {
697  i=N2+1;
698  ok=false;
699  }
700 
701  // event shape variable?
702  if (fSTVarmap.find(toks2[0]) != fSTVarmap.end())
703  {
704  cs.varid[i] = fSTVarmap[toks2[0]];
705  cs.op[i] = op;
706  cs.cutval[i] = toks2[1].Atof();
707  }
708  // candidate variable?
709  else if (CodeVariable(toks2[0])>0)
710  {
711  cs.varid[i] = CodeVariable(toks2[0]);
712  cs.op[i] = op;
713  cs.cutval[i] = toks2[1].Atof();
714 
715  }
716  else // fail
717  {
718  cout <<"[PndSoftTriggerTask] **** Unmapped cut variable: "<<toks2[0].Data()<<". Skipping mode "<<mcode<<"."<<endl;
719  i=N2+1;
720  ok=false;
721  }
722  }
723 
724  // *** do we have a TMVA Selector defined?
725  if (N1==3)
726  {
727  // now split the TMVA string into single tokens
728  int N3 = SplitString(toks[2]," ",toks2,200);
729 
730  if (N3<4) // at least four elements expected; fail
731  {
732  cout <<"[PndSoftTriggerTask] **** Invalid TMVA configuration for mode "<<mcode<<". Skipping."<<endl;
733  ok = false;
734  }
735 
736  if (ok)
737  {
738  // first token is the weightfile name without leading path and ending '.weights.xml'
739  TString wfile = TString(getenv("VMCWORKDIR"))+"/softrig/weights/"+toks2[0]+".weights.xml";
740  cs.tmvameth = toks2[0](toks2[0].Length()-3,3);
741 
742  // last token is the cut to the TMVA output
743  cs.tmvacut = toks2[N3-1].Atof();
744 
745  // tokens inbetween are the variables
746  cs.tmvanvar = N3 - 2;
747 
748  // create the reader
749  cs.tmvaread = new TMVA::Reader("Silent");
750 
751  for (int i=1; i<N3-1; ++i)
752  {
753  if (fSTVarmap.find(toks2[i]) != fSTVarmap.end()) cs.tmvavarid[i-1] = fSTVarmap[toks2[i]];
754  else if (CodeVariable(toks2[i])>0) cs.tmvavarid[i-1] = CodeVariable(toks2[i]);
755  else // fail
756  {
757  cout <<"[PndSoftTriggerTask] **** Unmapped TMVA variable: "<<toks2[i].Data()<<". Skipping mode "<<mcode<<"."<<endl;
758  i=N2+1;
759  ok=false;
760  }
761  if (ok) cs.tmvaread->AddVariable(toks2[i], &fSTInputValues[i-1]);
762  }
763  if (ok) cs.tmvaread->BookMVA(cs.tmvameth,wfile);
764  }
765  }
766 
767  if (ok) fSTSelmap[mcode] = cs;
768 
769  if (fVerbose)
770  {
771  cout <<"mcode="<<mcode<<endl;
772  STCutSet cts = fSTSelmap[mcode];
773  cout <<" Ncut="<<cts.ncut<<" (";
774  for (int k=0;k<cts.ncut;++k) cout <<cts.varid[k]<<" ";
775  cout <<")"<<endl;
776  if (cts.tmvaread!=0)
777  {
778  cout <<" TMVAMeth="<<cts.tmvameth<<"; Nvar="<<cts.tmvanvar<<"; Cut="<<cts.tmvacut<<"; TMVAobject="<<cts.tmvaread<<" (";
779  for (int k=0;k<cts.tmvanvar;++k) cout <<cts.tmvavarid[k]<<" ";
780  cout <<")"<<endl;
781  }
782  }
783  }
784  return true;
785 }
786 
787 // -------------------------------------------------------------------------
788 
789 int PndSoftTriggerTask::SplitString(TString s, TString delim, TString *toks, int maxtoks)
790 {
791  TObjArray *tok = s.Tokenize(delim);
792  int N = tok->GetEntries();
793  for (int i=0;i<N;++i)
794  if (i<maxtoks)
795  {
796  toks[i] = ((TObjString*)tok->At(i))->String();
797  toks[i].ReplaceAll("\t","");
798  toks[i] = toks[i].Strip(TString::kBoth);
799  }
800  return N;
801 }
802 
803 // -------------------------------------------------------------------------
804 
806 {
807  // Get run and runtime database
808  FairRun* run = FairRun::Instance();
809  if ( ! run ) Fatal("SetParContainers", "No analysis run");
810 }
811 
812 // -------------------------------------------------------------------------
813 
814 
815 // ----- Public method Exec --------------------------------------------
817 {
818  // *** prepare TCA for writing OnlineFilterInfo
819  if ( fTcaOnlineFilterInfo->GetEntriesFast() != 0 ) fTcaOnlineFilterInfo->Delete();
820 
821  // *** some variables
822  //int i=0,j=0, k=0; //[R.K. 01/2017] unused variable
823 
824  if (!(++fEvtCount%100)) cout << "[PndSoftTriggerTask] evt "<<fEvtCount<<endl;
825 
826  // *** read the next event
828 
829  // *** fill all lists necessary for combinatorics
830  FillGlobalLists();
831 
832  // *** go through MC truth list and codify the recoil mode (pi+-, pi0, eta, gamma, K+-)
834 
835  // *** estimate primary vertex
837 
838  // *** setup eventshape object
840  fEventShape = &evtShape;
841 
843 
844  int tag_glob = 0, tag_glob_pre = 0, tag_pre = 0;
845 
846  PndOnlineFilterInfo* info=new ( (*fTcaOnlineFilterInfo)[0] ) PndOnlineFilterInfo();
847 
848  // *** loop through all channels and tag
849  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
850  {
851  PndSoftTriggerLine *tl = it->second;
852 
853  if (tl->GetThreshold()<=fEcm && tl->GetTagActive() )
854  {
855  int ntag = TagMode( tl , tag_pre);
856  tl->SetNTagged( ntag );
857  info->SetNTag(it->first, ntag);
858  tag_glob += ntag;
859  tag_glob_pre += tag_pre;
860  }
861  }
862 
863  Int_t tagged = (tag_glob>0);
864  Int_t taggedm = (tag_glob_pre>0); // tagged by mass cut only
865 
866  // *** write common information
867  if (fQAEvent)
868  {
869  ntp->Column("ev", (Int_t) fEvtCount, 0);
870  ntp->Column("run", (Int_t) fRunNum, 0);
871  ntp->Column("mode", (Int_t) fMode, 0);
872  ntp->Column("recmode", (Int_t) fRecoilMode, 0);
873  ntp->Column("reccnt", (Int_t) fRecoilCnt, 0);
874 
875  fQA->qaP4("beam", fIniP4, ntp);
876  ntp->Column("primvx", (Float_t) fPrimVtx.X(), 0.0f);
877  ntp->Column("primvy", (Float_t) fPrimVtx.Y(), 0.0f);
878  ntp->Column("primvz", (Float_t) fPrimVtx.Z(), 0.0f);
879  ntp->Column("primvqa", (Float_t) fPrimVtxQa , 0.0f);
880 
881  // write tags of individual modes
882  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
883  {
884  PndSoftTriggerLine *tl = it->second;
885  if ( tl->GetThreshold()<=fEcm && tl->GetTagActive() )
886  {
887  //ntp->Column("tag_"+tl->GetName(), (Int_t) tl->GetNTagged(), 0);
888  ntp->Column(TString::Format("tag%d",tl->GetModeCode()), (Int_t) tl->GetNTagged(), 0);
889  }
890  }
891 
892  ntp->Column("tagall", (Int_t) tag_glob, 0);
893  ntp->Column("tagpre", (Int_t) tag_glob_pre, 0);
894  ntp->Column("tag", (Int_t) tagged, 0);
895  ntp->Column("tagm", (Int_t) taggedm, 0);
896 
897  fQA->qaEventShape("es", fEventShape, ntp);
898 
899  // replace PID mult values from fEventShape by actual counts with individual algorithms
900  ntp->Column("eslnpide", (Int_t) fPidMult_025[0], 0 );
901  ntp->Column("eslnpidmu",(Int_t) fPidMult_025[1], 0 );
902  ntp->Column("eslnpidpi",(Int_t) fPidMult_025[2], 0 );
903  ntp->Column("eslnpidk", (Int_t) fPidMult_025[3], 0 );
904  ntp->Column("eslnpidp", (Int_t) fPidMult_025[4], 0 );
905 
906  ntp->DumpData();
907  }
908 
909  // *** write out MC info
910  if (fQAMc) {fQA->qaMcList("",fMcTruth,nmc); nmc->DumpData();}
911 
912 }
913 
914 
916 {
917  if (ntp) ntp->GetInternalTree()->Write(); // overall info
918  if (nks0) nks0->GetInternalTree()->Write(); // Ks0 QA
919  if (npi0) npi0->GetInternalTree()->Write(); // pi0 QA
920  if (neta) neta->GetInternalTree()->Write(); // eta QA
921  if (nmc) nmc->GetInternalTree()->Write(); // MC QA
922 
923  for (TrigIt it=fSTTriggers.begin(); it!=fSTTriggers.end(); ++it)
924  {
925  PndSoftTriggerLine *tl = it->second;
926  if (tl->GetRhoTuple())
927  tl->GetRhoTuple()->GetInternalTree()->Write();
928  }
929 }
930 
931 // -------------------------------------------------------------------------
932 // determine recoil mode (N_gamma + 10*N_pi0 + 100*N_pi+- + 1000*N_K+- + 10000*N_K0 + 100000*N_eta)
934 {
935  mode = 0;
936  int code=-1;
937 
938  if ( fMcTruth.GetLength()==0 ) return -1;
939  if ( fMcTruth[0]==0 ) return -1;
940 
941  for (int i=0; i<fMcTruth[0]->NDaughters(); ++i)
942  {
943  int dpdg = abs(fMcTruth[0]->Daughter(i)->PdgCode());
944 
945  switch (dpdg)
946  {
947  case 22: mode += 1; break; // gamma
948  case 111: mode += 10; break; // pi0
949  case 211: mode += 100; break; // pi+
950  case 321: mode += 1000; break; // K+
951  case 130: mode += 10000; break; // K0
952  case 310: mode += 10000; break; // K0
953  case 221: mode += 100000; break; // eta
954  }
955  }
956 
957  switch (mode)
958  {
959  case 0 : code = 0; break; // no recoil
960  case 1 : code = 1; break; // gamma
961  case 10 : code = 2; break; // pi0
962  case 100000 : code = 3; break; // eta
963  case 20 : code = 4; break; // pi0 pi0
964  case 200 : code = 5; break; // pi+ pi-
965  case 2000 : code = 6; break; // K+ K-
966  case 20000 : code = 7; break; // K0 K0b
967  case 200000 : code = 8; break; // eta eta
968  case 210 : code = 9; break; // pi+ pi- pi0
969  default : code =-1; break;
970  }
971 
972  return code;
973 }
974 
975 // -------------------------------------------------------------------------
976 
978 {
979  // *** basic KS0 Reco
980  //fKs0Cands.Combine(fPionPlus, fPionMinus);
982  // *** pre select for qa
984  fKs0Cands.SetType(310);
985 
986  if (n)
987  {
988  for (int i=0; i<fKs0Cands.GetLength();++i)
989  {
990  Float_t nsig = (Float_t) fabs(fKs0Cands[i]->Mass()-fKs0Mean)/fKs0Sigma;
991  Float_t tag = nsig<fNsigAux; // instead of fNsigTag, fNsigAux cut is applied for KS
992 
993  n->Column("tag", (Int_t) tag, 0);
994  n->Column("nsig", (Float_t) nsig, 0.0f);
995  n->Column("ev", (Int_t) fEvtCount, 0);
996  n->Column("run", (Int_t) fRunNum, 0);
997  n->Column("num", (Int_t) i, 0);
998  n->Column("mode", (Int_t) fMode, 0);
999  n->Column("xmean", (Float_t) fKs0Mean, 0.0f);
1000  n->Column("xsig", (Float_t) fKs0Sigma, 0.0f);
1001 
1002  fQA->qaP4("beam", fIniP4, n);
1003  fQA->qaKs0("x", fKs0Cands[i], n);
1004 
1005  n->DumpData();
1006  }
1007  }
1008 
1009  // *** final selection
1011 
1012  return fKs0Cands.GetLength();
1013 }
1014 
1015 // -------------------------------------------------------------------------
1016 
1017 
1019 {
1020  int i;
1021 
1022  // *** fill standard lists
1023  fAnalysis->FillList( fMcTruth, "McTruth" );
1026  fAnalysis->FillList( fNeutralCands, "Neutral" );
1027 
1028  // *** fill PID lists
1029 
1030  fAnalysis->FillList( fPidList[0], "ElectronAllPlus", fAlgoElectron );
1031  fAnalysis->FillList( fPidList[1], "ElectronAllMinus", fAlgoElectron );
1032 
1033  fAnalysis->FillList( fPidList[2], "MuonAllPlus", fAlgoMuon );
1034  fAnalysis->FillList( fPidList[3], "MuonAllMinus", fAlgoMuon );
1035 
1036  fAnalysis->FillList( fPidList[4], "PionAllPlus", fAlgoPion );
1037  fAnalysis->FillList( fPidList[5], "PionAllMinus", fAlgoPion );
1038 
1039  fAnalysis->FillList( fPidList[6], "KaonAllPlus", fAlgoKaon );
1040  fAnalysis->FillList( fPidList[7], "KaonAllMinus", fAlgoKaon );
1041 
1042  fAnalysis->FillList( fPidList[8], "ProtonAllPlus", fAlgoProton );
1043  fAnalysis->FillList( fPidList[9], "ProtonAllMinus", fAlgoProton );
1044 
1045  // *** apply minimum PID cut if required
1046  if (fIniPidCut>0) for (i=0;i<10;++i) SelectPidProb(fPidList[i],i/2, fIniPidCut);
1047 
1048  // *** select on lists
1050 
1051  // *** apply minimum momentum cut if required
1052  if (fTrackMinP>0) for (i=0;i<10;++i) fPidList[i].Select(fMomentumSel);
1053 
1054  // *** count PID multiplicities
1055  for (i=0;i<5;++i) fPidMult_025[i] = MultPidProb(fPidList[2*i], i, 0.25) + MultPidProb(fPidList[2*i+1], i, 0.25);
1056 
1057  // *** Ks reco
1059 
1060  // *** pi0 reco
1063  fPi0Cands.SetType(111);
1064  if (fQAPi0)
1065  {
1066  for (i=0; i<fPi0Cands.GetLength();++i)
1067  {
1068  Float_t nsig = (Float_t) fabs(fPi0Cands[i]->Mass()-fPi0Mean)/fPi0Sigma;
1069  Float_t tag = nsig<fNsigAux; // instead of fNsigTag, a fNsigAux cut is applied for pi0
1070 
1071  npi0->Column("tag", (Int_t) tag, 0);
1072  npi0->Column("nsig", (Float_t) nsig, 0.0f);
1073  npi0->Column("ev", (Int_t) fEvtCount, 0);
1074  npi0->Column("run", (Int_t) fRunNum, 0);
1075  npi0->Column("num", (Int_t) i, 0);
1076  npi0->Column("mode", (Int_t) fMode, 0);
1077  npi0->Column("xmean", (Float_t) fPi0Mean, 0.0f);
1078  npi0->Column("xsig", (Float_t) fPi0Sigma, 0.0f);
1079 
1080  fQA->qaP4("beam", fIniP4, npi0);
1081  fQA->qaPi0("x", fPi0Cands[i], npi0);
1082  npi0->DumpData();
1083  }
1084  }
1086 
1087  // *** era reco
1090  fEtaCands.SetType(221);
1091  if (fQAEta)
1092  {
1093  for (i=0; i<fEtaCands.GetLength();++i)
1094  {
1095  Float_t nsig = (Float_t) fabs(fEtaCands[i]->Mass()-fEtaMean)/fEtaSigma;
1096  Float_t tag = nsig<fNsigAux; // instead of fNsigTag, a fNsigAux cut is applied for eta
1097 
1098  neta->Column("tag", (Int_t) tag, 0);
1099  neta->Column("nsig", (Float_t) nsig, 0.0f);
1100  neta->Column("ev", (Int_t) fEvtCount, 0);
1101  neta->Column("run", (Int_t) fRunNum, 0);
1102  neta->Column("num", (Int_t) i, 0);
1103  neta->Column("mode", (Int_t) fMode, 0);
1104  neta->Column("xmean", (Float_t) fEtaMean, 0.0f);
1105  neta->Column("xsig", (Float_t) fEtaSigma, 0.0f);
1106 
1107  fQA->qaP4("beam", fIniP4, neta);
1108  fQA->qaPi0("x", fEtaCands[i], neta);
1109 
1110  neta->DumpData();
1111  }
1112  }
1114 
1115  // *** copy all lists to indexed ones for combinatorics
1116 
1117  fPidList[10] = fGammaCands;
1118 
1119  fPidList[11] = fPi0Cands;
1120  fPidList[12] = fKs0Cands;
1121  fPidList[13] = fEtaCands;
1122 }
1123 
1124 // -------------------------------------------------------------------------
1125 // Fill the array of event shape variables (called once per event)
1126 // here the connection between the variable indices and the values is made
1127 
1129 {
1130  // *** Event Shape Vars
1131 
1132  // ( 0) eslnpide ( 1) eslnpidmu ( 2) eslnpidpi ( 3) eslnpidk ( 4) eslnpidp
1133  // ( 5) esfw1 ( 6) esfw2 ( 7) esfw3 ( 8) esfw4 ( 9) esfw5
1134  // ( 10) esapl ( 11) escir ( 12) espla ( 13) esnneut ( 14) esnpart
1135  // ( 15) esnchrg ( 16) espmax ( 17) espmaxl ( 18) espmin ( 19) espminl
1136  // ( 20) esprapmax ( 21) esptmax ( 22) essumen ( 23) essumenl ( 24) essumetn
1137  // ( 25) essumpc ( 26) essumpc05 ( 27) essumpcl ( 28) essumpt ( 29) essumptc
1138  // ( 30) esthr ( 31) essph ( 32) esptmin ( 33) esncp10l ( 34) esnne10l
1139 
1140  // don't use PID mult values from fEventShape (based on AllCands and only one algo)
1141  for (int i=0;i<5;++i) fSTVarEvArray[i] = fPidMult_025[i];
1142 
1148 
1154 
1160 
1162  fSTVarEvArray[21] = fEventShape->Ptmax();
1166 
1172 
1173  fSTVarEvArray[30] = fEventShape->Thrust();
1175  fSTVarEvArray[32] = fEventShape->Ptmin();
1178 
1179 }
1180 
1181 // -------------------------------------------------------------------------
1182 
1183 TLorentzVector PndSoftTriggerTask::BoostCms(TLorentzVector l)
1184 {
1185  TLorentzVector lcm = l;
1186  lcm.Boost(-fIniP4.BoostVector());
1187  return lcm;
1188 }
1189 
1190 // -------------------------------------------------------------------------
1191 void PndSoftTriggerTask::GetAngles(RhoCandidate *c, double &oang, double &decang)
1192 {
1193  oang = -999.;
1194  decang = -999.;
1195 
1196  if (c->NDaughters()!=2) return;
1197 
1198  RhoCandidate *d0 = c->Daughter(0);
1199  RhoCandidate *d1 = c->Daughter(1);
1200 
1201  // opening angle lab
1202  oang = d0->P3().Angle(d1->P3());
1203 
1204  // decay angle
1205  TLorentzVector d_cms = d0->P4();
1206  d_cms.Boost(-(c->P4().BoostVector()));
1207  decang = d_cms.Vect().Angle(c->P3());
1208 }
1209 
1210 
1211 // -------------------------------------------------------------------------
1212 // *** Get Poca vertex info (for secondary vertex cuts)
1213 double PndSoftTriggerTask::GetPocaVtx(RhoCandidate* c, double &dist, double &ctau)
1214 {
1215  // *** simple vtx finder
1216  TVector3 vtx, altvtx, primvtx;
1217  double qavtx = fPocaVertexer->GetPocaVtx(vtx, c);
1218 
1219  // *** determine poca of rest of tracks
1220  RhoCandList l;
1221  fAnalysis->FillList(l, "Charged");
1222 
1223  // remove c's daughters from charged tracks list
1224  l.RemoveFamily(c);
1225 
1226  // if at least 2 tracks left, compute the poca vertex of the residual tracks
1227  if (l.GetLength()>1) fPocaVertexer->GetPocaVtx(altvtx, l);
1228  else altvtx = fPrimVtx; // else use the primary vtx found before
1229 
1230  dist = 999.;
1231  // *** if vertex of either all charged or residual tracks was found, compute distance and ctau
1232  if (altvtx.Mag()>0.) dist = (vtx-altvtx).Mag();
1233 
1234  ctau = dist*c->M()/c->P();
1235 
1236  return qavtx;
1237 }
1238 
1239 
1240 // -------------------------------------------------------------------------
1241 // transform a variable string to an integer code
1242 // code conventions are for code ABCDEF
1243 // - ABCDEF < 100: special variable directly filled (e.g. mmiss = 80, mdif = 81, tag = 99)
1244 // - ABCD defines candidate (e.g. x=9999, xd0=9099, xd1=9199, xd0d1=9019, xd0d0g2=9001 etc)
1245 // - EF defines variable:
1246 // - E: 1=kin, 2=pid, 3=vtx, 4=ang
1247 //
1248 // E\F | 0 | 1 | 2 | 3 | 4 | 5
1249 // ----------------------------------------------------------------------
1250 // 1 (kin) | p | pt | tht | pcm | thtcm |
1251 // 2 (pid) | pide | pidmu | pidpi | pidk | pidp | pidmax
1252 // 3 (vtx) | pocqa | pocdist | pocctau | | |
1253 // 4 (ang) | oang | decang | cdecang | | |
1254 
1256 {
1257  int i=0, code=0, A=9, B=9, C=9, D=9, EF=0;
1258 
1259  if (v=="mmiss") return 80;
1260  if (v=="xmdif") return 82;
1261  if (v=="tag") return 99;
1262 
1263  for (i=0; i<fSTNQuant; ++i)
1264  if (v.EndsWith(fSTVarSuff[i]))
1265  {
1266  EF=fSTQuantCode[i];
1267  v=v(0,v.Length()-fSTVarSuff[i].Length());
1268  i=1000;
1269  }
1270 
1271  if (i<1000) {cout <<"[PndSoftTriggerTask] ***** Unknown variable suffix in variable '"<<v<<"'"<<endl; return -1;}
1272 
1273  // get rid of special daughter namings for pi0 gammas and KS pi+-
1274  v.ReplaceAll("pi1","d0");
1275  v.ReplaceAll("pi2","d1");
1276  v.ReplaceAll("g1","d0");
1277  v.ReplaceAll("g2","d1");
1278 
1279  TRegexp r1("xd[0-4]$");
1280  TRegexp r2("xd[0-4]d[0-4]$");
1281  TRegexp r3("xd[0-4]d[0-4]d[0-4]$");
1282 
1283  if (v!="x" && v(r1)=="" && v(r2)=="" && v(r3)=="") {cout <<"[PndSoftTriggerTask] ***** Unknown candidate prefix in variable '"<<v<<"'"<<endl; return -1;}
1284 
1285  int l = v.Length();
1286 
1287  if (l>1) B = TString(v(2,1)).Atoi();
1288  if (l>3) C = TString(v(4,1)).Atoi();
1289  if (l>5) D = TString(v(6,1)).Atoi();
1290 
1291  code=EF+100*D+1000*C+10000*B+100000*A;
1292 
1293  return code;
1294 }
1295 
1296 // -------------------------------------------------------------------------
1297 // Fill individual candidate specific variables (called for every candidate)
1298 // here the connection between the variable indices and the values is made
1299 // in case variables are connected with more computing effort, they are filled simultaneously
1300 
1302 {
1303  // previously filled event shape var?
1304  if (id>=0 && id<35) return fSTVarEvArray[id];
1305 
1306  // is variable already in cache?
1307  if (fSTVarCandArray.find(id) != fSTVarCandArray.end()) return fSTVarCandArray[id];
1308 
1309  // *** now compute variable values based on code 'id'
1310 
1311  // these are special cases with simple code numbers
1312  if (id==80) return (fIniP4 - c->P4()).M(); // mmiss
1313  if (id==82) return (c->M() - c->Daughter(0)->M()); // D* mdif
1314  if (id==99) return 1.0; // artificial value for tag>0 (= empty cut)
1315 
1316  // determine quantity and candidate from id code
1317  int quant = id%100;
1318  int cand = id/100;
1319  int B = (id/10000)%10, C = (id/1000)%10, D = (id/100)%10;
1320 
1321  // find particle under consideration
1322  RhoCandidate *cnd=0;
1323 
1324  if (B==9) cnd = c;
1325  else if (C==9) cnd = c->Daughter(B);
1326  else if (D==9) cnd = c->Daughter(B)->Daughter(C);
1327  else cnd = c->Daughter(B)->Daughter(C)->Daughter(D);
1328 
1329  // coding for variable suffix (to translate variable names to code)
1330  // E\F | 0 | 1 | 2 | 3 | 4 | 5
1331  // ----------------------------------------------------------------------
1332  // 1 (kin) | p | pt | tht | pcm | thtcm |
1333  // 2 (pid) | pide | pidmu | pidpi | pidk | pidp | pidmax
1334  // 3 (vtx) | pocqa | pocdist | pocctau | | |
1335  // 4 (ang) | oang | decang | cdecang | | |
1336 
1337  //double val = -999.; //[R.K. 01/2017] unused variable
1338  double pocdist, pocctau, oang, decang, cdecang, pidmax; //pocqa, //[R.K. 01/2017] unused variable
1339  TLorentzVector bl;
1340 
1341  switch (quant)
1342  {
1343  case 10 : // p
1344  fSTVarCandArray[id] = cnd->P(); break;
1345 
1346  case 11 : // pt
1347  fSTVarCandArray[id] = cnd->Pt(); break;
1348 
1349  case 12 : // tht
1350  fSTVarCandArray[id] = cnd->P4().Theta(); break;
1351 
1352  case 13 : case 14 : // pcm, thtcm
1353  bl = BoostCms(cnd->P4());
1354  fSTVarCandArray[cand*100+13] = bl.P();
1355  fSTVarCandArray[cand*100+14] = bl.Theta();
1356  break;
1357 
1358  case 20 : case 21 : case 22 : case 23 : case 24 : // pide,..., pidp
1359  fSTVarCandArray[id] = cnd->GetPidInfo(quant%10);
1360  break;
1361 
1362  case 25 : // pidmax
1363  pidmax = 0;
1364  for (int i=0;i<5;++i) if (cnd->GetPidInfo(i)>pidmax) pidmax = cnd->GetPidInfo(i);
1365  fSTVarCandArray[id] = pidmax;
1366  break;
1367 
1368  case 30 : case 31 : case 32 : // pocqa, pocdist, pocctau
1369  fSTVarCandArray[cand*100+30] = GetPocaVtx(cnd, pocdist, pocctau);
1370  fSTVarCandArray[cand*100+31] = pocdist;
1371  fSTVarCandArray[cand*100+32] = pocctau;
1372  break;
1373 
1374  case 40 : case 41 : case 42 : // oang, decang, cdecang
1375  GetAngles(cnd, oang, decang);
1376  cdecang = cos(decang);
1377  fSTVarCandArray[cand*100+40] = oang;
1378  fSTVarCandArray[cand*100+41] = decang;
1379  fSTVarCandArray[cand*100+42] = cdecang;
1380  break;
1381  default:
1382  cout <<"[PndSoftTriggerTask] ***** wrong variable code '"<<quant<<"'"<<endl;
1383  }
1384 
1385  return fSTVarCandArray[id];
1386 }
1387 
1388 // -------------------------------------------------------------------------
1389 // Fill the array of candidate specific variables (called for every candidate)
1390 // here the connection between the variable indices and the values is made
1391 
1393 {
1394  int i=0;
1395 
1396  fSTVarCandArray.clear(); // reset cand var cache (used during fill of values)
1397 
1398  // fill all requested variables
1399  if (tmva)
1400  for (i=0; i<fSTSelmap[mcode].tmvanvar; ++i)
1401  fSTInputValues[i] = GetVarValue(c, fSTSelmap[mcode].tmvavarid[i]); // fill variables for tmva
1402  else
1403  for (i=0; i<fSTSelmap[mcode].ncut; ++i)
1404  fSTInputValues[i] = GetVarValue(c, fSTSelmap[mcode].varid[i]); // fill variables for cut selection
1405 }
1406 
1407 // -------------------------------------------------------------------------
1408 // *** In case a particle is a D*0, D*+, D_s*+, apply cut on mass difference
1409 // *** for QA NTuple writeout (smaller file size)
1411 {
1412  int pdg = abs(c->PdgCode());
1413 
1414  if ( (pdg==413 || pdg==423 || pdg==433) && fabs(c->M() - c->Daughter(0)->M() - 0.143) > fDstMDiffCut ) return false;
1415 
1416  return true;
1417 }
1418 
1419 // -------------------------------------------------------------------------
1420 // *** Apply full selection to a candidate
1422 {
1423  int mcode = fSTencode[fSTModeIndex]*1000+mode;
1424 
1425  if (fVerbose>1) cout <<"AcceptCandidate : mode="<<mode<<" mcode="<<mcode;
1426 
1427  // no selection defined for this mode
1428  if ( fSTSelmap.find(mcode) == fSTSelmap.end() ) {if (fVerbose>1) cout <<endl;return false;}
1429  if (fVerbose>1) cout <<" found cut set";
1430 
1431  // not accepted by final (mass) selector
1432  if ( sel && !(sel->Accept(c)) ) {if (fVerbose>1) cout <<endl;return false;}
1433  // not accepted by D* mass diff cut
1434  if (!AcceptDstCut(c)) return false;
1435 
1436  if (fVerbose>1) cout <<" accepted by precuts"<<endl;
1437 
1438  // ***
1439  // *** check for cut selection
1440  // ***
1441  STCutSet cs = fSTSelmap[mcode];
1442  FillVarArray(c, mcode);
1443 
1444  bool acc = true;
1445 
1446  for (int i=0;i<cs.ncut;++i)
1447  {
1448  if (fVerbose>1) cout <<" -> checking var["<<cs.varid[i]<<"] ("<<fSTInputValues[i]<<") "<<fSTOps[cs.op[i]]<<" "<<cs.cutval[i]<<endl;
1449 
1450  switch (cs.op[i]) // which operator is used for this cut?
1451  {
1452  case 0: // check var >= value
1453  if ( !(fSTInputValues[i]>=cs.cutval[i]) ) acc=false;
1454  break;
1455  case 1: // check var > value
1456  if ( !(fSTInputValues[i]>cs.cutval[i]) ) acc=false;
1457  break;
1458  case 2: // check var == value
1459  if ( !(fSTInputValues[i]==cs.cutval[i]) ) acc=false;
1460  break;
1461  case 3: // check var <= value
1462  if ( !(fSTInputValues[i]<=cs.cutval[i]) ) acc=false;
1463  break;
1464  case 4: // check var < value
1465  if ( !(fSTInputValues[i]<cs.cutval[i]) ) acc=false;
1466  break;
1467  default :
1468  acc=false;
1469  break;
1470  }
1471  }
1472 
1473  // ***
1474  // *** check for TMVA selector in addition
1475  // ***
1476  if (cs.tmvaread)
1477  {
1478  FillVarArray(c, mcode, true);
1479  if (fVerbose>1) cout <<" -> checking TMVA "<<cs.tmvameth<<endl;
1480  float mvaout = cs.tmvaread->EvaluateMVA(cs.tmvameth);
1481  if (mvaout<cs.tmvacut) acc=false;
1482  }
1483 
1484 
1485  if (fVerbose>1)
1486  {
1487  if (acc) cout <<" final accept"<<endl;
1488  else cout <<endl;
1489  }
1490 
1491  return acc;
1492 }
1493 
1494 
1495 // -------------------------------------------------------------------------
1496 // *** Determine pdg code of anti particle
1498 {
1499  int antipdg = pdg;
1500  TParticlePDG *part = fPdg->GetParticle(pdg)->AntiParticle();
1501  if (part) antipdg = part->PdgCode();
1502 
1503  return antipdg;
1504 }
1505 
1506 // -------------------------------------------------------------------------
1507 // *** create a list based on indices
1508 void PndSoftTriggerTask::CombineList(RhoCandList &l, int mothpdg, int amothpdg, std::vector<int> &idx, std::vector<int> &aidx, bool cc)
1509 {
1510  l.Cleanup();
1511 
1512  int nd = idx.size();
1513 
1514  switch (nd)
1515  {
1516  case 2:
1517  l.Combine(fPidList[idx[0]], fPidList[idx[1]], mothpdg);
1518  if (cc) l.CombineAndAppend(fPidList[aidx[0]], fPidList[aidx[1]], amothpdg);
1519  break;
1520 
1521  case 3:
1522  l.Combine(fPidList[idx[0]], fPidList[idx[1]], fPidList[idx[2]], mothpdg);
1523  if (cc) l.CombineAndAppend(fPidList[aidx[0]], fPidList[aidx[1]], fPidList[aidx[2]], amothpdg);
1524  break;
1525 
1526  case 4:
1527  l.Combine(fPidList[idx[0]], fPidList[idx[1]], fPidList[idx[2]], fPidList[idx[3]], mothpdg);
1528  if (cc) l.CombineAndAppend(fPidList[aidx[0]], fPidList[aidx[1]], fPidList[aidx[2]], fPidList[aidx[3]], amothpdg);
1529  break;
1530 
1531  case 5:
1532  l.Combine(fPidList[idx[0]], fPidList[idx[1]], fPidList[idx[2]], fPidList[idx[3]], fPidList[idx[4]], mothpdg);
1533  if (cc) l.CombineAndAppend(fPidList[aidx[0]], fPidList[aidx[1]], fPidList[aidx[2]], fPidList[aidx[3]], fPidList[aidx[4]], amothpdg);
1534  break;
1535  default: return;
1536  }
1537 
1538 }
1539 
1540 // -------------------------------------------------------------------------
1541 // *** General common combinatorics based on PndSoftTriggerLine input
1543 {
1544  l.Cleanup();
1545  fPidList[14].Cleanup();
1546  fPidList[15].Cleanup();
1547 
1548  int mothpdg = tl->GetMotherPdg();
1549  if (mothpdg<0) {cout << "[PndSoftTriggerTask] **** Invalid mother for combinatorics."<<endl; return 0;}
1550 
1551  int nd = tl->GetNDaughters(); // cache for mapping: pdgcode -> fPdgList indices
1552 
1553  std::vector<int> idx, aidx, auxidx, auxaidx;
1554 
1555 
1556  // fetch pdg code of antiparticle for mother
1557  int amothpdg = AntiPdg(mothpdg);
1558 
1559  // do we need to add charged conjugate combinatorics?
1560  bool cc = tl->GetCC();
1561 
1562  if (fVerbose>1)
1563  {
1564  for (int i=0; i<nd; ++i) cout <<tl->GetDaughterPdg(i)<<" ";
1565  if (cc) cout <<" +cc";
1566  cout <<endl;
1567  }
1568 
1569  // flag whether pdg goes to idx list or auxidx list
1570  bool auxmode=false;
1571 
1572  // idx for the mothers of the aux list
1573  int auxmothpdg, auxamothpdg;
1574 
1575  // cache the indices of the fPidList according to the pdg codes
1576  for (int i=0; i<nd; ++i)
1577  {
1578  int dpdg = tl->GetDaughterPdg(i);
1579 
1580  // skip the codes for '[' and ']'
1581  if (dpdg==-99) continue;
1582  if (dpdg==-98) {auxmode=false; continue;} // filling auf aux list finished
1583 
1584  // if no list for particle with code 'dpdg'
1585  if ( fSTMapListIndex.find(dpdg) == fSTMapListIndex.end() )
1586  {
1587  // is it start of aux list definition? (e.g. D*0 -> D0 [K- pi+] pi0)
1588  if (tl->GetDaughterPdg(i+1) != -99)
1589  {
1590  // if not, throw error!
1591  cout << "[PndSoftTriggerTask] **** Invalid daughter pdg code: "<<dpdg<<endl;
1592  return 0;
1593  }
1594  else // start aux mode and decide about cc of aux list
1595  {
1596  auxmode = true;
1597  auxmothpdg = dpdg;
1598  auxamothpdg = AntiPdg(dpdg);
1599  idx.push_back(14);
1600  if (auxamothpdg==auxmothpdg) aidx.push_back(14);
1601  else aidx.push_back(15);
1602  continue;
1603  }
1604  }
1605  if (auxmode)
1606  {
1607  auxidx.push_back(fSTMapListIndex[dpdg]);
1608  auxaidx.push_back(fSTMapListIndex[AntiPdg(dpdg)]);
1609  }
1610  else
1611  {
1612  idx.push_back(fSTMapListIndex[dpdg]);
1613  aidx.push_back(fSTMapListIndex[AntiPdg(dpdg)]);
1614  }
1615  }
1616 
1617  // does one of the lists have too many daughters?
1618  if (idx.size()>5) {cout << "[PndSoftTriggerTask] **** Too many daughters ("<<idx.size()<<") for combinatorics of target resonance."<<endl;return 0;}
1619  if (auxidx.size()>5) {cout << "[PndSoftTriggerTask] **** Too many daughters ("<<auxidx.size()<<") for combinatorics of auxiliary resonance."<<endl;return 0;}
1620 
1621  if (fVerbose>1)
1622  {
1623  cout <<"("<<mothpdg<<"/"<<amothpdg<<") : ";
1624  for (size_t i=0; i<idx.size(); ++i) cout <<"("<<idx[i]<<"/"<<aidx[i]<<") ";
1625  cout <<endl<<endl;
1626  if (auxidx.size()>0)
1627  {
1628  cout <<"("<<auxmothpdg<<"/"<<auxamothpdg<<") : ";
1629  for (size_t i=0; i<auxidx.size(); ++i) cout <<"("<<auxidx[i]<<"/"<<auxaidx[i]<<") "; cout <<endl;
1630  }
1631  }
1632 
1633  // do the combinatorics for aux list and final list
1634  if (tl->GetAuxNeeded())
1635  {
1636  // create different aux lists is pdg!=anti-pdg (e.g. D0 -> K- pi+ and anti-D0 -> K+ pi-)
1637  if (auxmothpdg!=auxamothpdg)
1638  {
1639  CombineList(fPidList[14], auxmothpdg, 0, auxidx, auxidx, false);
1640  CombineList(fPidList[15], auxamothpdg, 0, auxaidx, auxaidx, false);
1641  }
1642  // create one aux list e.g. for eta_c -> ...
1643  else
1644  {
1645  // check whether final state is its anti-final state (e.g. pi+ pi- pi0 = pi- pi+ pi0)
1646  // or not (e.g. eta_c -> KS K- pi+ != KS K+ pi-)
1647  // this can be done by summing all pdg codes for particles with pdg!=anti-pdg
1648  // if sum = 0, final state and anti-final state are the same (i.e. FS has for each particle the according anti-particle)
1649  int pdgsum = 0;
1650  for (size_t i=0;i<auxidx.size();++i) if (auxidx[i]<10) pdgsum+=fSTPidIndex[auxidx[i]];
1651 
1652  // if FS = anti-FS, just do combinatorics once (not adding the composites from anti-FS list)
1653  if (pdgsum==0)
1654  CombineList(fPidList[14], auxmothpdg, 0, auxidx, auxidx, false);
1655  // if FS != anti-FS, append also comb from anti-FS list (e.g. etac from (KS K- pi+) + (KS K+ pi-))
1656  else
1657  CombineList(fPidList[14], auxmothpdg, auxamothpdg, auxidx, auxaidx, true);
1658  }
1659  }
1660 
1661  // create the final list using (or not using) the aux list
1662  CombineList(l, mothpdg, amothpdg, idx, aidx, cc);
1663 
1664  return l.GetLength();
1665 }
1666 
1667 
1668 // -------------------------------------------------------------------------
1669 // *** General common tagging methode
1671 {
1672  // *** counter for full selection accepted cands
1673  int nacc = 0;
1674  npre = 0;
1675 
1676  if ( !tl->GetTagActive() ) return 0;
1677  if ( fEcm < tl->GetThreshold() ) return 0;
1678 
1679  int mode = tl->GetModeCode();
1680  RhoTuple *n = tl->GetRhoTuple();
1681  TString prefix = tl->GetPrefix();
1682  double mean = tl->GetMean();
1683  double sigma = tl->GetSigma();
1684 
1685  RhoMassParticleSelector *sel = tl->GetSelector();
1686 
1687  RhoCandList l;
1688 
1689  // ** combinatorics
1690  DoCombinatorics(l, tl);
1691 
1692  // ** preselection
1693  l.Select(tl->GetQASelector());
1694 
1695  // *** store QA
1696  for (int i=0;i<l.GetLength();++i)
1697  {
1698  bool acc = false;
1699  bool tag = sel->Accept(l[i]) && AcceptDstCut(l[i]);
1700 
1701  // full selection?
1702  if (fApplyFullSelection>0)
1703  acc = AcceptCandidate(mode, l[i], sel);
1704  // if not, simple mass window selection (and D* cut if applicable)
1705  else acc = tag;
1706 
1707  // for full selection in open mode (=2), trigger w/o detailed cuts are accepted based on mass window only
1708  if ( !acc && fApplyFullSelection==2 && fSTSelmap.find(fSTencode[fSTModeIndex]*1000+mode) == fSTSelmap.end() )
1709  acc = tag;
1710 
1711  if (acc) nacc++;
1712  if (tag) npre++;
1713 
1714  fAnalysis->McTruthMatch(l[i]);
1715 
1716  if (n && AcceptDstCut(l[i]))
1717  {
1718  // in case we only want to store MC truth matched candidates, check!
1719  if (fQAMctOnly && (l[i]->GetMcTruth())==0) continue;
1720 
1721  fQA->qaComp(prefix, l[i], n);
1722  fQA->qaEventShapeShort("es", fEventShape, n);
1723 
1724  // replace PID mult values from event shape by actual counts with individual algos
1725  n->Column("eslnpide", (Int_t) fPidMult_025[0], 0);
1726  n->Column("eslnpidmu",(Int_t) fPidMult_025[1], 0);
1727  n->Column("eslnpidpi",(Int_t) fPidMult_025[2], 0);
1728  n->Column("eslnpidk", (Int_t) fPidMult_025[3], 0);
1729  n->Column("eslnpidp", (Int_t) fPidMult_025[4], 0);
1730 
1731  fQA->qaP4("beam", fIniP4, n);
1732  n->Column("primvx", (Float_t) fPrimVtx.X(), 0.0f);
1733  n->Column("primvy", (Float_t) fPrimVtx.Y(), 0.0f);
1734  n->Column("primvz", (Float_t) fPrimVtx.Z(), 0.0f);
1735  n->Column("primvqa", (Float_t) fPrimVtxQa , 0.0f);
1736 
1737  Float_t mmiss = (fIniP4-(l[i]->P4())).M();
1738  Float_t nsig = (Float_t) fabs(l[i]->Mass()-mean)/sigma;
1739 
1740  n->Column("ev", (Int_t) fEvtCount, 0);
1741  n->Column("num", (Int_t) i, 0);
1742  n->Column("run", (Int_t) fRunNum, 0);
1743  n->Column("mode", (Int_t) fMode, 0);
1744  n->Column("recmode",(Int_t) fRecoilMode, 0);
1745  n->Column("reccnt", (Int_t) fRecoilCnt, 0);
1746  n->Column("mmiss", (Float_t) mmiss, 0.0f);
1747  n->Column("nsig", (Float_t) nsig, 0.0f);
1748  n->Column(prefix+"mean",(Float_t) mean, 0.0f);
1749  n->Column(prefix+"sig", (Float_t) sigma, 0.0f);
1750  n->Column("tag", (Int_t) tag, 0);
1751  n->Column("acc", (Int_t) acc, 0);
1752 
1753  n->DumpData();
1754  }
1755  }
1756 
1757  return nacc;
1758 }
1759 
void SetQAMassWindow(double min, double max)
RhoMassParticleSelector * fEtaSel
int GetDaughterPdg(int idx)
RhoMassParticleSelector * fKs0Sel
double DbMass(TString name)
float fSTInputValues[STMAXCUT]
void qaPi0(TString pre, RhoCandidate *c, RhoTuple *n)
std::map< int, PndSoftTriggerLine * > fSTTriggers
int fVerbose
Definition: poormantracks.C:24
double ChrgPtSumCms() const
Definition: PndEventShape.h:45
RhoMassParticleSelector * GetQASelector()
double ChrgPSumCms() const
Definition: PndEventShape.h:46
friend F32vec4 cos(const F32vec4 &a)
Definition: P4_F32vec4.h:112
RhoEnergyParticleSelector * fEnergySel
const int STMAXCUT
void SetPidAlgoPion(TString algo)
std::map< int, PndSoftTriggerLine * >::iterator TrigIt
Int_t run
Definition: autocutx.C:47
void Cleanup()
Definition: RhoCandList.cxx:62
double Sphericity()
Int_t i
Definition: run_full.C:25
RhoTuple * GetRhoTuple()
TFile * file
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
int tmvavarid[STMAXCUT]
double ChrgPSumLab() const
Definition: PndEventShape.h:39
Double_t val[nBoxes][nFEBox]
Definition: createCalib.C:11
void SetPidAlgoMuon(TString algo)
double NeutEtSumCms() const
Definition: PndEventShape.h:43
std::map< int, double > fSTVarCandArray
int SelectTruePid(RhoCandList &l)
Int_t tag
Definition: crosstag.C:23
Double_t Pt() const
Definition: RhoCandidate.h:211
TLorentzVector s
Definition: Pnd2DStar.C:50
Int_t GetLength() const
Definition: RhoCandList.h:46
TMVA::Reader * tmvaread
int n
Bool_t FillList(RhoCandList &l, TString listkey="All", TString pidTcaNames="", int trackHypothesis=-1)
virtual Bool_t Accept(RhoCandidate *)=0
void SetConfigurationFile(TString fname)
double SumChrgPminCms(double pmin)
void SetPidAlgoElectron(TString algo)
void SetKs0SignalParams(double mean, double sigma)
RhoCandidate * Daughter(Int_t n)
RhoMassParticleSelector * fPi0PreSel
TString cuts[MAX]
Definition: autocutx.C:35
double Circularity()
Double_t sigma[nsteps]
Definition: dedx_bands.C:65
void SetPidAlgoKaon(TString algo)
std::map< int, int > fSTMapListIndex
double r1
double PtSumLab() const
Definition: PndEventShape.h:35
void SetThreshold(double thr)
void qaP4(TString pre, TLorentzVector c, RhoTuple *n, bool skip=false)
void qaKs0(TString pre, RhoCandidate *c, RhoTuple *n)
int NParticles() const
Definition: PndEventShape.h:17
void SetPi0SignalParams(double mean, double sigma)
int SelectPidProb(RhoCandList &l, int pididx, double cut)
__m128 v
Definition: P4_F32vec4.h:4
virtual void Exec(Option_t *opt)
static const double mp
Definition: mzparameters.h:11
void qaMcList(TString pre, RhoCandList &l, RhoTuple *n, int max=10000)
Double_t d0
Definition: checkhelixhit.C:59
int Pic_FED Eff_lEE C()
double cut[MAX]
Definition: autocutx.C:36
int DoCombinatorics(RhoCandList &l, PndSoftTriggerLine *tl)
Int_t mode
Definition: autocutx.C:47
void SetNTag(int mode, const int tag)
double cutval[STMAXCUT]
int idx[MAX]
Definition: autocutx.C:38
void CombineAndAppend(RhoCandList &l1, RhoCandList &l2)
void SetTagNSig(double nsig)
void GetEventInTask()
void Combine(RhoCandList &l1, RhoCandList &l2)
std::map< int, STCutSet > fSTSelmap
void FillVarArray(RhoCandidate *c, int id, Bool_t tmva=false)
Double_t GetPocaVtx(TVector3 &vertex, RhoCandidate *composite)
Definition: RhoVtxPoca.cxx:27
RhoMassParticleSelector * fKs0PreSel
void SetEtaSignalParams(double mean, double sigma)
double Planarity()
bool AcceptCandidate(int mode, RhoCandidate *c, RhoParticleSelectorBase *sel=0)
float_m ok
double GetPocaVtx(RhoCandidate *c, double &dist, double &ctau)
void qaComp(TString pre, RhoCandidate *c, RhoTuple *n, bool covs=false, bool pulls=false)
double PRapmax() const
Definition: PndEventShape.h:28
std::map< TString, int > fSTVarmap
TLorentzVector BoostCms(TLorentzVector in)
void Select(RhoParticleSelectorBase *pidmgr)
TFile * f
Definition: bump_analys.C:12
double NeutESumLab() const
Definition: PndEventShape.h:37
void qaEventShape(TString pre, PndEventShape *evsh, RhoTuple *n)
TLorentzVector P4() const
Definition: RhoCandidate.h:195
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)
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
Double_t M() const
Int_t Remove(RhoCandidate *)
TString name
void SetTagMode(int mode, bool tag=true)
double fSTVarEvArray[STMAXEVVARS]
Double_t P() const
void SetQAMode(int mode, bool qa=true)
void DumpData()
Definition: RhoTuple.cxx:391
void qaEventShapeShort(TString pre, PndEventShape *evsh, RhoTuple *n)
Int_t RemoveFamily(RhoCandidate *)
double PmaxLab() const
Definition: PndEventShape.h:22
void McMatchAllowPhotos(int maxn=1, double thresh=0.05)
Definition: PndAnalysis.h:62
double PminCms() const
Definition: PndEventShape.h:25
int nsig
Definition: toy_core.C:46
double Thrust(int Nmax=4)
int TagMode(PndSoftTriggerLine *tl, int &npre)
TString fSTenames[]
RhoMassParticleSelector * fEtaPreSel
void SetTagNSigMode(int mode, double nsig)
double Ptmin() const
Definition: PndEventShape.h:27
Int_t NDaughters() const
int MultNeutEminLab(double emin)
int MultChrgPminLab(double pmin)
ClassImp(PndAnaContFact)
double FoxWolfMomR(int order)
RhoMassParticleSelector * fPi0Sel
RhoMomentumParticleSelector * fMomentumSel
TTree * GetInternalTree()
Definition: RhoTuple.h:207
int MultPidProb(RhoCandList &l, int pididx, double prob)
Bool_t McTruthMatch(RhoCandidate *cand, Int_t level=2, bool verbose=false)
Double_t mean[nsteps]
Definition: dedx_bands.C:65
double NeutESumCms() const
Definition: PndEventShape.h:44
TString fSTOps[5]
RhoMassParticleSelector * GetSelector()
void SetInitialPidCut(double cut)
void SetPidAlgoProton(TString algo)
int NNeutral() const
Definition: PndEventShape.h:19
int NCharged() const
Definition: PndEventShape.h:18
bool AcceptDstCut(RhoCandidate *c)
double PminLab() const
Definition: PndEventShape.h:24
int fSTPidIndex[16]
std::vector< int > fSTencode
void SetTagActive(bool ac=true)
TString fSTVarSuff[]
void SetWriteQA(bool qa=true)
if(fWindowIsBox)
void SetMeanSigma(double mean, double sig)
void CombineList(RhoCandList &l, int mothpdg, int amothpdg, std::vector< int > &idx, std::vector< int > &aidx, bool cc=false)
double GetVarValue(RhoCandidate *c, int id)
Double_t mult
const int STMAXEVVARS
double GetPidInfo(int hypo)
double r2
TVector3 P3() const
Definition: RhoCandidate.h:199
double PmaxCms() const
Definition: PndEventShape.h:23
void ApplyFullSelection(int selmode=1)
double Ptmax() const
Definition: PndEventShape.h:26
virtual Bool_t Accept(RhoCandidate *b)
int SplitString(TString s, TString delim, TString *toks, int maxtoks)
void GetAngles(RhoCandidate *c, double &oang, double &decang)
double Aplanarity()