FairRoot/PandaRoot
Functions
QA/dpm3/ana_complete.C File Reference
#include "../auxi.C"

Go to the source code of this file.

Functions

int SelectTruePid (PndAnalysis *ana, RhoCandList &l)
 
void printCand (RhoCandidate *c)
 
void countDoubles (RhoCandList &l, int &n1, int &n2, int &n3)
 
int ana_complete (int nevts=0)
 

Function Documentation

int ana_complete ( int  nevts = 0)

Definition at line 57 of file QA/dpm3/ana_complete.C.

References RhoKinFitter::Add4MomConstraint(), RhoKinFitter::AddMassConstraint(), CloseGeoManager(), RhoCandList::Combine(), countDoubles(), ctime, Double_t, PndAnalysis::FillList(), RhoKinFitter::Fit(), RhoFitterBase::Fit(), fRun, RhoFitterBase::GetChi2(), PndAnalysis::GetEntries(), PndAnalysis::GetEvent(), RhoCandList::GetLength(), RhoFitterBase::GetProb(), hvpos, i, RhoCandidate::M(), PndAnalysis::McTruthMatch(), muminus, muplus, out, P, piminus, piplus, PndAnalysis::PndAnalysis(), RhoCandidate::Pos(), rtdb, rtime, RhoCandList::Select(), SelectTruePid(), RhoCandList::SetType(), and TString.

58 {
59  TDatabasePDG::Instance()->AddParticle("pbarpSystem","pbarpSystem",1.9,kFALSE,0.1,0,"",88888);
60  TStopwatch fTimer;
61  // *** some variables
62  int i=0,j=0, k=0, l=0;
63  gStyle->SetOptFit(1011);
64 
65  // *** the output file for FairRunAna
66  TString OutFile="output.root";
67 
68  // *** the files coming from the simulation
69  TString inPidFile = "dpm_qa_pid.root"; // this file contains the PndPidCandidates and McTruth
70  TString inParFile = "dpm_qa_par.root";
71 
72  // *** PID table with selection thresholds; can be modified by the user
73  TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";
74 
75  // *** initialization
76  FairLogger::GetLogger()->SetLogToFile(kFALSE);
77  FairRunAna* fRun = new FairRunAna();
78  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
79  fRun->SetInputFile(inPidFile);
80 
81  // *** setup parameter database
82  FairParRootFileIo* parIO = new FairParRootFileIo();
83  parIO->open(inParFile);
84  FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo();
85  parIOPid->open(pidParFile.Data(),"in");
86 
87  rtdb->setFirstInput(parIO);
88  rtdb->setSecondInput(parIOPid);
89  rtdb->setOutput(parIO);
90 
91  fRun->SetOutputFile(OutFile);
92  fRun->Init();
93 
94  // *** create an output file for all histograms
95  TFile *out = TFile::Open("output_ana.root","RECREATE");
96 
97  // *** create some histograms
98  TH1F *hmomtrk = new TH1F("hmomtrk","track momentum (all)",200,0,5);
99  TH1F *hthttrk = new TH1F("hthttrk","track theta (all)",200,0,3.1415);
100 
101  TH1F *hjpsim_all = new TH1F("hjpsim_all","J/#psi mass (all)",200,0,4.5);
102  TH1F *hpsim_all = new TH1F("hpsim_all","#psi(2S) mass (all)",200,0,5);
103 
104  TH1F *hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4.5);
105  TH1F *hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5);
106 
107  TH1F *hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (tight pid)",200,0,4.5);
108  TH1F *hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (tight pid)",200,0,5);
109 
110  TH1F *hjpsim_trpid = new TH1F("hjpsim_trpid","J/#psi mass (true pid)",200,0,4.5);
111  TH1F *hpsim_trpid = new TH1F("hpsim_trpid","#psi(2S) mass (true pid)",200,0,5);
112 
113 
114  TH1F *hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4.5);
115  TH1F *hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass (full truth match)",200,0,5);
116 
117  TH1F *hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4.5);
118  TH1F *hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5);
119 
120  TH1F *hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2);
121  TH1F *hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2);
122 
123 
124  TH1F *hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass (vertex fit)",200,0,4.5);
125  TH1F *hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4.5);
126  TH1F *hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (mass constraint fit)",200,0,4.5);
127 
128  TH1F *hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf", "J/#psi: #chi^{2} vertex fit",100,0,10);
129  TH1F *hpsi_chi2_4c = new TH1F("hpsi_chi2_4c", "#psi(2S): #chi^{2} 4C fit",100,0,250);
130  TH1F *hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf", "J/#psi: #chi^{2} mass fit",100,0,10);
131 
132  TH1F *hjpsi_prob_vf = new TH1F("hjpsi_prob_vf", "J/#psi: Prob vertex fit",100,0,1);
133  TH1F *hpsi_prob_4c = new TH1F("hpsi_prob_4c", "#psi(2S): Prob 4C fit",100,0,1);
134  TH1F *hjpsi_prob_mf = new TH1F("hjpsi_prob_mf", "J/#psi: Prob mass fit",100,0,1);
135 
136  TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2);
137 
138  //
139  // Now the analysis stuff comes...
140  //
141 
142 
143  // *** the data reader object
144  PndAnalysis* theAnalysis = new PndAnalysis();
145  if (nevts==0) nevts= theAnalysis->GetEntries();
146 
147  // *** RhoCandLists for the analysis
148  RhoCandList chrg, muplus, muminus, piplus, piminus, jpsi, psi2s;
149 
150  // *** Mass selector for the jpsi cands
151  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
152  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
153 
154  // *** the lorentz vector of the initial psi(2S)
155  TLorentzVector ini(0, 0, 6.231552, 7.240065);
156 
157  // ***
158  // the event loop
159  // ***
160 
161  int cntdbltrk=0, cntdblmc=0, cntdblboth=0, cnttrk=0, cnt_dbl_jpsi=0, cnt_dbl_psip=0;
162 
163  while (theAnalysis->GetEvent() && i++<nevts)
164  {
165  if ((i%100)==0) cout<<"evt " << i << endl;
166 
167  // *** Select with no PID info ('All'); type and mass are set
168  theAnalysis->FillList(chrg, "Charged");
169  theAnalysis->FillList(muplus, "MuonAllPlus");
170  theAnalysis->FillList(muminus, "MuonAllMinus");
171  theAnalysis->FillList(piplus, "PionAllPlus");
172  theAnalysis->FillList(piminus, "PionAllMinus");
173 
174  // *** momentum and theta histograms
175  for (j=0;j<muplus.GetLength();++j)
176  {
177  hmomtrk->Fill(muplus[j]->P());
178  hthttrk->Fill(muplus[j]->P4().Theta());
179  }
180  for (j=0;j<muminus.GetLength();++j)
181  {
182  hmomtrk->Fill(muminus[j]->P());
183  hthttrk->Fill(muminus[j]->P4().Theta());
184  }
185 
186  cnttrk += chrg.GetLength();
187 
188  int n1, n2, n3;
189 
190  countDoubles(chrg,n1,n2,n3);
191  cntdbltrk += n1;
192  cntdblmc += n2;
193  cntdblboth += n3;
194 
195  // *** combinatorics for J/psi -> mu+ mu-
196  jpsi.Combine(muplus, muminus);
197 
198 
199  // ***
200  // *** do the TRUTH MATCH for jpsi
201  // ***
202  jpsi.SetType(443);
203 
204  int nm = 0;
205  for (j=0;j<jpsi.GetLength();++j)
206  {
207  hjpsim_all->Fill( jpsi[j]->M() );
208 
209  if (theAnalysis->McTruthMatch(jpsi[j]))
210  {
211  nm++;
212  hjpsim_ftm->Fill( jpsi[j]->M() );
213  hjpsim_diff->Fill( jpsi[j]->GetMcTruth()->M() - jpsi[j]->M() );
214  }
215  else
216  hjpsim_nm->Fill( jpsi[j]->M() );
217  }
218 
219  if (nm>1) cnt_dbl_jpsi++;
220  // ***
221  // *** do VERTEX FIT (J/psi)
222  // ***
223  for (j=0;j<jpsi.GetLength();++j)
224  {
225  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
226  vtxfitter.Fit();
227 
228  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
229  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
230  hjpsi_chi2_vf->Fill(chi2_vtx);
231  hjpsi_prob_vf->Fill(prob_vtx);
232 
233  if ( prob_vtx > 0.01 ) // when good enough, fill some histos
234  {
235  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
236  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
237 
238  hjpsim_vf->Fill(jfit->M());
239  hvpos->Fill(jVtx.X(),jVtx.Y());
240  }
241  }
242 
243  // *** some rough mass selection
244  jpsi.Select(jpsiMassSel);
245 
246  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
247  psi2s.Combine(jpsi, piplus, piminus);
248 
249 
250  // ***
251  // *** do the TRUTH MATCH for psi(2S)
252  // ***
253  psi2s.SetType(88888);
254 
255  nm = 0;
256  for (j=0;j<psi2s.GetLength();++j)
257  {
258  hpsim_all->Fill( psi2s[j]->M() );
259 
260  if (theAnalysis->McTruthMatch(psi2s[j]))
261  {
262  nm++;
263  hpsim_ftm->Fill( psi2s[j]->M() );
264  hpsim_diff->Fill( psi2s[j]->GetMcTruth()->M() - psi2s[j]->M() );
265  }
266  else
267  hpsim_nm->Fill( psi2s[j]->M() );
268  }
269  if (nm>1) cnt_dbl_psip++;
270 
271 
272  // ***
273  // *** do 4C FIT (initial psi(2S) system)
274  // ***
275  for (j=0;j<psi2s.GetLength();++j)
276  {
277  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
278  fitter.Add4MomConstraint(ini); // set 4 constraint
279  fitter.Fit(); // do fit
280 
281  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
282  double prob_4c = fitter.GetProb(); // access probability of fit
283  hpsi_chi2_4c->Fill(chi2_4c);
284  hpsi_prob_4c->Fill(prob_4c);
285 
286  if ( prob_4c > 0.01 ) // when good enough, fill some histo
287  {
288  RhoCandidate *jfit = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
289 
290  hjpsim_4cf->Fill(jfit->M());
291  }
292  }
293 
294 
295  // ***
296  // *** do MASS CONSTRAINT FIT (J/psi)
297  // ***
298  for (j=0;j<jpsi.GetLength();++j)
299  {
300  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
301  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
302  mfitter.Fit(); // do fit
303 
304  double chi2_m = mfitter.GetChi2(); // get chi2 of fit
305  double prob_m = mfitter.GetProb(); // access probability of fit
306  hjpsi_chi2_mf->Fill(chi2_m);
307  hjpsi_prob_mf->Fill(prob_m);
308 
309  if ( prob_m > 0.01 ) // when good enough, fill some histo
310  {
311  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
312  hjpsim_mcf->Fill(jfit->M());
313  }
314  }
315 
316 
317  // ***
318  // *** TRUE PID combinatorics
319  // ***
320 
321  // *** do MC truth match for PID type
322  SelectTruePid(theAnalysis, muplus);
323  SelectTruePid(theAnalysis, muminus);
324  SelectTruePid(theAnalysis, piplus);
325  SelectTruePid(theAnalysis, piminus);
326 
327  // *** all combinatorics again with true PID
328  jpsi.Combine(muplus, muminus);
329  for (j=0;j<jpsi.GetLength();++j) hjpsim_trpid->Fill( jpsi[j]->M() );
330  jpsi.Select(jpsiMassSel);
331 
332  psi2s.Combine(jpsi, piplus, piminus);
333  for (j=0;j<psi2s.GetLength();++j) hpsim_trpid->Fill( psi2s[j]->M() );
334 
335 
336  // ***
337  // *** LOOSE PID combinatorics
338  // ***
339 
340  // *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection
341  theAnalysis->FillList(muplus, "MuonLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc;PidAlgoMdtHardCuts");
342  theAnalysis->FillList(muminus, "MuonLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc;PidAlgoMdtHardCuts");
343  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc");
344  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc;PidAlgoDisc");
345 
346  jpsi.Combine(muplus, muminus);
347  for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j]->M() );
348  jpsi.Select(jpsiMassSel);
349 
350  psi2s.Combine(jpsi, piplus, piminus);
351  for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j]->M() );
352 
353 
354  // ***
355  // *** TIGHT PID combinatorics
356  // ***
357 
358  // *** and again with PidAlgoMvd;PidAlgoStt and tight selection
359  theAnalysis->FillList(muplus, "MuonTightPlus", "PidAlgoMdtHardCuts");
360  theAnalysis->FillList(muminus, "MuonTightMinus", "PidAlgoMdtHardCuts");
361  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
362  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
363 
364  jpsi.Combine(muplus, muminus);
365  for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j]->M() );
366  jpsi.Select(jpsiMassSel);
367 
368  psi2s.Combine(jpsi, piplus, piminus);
369  for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j]->M() );
370 
371  }
372 
373  // *** write out all the histos
374  out->cd();
375 
376  hmomtrk->Write();
377  hthttrk->Write();
378 
379  hjpsim_all->Write();
380  hpsim_all->Write();
381  hjpsim_lpid->Write();
382  hpsim_lpid->Write();
383  hjpsim_tpid->Write();
384  hpsim_tpid->Write();
385  hjpsim_trpid->Write();
386  hpsim_trpid->Write();
387 
388  hjpsim_ftm->Write();
389  hpsim_ftm->Write();
390  hjpsim_nm->Write();
391  hpsim_nm->Write();
392 
393  hpsim_diff->Write();
394  hjpsim_diff->Write();
395 
396  hjpsim_vf->Write();
397  hjpsim_4cf->Write();
398  hjpsim_mcf->Write();
399 
400  hjpsi_chi2_vf->Write();
401  hpsi_chi2_4c->Write();
402  hjpsi_chi2_mf->Write();
403 
404  hjpsi_prob_vf->Write();
405  hpsi_prob_4c->Write();
406  hjpsi_prob_mf->Write();
407 
408  hvpos->Write();
409 
410  out->Save();
411 
412  // Extract the maximal used memory an add is as Dart measurement
413  // This line is filtered by CTest and the value send to CDash
414  FairSystemInfo sysInfo;
415  Float_t maxMemory=sysInfo.GetMaxMemory();
416  cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
417  cout << maxMemory;
418  cout << "</DartMeasurement>" << endl;
419 
420  fTimer.Stop();
421  Double_t rtime = fTimer.RealTime();
422  Double_t ctime = fTimer.CpuTime();
423 
424  Float_t cpuUsage=ctime/rtime;
425  cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
426  cout << cpuUsage;
427  cout << "</DartMeasurement>" << endl;
428 
429  cout << endl;
430  cout << "Real time " << rtime << " s, CPU time " << ctime
431  << "s" << endl;
432  cout << "CPU usage " << cpuUsage*100. << "%" << endl;
433  cout << "Max Memory " << maxMemory << " MB" << endl;
434 
435  cout << "Macro finished successfully." << endl;
436 
437  CloseGeoManager();
438  return 0;
439 
440 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
countDoubles
void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
Definition: QA/dpm3/ana_complete.C:32
i
Int_t i
Definition: run_full.C:25
RhoCandidate::Pos
TVector3 Pos() const
Definition: RhoCandidate.h:186
RhoCandList::GetLength
Int_t GetLength() const
Definition: RhoCandList.h:46
PndAnalysis::FillList
Bool_t FillList(RhoCandList &l, TString listkey="All", TString pidTcaNames="", int trackHypothesis=-1)
Definition: PndAnalysis.cxx:319
piplus
CandList piplus
Definition: softtrigger_kin5.C:144
SelectTruePid
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
Definition: QA/dpm3/ana_complete.C:9
CloseGeoManager
void CloseGeoManager()
Definition: QA/auxi.C:11
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
muplus
CandList muplus
Definition: softtrigger_kin5.C:144
RhoCandList::Combine
void Combine(RhoCandList &l1, RhoCandList &l2)
Definition: RhoCandList.cxx:289
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
RhoCandList::Select
void Select(RhoParticleSelectorBase *pidmgr)
Definition: RhoCandList.cxx:1359
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
RhoCandList::SetType
void SetType(const TParticlePDG *pdt, int start=0)
Definition: RhoCandList.cxx:264
out
TFile * out
Definition: reco_muo.C:20
RhoCandidate::M
Double_t M() const
Definition: RhoCandidate.cxx:1636
ctime
Double_t ctime
Definition: hit_dirc.C:114
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
piminus
CandList piminus
Definition: softtrigger_kin5.C:144
PndAnalysis::McTruthMatch
Bool_t McTruthMatch(RhoCandidate *cand, Int_t level=2, bool verbose=false)
Definition: PndAnalysis.cxx:932
PndAnalysis::GetEvent
Int_t GetEvent(Int_t n=-1)
Definition: PndAnalysis.cxx:257
muminus
CandList muminus
Definition: softtrigger_kin5.C:144
rtime
Double_t rtime
Definition: hit_dirc.C:113
hvpos
TH2F * hvpos
Definition: plot_eta_c_stt.C:36
void countDoubles ( RhoCandList &  l,
int &  n1,
int &  n2,
int &  n3 
)

Definition at line 32 of file QA/dpm3/ana_complete.C.

References d, fabs(), RhoCandList::GetLength(), and i.

Referenced by ana_complete().

33 {
34  int n_smc = 0;
35  int n_strk = 0;
36  int n_both = 0;
37  double d = 0.00001;
38 
39  for (int i=0;i<l.GetLength()-1;++i)
40  {
41  for (int j=i+1;j<l.GetLength();++j)
42  {
43  TLorentzVector dl = l[i]->P4() - l[j]->P4();
44 
45  bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth());
46  bool chktrk = (fabs(dl.X())<d) && (fabs(dl.Y())<d) && (fabs(dl.Z())<d) && (fabs(dl.E())<d);
47  if (chkmc) n_smc++;
48  if (chktrk) n_strk++;
49  if (chktrk && chkmc) n_both++;
50  }
51  }
52  n1 = n_strk;
53  n2 = n_smc;
54  n3 = n_both;
55 }
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
i
Int_t i
Definition: run_full.C:25
RhoCandList::GetLength
Int_t GetLength() const
Definition: RhoCandList.h:46
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
void printCand ( RhoCandidate *  c)

Definition at line 25 of file QA/dpm3/ana_complete.C.

References RhoCandidate::P4(), and RhoCandidate::PdgCode().

26 {
27  TLorentzVector lv=c->P4();
28 
29  cout <<c->PdgCode()<<" ("<<lv.X()<<"/"<<lv.Y()<<"/"<<lv.Z()<<"/"<<lv.E()<<")"<<endl;
30 }
RhoCandidate::P4
TLorentzVector P4() const
Definition: RhoCandidate.h:195
int SelectTruePid ( PndAnalysis *  ana,
RhoCandList &  l 
)

Definition at line 9 of file QA/dpm3/ana_complete.C.

References RhoCandList::GetLength(), PndAnalysis::McTruthMatch(), and RhoCandList::Remove().

Referenced by ana_complete().

10 {
11  int removed = 0;
12 
13  for (int ii=l.GetLength()-1;ii>=0;--ii)
14  {
15  if ( !(ana->McTruthMatch(l[ii])) )
16  {
17  l.Remove(l[ii]);
18  removed++;
19  }
20  }
21 
22  return removed;
23 }
RhoCandList::GetLength
Int_t GetLength() const
Definition: RhoCandList.h:46
RhoCandList::Remove
Int_t Remove(RhoCandidate *)
Definition: RhoCandList.cxx:131
PndAnalysis::McTruthMatch
Bool_t McTruthMatch(RhoCandidate *cand, Int_t level=2, bool verbose=false)
Definition: PndAnalysis.cxx:932
Generated on Sun Aug 18 2019 03:02:20 for FairRoot/PandaRoot by   doxygen 1.8.8