FairRoot/PandaRoot
Functions
production/dayone-2017/fullsim/ana_day1.C File Reference

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_day1 (int nevts=0)
 

Function Documentation

int ana_day1 ( int  nevts = 0)

Definition at line 56 of file production/dayone-2017/fullsim/ana_day1.C.

References RhoKinFitter::Add4MomConstraint(), PndMasterRunAna::AddFriend(), RhoKinFitter::AddMassConstraint(), RhoCandList::Combine(), countDoubles(), PndAnalysis::FillList(), PndMasterRunAna::Finish(), RhoKinFitter::Fit(), RhoFitterBase::Fit(), fRun, RhoFitterBase::GetChi2(), PndAnalysis::GetEntries(), PndAnalysis::GetEvent(), RhoCandList::GetLength(), RhoFitterBase::GetProb(), gGeoManager, hvpos, i, RhoCandidate::M(), PndAnalysis::McTruthMatch(), muminus, muplus, out, output, P, piminus, piplus, PndAnalysis::PndAnalysis(), RhoCandidate::Pos(), RhoCandList::Select(), SelectTruePid(), PndMasterRunAna::SetInput(), PndMasterRunAna::SetOutput(), PndMasterRunAna::SetParamAsciiFile(), RhoCandList::SetType(), PndMasterRunAna::Setup(), and TString.

57 {
58  //-----User Settings:------------------------------------------------------
59  TString parAsciiFile = "all.par";
60  TString prefix = "evtday1";
61  TString input = "psi2s_Jpsi2pi_Jpsi_mumu.dec";
62  TString output = "ana";
63  TString friend1 = "pidideal";
64  TString friend2 = "";
65  TString friend3 = "";
66  TString friend4 = "";
67 
68  // ----- Initial Settings --------------------------------------------
69  PndMasterRunAna *fRun= new PndMasterRunAna();
70  fRun->SetInput(input);
71  fRun->SetOutput(output);
72  fRun->AddFriend(friend1);
73  fRun->AddFriend(friend2);
74  fRun->AddFriend(friend3);
75  fRun->AddFriend(friend4);
76  fRun->SetParamAsciiFile(parAsciiFile);
77  fRun->Setup(prefix);
78 
79  // *** some variables
80  int i=0,j=0, k=0, l=0;
81  gStyle->SetOptFit(1011);
82 
83  fRun->Init();
84 
85  // *** create an output file for all histograms
86  TFile *out = TFile::Open(prefix+"_output_ana.root","RECREATE");
87 
88  // *** create some histograms
89  TH1F *hmomtrk = new TH1F("hmomtrk","track momentum (all)",200,0,5);
90  TH1F *hthttrk = new TH1F("hthttrk","track theta (all)",200,0,3.1415);
91 
92  TH1F *hjpsim_all = new TH1F("hjpsim_all","J/#psi mass (all)",200,0,4.5);
93  TH1F *hpsim_all = new TH1F("hpsim_all","#psi(2S) mass (all)",200,0,5);
94 
95  TH1F *hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4.5);
96  TH1F *hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5);
97 
98  TH1F *hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (tight pid)",200,0,4.5);
99  TH1F *hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (tight pid)",200,0,5);
100 
101  TH1F *hjpsim_trpid = new TH1F("hjpsim_trpid","J/#psi mass (true pid)",200,0,4.5);
102  TH1F *hpsim_trpid = new TH1F("hpsim_trpid","#psi(2S) mass (true pid)",200,0,5);
103 
104 
105  TH1F *hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4.5);
106  TH1F *hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass (full truth match)",200,0,5);
107 
108  TH1F *hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4.5);
109  TH1F *hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5);
110 
111  TH1F *hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2);
112  TH1F *hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2);
113 
114 
115  TH1F *hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass (vertex fit)",200,0,4.5);
116  TH1F *hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4.5);
117  TH1F *hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (mass constraint fit)",200,0,4.5);
118 
119  TH1F *hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf", "J/#psi: #chi^{2} vertex fit",100,0,10);
120  TH1F *hpsi_chi2_4c = new TH1F("hpsi_chi2_4c", "#psi(2S): #chi^{2} 4C fit",100,0,250);
121  TH1F *hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf", "J/#psi: #chi^{2} mass fit",100,0,10);
122 
123  TH1F *hjpsi_prob_vf = new TH1F("hjpsi_prob_vf", "J/#psi: Prob vertex fit",100,0,1);
124  TH1F *hpsi_prob_4c = new TH1F("hpsi_prob_4c", "#psi(2S): Prob 4C fit",100,0,1);
125  TH1F *hjpsi_prob_mf = new TH1F("hjpsi_prob_mf", "J/#psi: Prob mass fit",100,0,1);
126 
127  TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2);
128 
129  //
130  // Now the analysis stuff comes...
131  //
132 
133 
134  // *** the data reader object
135  PndAnalysis* theAnalysis = new PndAnalysis();
136  if (nevts==0) nevts= theAnalysis->GetEntries();
137 
138  // *** RhoCandLists for the analysis
139  RhoCandList chrg, muplus, muminus, piplus, piminus, jpsi, psi2s;
140 
141  // *** Mass selector for the jpsi cands
142  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
143  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
144 
145  // *** the lorentz vector of the initial psi(2S)
146  TLorentzVector ini(0, 0, 6.231552, 7.240065);
147 
148  // ***
149  // the event loop
150  // ***
151 
152  int cntdbltrk=0, cntdblmc=0, cntdblboth=0, cnttrk=0, cnt_dbl_jpsi=0, cnt_dbl_psip=0;
153 
154  while (theAnalysis->GetEvent() && i++<nevts)
155  {
156  if ((i%100)==0) cout<<"evt " << i << endl;
157 
158  // *** Select with no PID info ('All'); type and mass are set
159  theAnalysis->FillList(chrg, "Charged");
160  theAnalysis->FillList(muplus, "MuonAllPlus");
161  theAnalysis->FillList(muminus, "MuonAllMinus");
162  theAnalysis->FillList(piplus, "PionAllPlus");
163  theAnalysis->FillList(piminus, "PionAllMinus");
164 
165  // *** momentum and theta histograms
166  for (j=0;j<muplus.GetLength();++j)
167  {
168  hmomtrk->Fill(muplus[j]->P());
169  hthttrk->Fill(muplus[j]->P4().Theta());
170  }
171  for (j=0;j<muminus.GetLength();++j)
172  {
173  hmomtrk->Fill(muminus[j]->P());
174  hthttrk->Fill(muminus[j]->P4().Theta());
175  }
176 
177  cnttrk += chrg.GetLength();
178 
179  int n1, n2, n3;
180 
181  countDoubles(chrg,n1,n2,n3);
182  cntdbltrk += n1;
183  cntdblmc += n2;
184  cntdblboth += n3;
185 
186  // *** combinatorics for J/psi -> mu+ mu-
187  jpsi.Combine(muplus, muminus);
188 
189 
190  // ***
191  // *** do the TRUTH MATCH for jpsi
192  // ***
193  jpsi.SetType(443);
194 
195  int nm = 0;
196  for (j=0;j<jpsi.GetLength();++j)
197  {
198  hjpsim_all->Fill( jpsi[j]->M() );
199 
200  if (theAnalysis->McTruthMatch(jpsi[j]))
201  {
202  nm++;
203  hjpsim_ftm->Fill( jpsi[j]->M() );
204  hjpsim_diff->Fill( jpsi[j]->GetMcTruth()->M() - jpsi[j]->M() );
205  }
206  else
207  hjpsim_nm->Fill( jpsi[j]->M() );
208  }
209 
210  if (nm>1) cnt_dbl_jpsi++;
211  // ***
212  // *** do VERTEX FIT (J/psi)
213  // ***
214  for (j=0;j<jpsi.GetLength();++j)
215  {
216  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
217  vtxfitter.Fit();
218 
219  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
220  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
221  hjpsi_chi2_vf->Fill(chi2_vtx);
222  hjpsi_prob_vf->Fill(prob_vtx);
223 
224  if ( prob_vtx > 0.01 ) // when good enough, fill some histos
225  {
226  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
227  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
228 
229  hjpsim_vf->Fill(jfit->M());
230  hvpos->Fill(jVtx.X(),jVtx.Y());
231  }
232  }
233 
234  // *** some rough mass selection
235  jpsi.Select(jpsiMassSel);
236 
237  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
238  psi2s.Combine(jpsi, piplus, piminus);
239 
240 
241  // ***
242  // *** do the TRUTH MATCH for psi(2S)
243  // ***
244  psi2s.SetType(88888);
245 
246  nm = 0;
247  for (j=0;j<psi2s.GetLength();++j)
248  {
249  hpsim_all->Fill( psi2s[j]->M() );
250 
251  if (theAnalysis->McTruthMatch(psi2s[j]))
252  {
253  nm++;
254  hpsim_ftm->Fill( psi2s[j]->M() );
255  hpsim_diff->Fill( psi2s[j]->GetMcTruth()->M() - psi2s[j]->M() );
256  }
257  else
258  hpsim_nm->Fill( psi2s[j]->M() );
259  }
260  if (nm>1) cnt_dbl_psip++;
261 
262 
263  // ***
264  // *** do 4C FIT (initial psi(2S) system)
265  // ***
266  for (j=0;j<psi2s.GetLength();++j)
267  {
268  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
269  fitter.Add4MomConstraint(ini); // set 4 constraint
270  fitter.Fit(); // do fit
271 
272  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
273  double prob_4c = fitter.GetProb(); // access probability of fit
274  hpsi_chi2_4c->Fill(chi2_4c);
275  hpsi_prob_4c->Fill(prob_4c);
276 
277  if ( prob_4c > 0.01 ) // when good enough, fill some histo
278  {
279  RhoCandidate *jfit = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
280 
281  hjpsim_4cf->Fill(jfit->M());
282  }
283  }
284 
285 
286  // ***
287  // *** do MASS CONSTRAINT FIT (J/psi)
288  // ***
289  for (j=0;j<jpsi.GetLength();++j)
290  {
291  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
292  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
293  mfitter.Fit(); // do fit
294 
295  double chi2_m = mfitter.GetChi2(); // get chi2 of fit
296  double prob_m = mfitter.GetProb(); // access probability of fit
297  hjpsi_chi2_mf->Fill(chi2_m);
298  hjpsi_prob_mf->Fill(prob_m);
299 
300  if ( prob_m > 0.01 ) // when good enough, fill some histo
301  {
302  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
303  hjpsim_mcf->Fill(jfit->M());
304  }
305  }
306 
307 
308  // ***
309  // *** TRUE PID combinatorics
310  // ***
311 
312  // *** do MC truth match for PID type
313  SelectTruePid(theAnalysis, muplus);
314  SelectTruePid(theAnalysis, muminus);
315  SelectTruePid(theAnalysis, piplus);
316  SelectTruePid(theAnalysis, piminus);
317 
318  // *** all combinatorics again with true PID
319  jpsi.Combine(muplus, muminus);
320  for (j=0;j<jpsi.GetLength();++j) hjpsim_trpid->Fill( jpsi[j]->M() );
321  jpsi.Select(jpsiMassSel);
322 
323  psi2s.Combine(jpsi, piplus, piminus);
324  for (j=0;j<psi2s.GetLength();++j) hpsim_trpid->Fill( psi2s[j]->M() );
325 
326 
327  // ***
328  // *** LOOSE PID combinatorics
329  // ***
330 
331  // *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection
332  theAnalysis->FillList(muplus, "MuonLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
333  theAnalysis->FillList(muminus, "MuonLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
334  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
335  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
336 
337  jpsi.Combine(muplus, muminus);
338  for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j]->M() );
339  jpsi.Select(jpsiMassSel);
340 
341  psi2s.Combine(jpsi, piplus, piminus);
342  for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j]->M() );
343 
344 
345  // ***
346  // *** TIGHT PID combinatorics
347  // ***
348 
349  // *** and again with PidAlgoMvd;PidAlgoStt and tight selection
350  theAnalysis->FillList(muplus, "MuonTightPlus", "PidAlgoMdtHardCuts");
351  theAnalysis->FillList(muminus, "MuonTightMinus", "PidAlgoMdtHardCuts");
352  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
353  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
354 
355  jpsi.Combine(muplus, muminus);
356  for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j]->M() );
357  jpsi.Select(jpsiMassSel);
358 
359  psi2s.Combine(jpsi, piplus, piminus);
360  for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j]->M() );
361 
362  }
363 
364  // *** write out all the histos
365  out->cd();
366 
367  hmomtrk->Write();
368  hthttrk->Write();
369 
370  hjpsim_all->Write();
371  hpsim_all->Write();
372  hjpsim_lpid->Write();
373  hpsim_lpid->Write();
374  hjpsim_tpid->Write();
375  hpsim_tpid->Write();
376  hjpsim_trpid->Write();
377  hpsim_trpid->Write();
378 
379  hjpsim_ftm->Write();
380  hpsim_ftm->Write();
381  hjpsim_nm->Write();
382  hpsim_nm->Write();
383 
384  hpsim_diff->Write();
385  hjpsim_diff->Write();
386 
387  hjpsim_vf->Write();
388  hjpsim_4cf->Write();
389  hjpsim_mcf->Write();
390 
391  hjpsi_chi2_vf->Write();
392  hpsi_chi2_4c->Write();
393  hjpsi_chi2_mf->Write();
394 
395  hjpsi_prob_vf->Write();
396  hpsi_prob_4c->Write();
397  hjpsi_prob_mf->Write();
398 
399  hvpos->Write();
400 
401  out->Save();
402 
403  fRun->Finish();
404  if (gROOT->GetVersionInt() >= 60602) {
405  gGeoManager->GetListOfVolumes()->Delete();
406  gGeoManager->GetListOfShapes()->Delete();
407  delete gGeoManager;
408  }
409  return 0;
410 
411 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
SelectTruePid
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
Definition: master/ana_complete.C:8
PndMasterRunAna
Class for the master reconstruction chain.
Definition: PndMasterRunAna.h:21
i
Int_t i
Definition: run_full.C:25
PndMasterRunAna::Setup
Bool_t Setup(TString outprefix="")
Initial setup.
Definition: PndMasterRunAna.cxx:48
RhoCandidate::Pos
TVector3 Pos() const
Definition: RhoCandidate.h:186
RhoCandList::GetLength
Int_t GetLength() const
Definition: RhoCandList.h:46
hvpos
TH2F * hvpos
Definition: plot_eta_c.C:47
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
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
PndMasterRunAna::AddFriend
void AddFriend(TString par)
Setter of friend root files.
Definition: PndMasterRunAna.h:123
output
FairParRootFileIo * output
Definition: sim_emc_apd.C:120
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
muplus
CandList muplus
Definition: softtrigger_kin5.C:144
PndMasterRunAna::SetInput
void SetInput(TString par)
Input of the macro.
Definition: PndMasterRunAna.h:103
RhoCandList::Combine
void Combine(RhoCandList &l1, RhoCandList &l2)
Definition: RhoCandList.cxx:289
RhoCandList::Select
void Select(RhoParticleSelectorBase *pidmgr)
Definition: RhoCandList.cxx:1359
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
countDoubles
void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
Definition: master/ana_complete.C:31
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
PndMasterRunAna::Finish
void Finish()
Final diagnostics.
Definition: PndMasterRunAna.cxx:205
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
PndMasterRunAna::SetOutput
void SetOutput(TString par)
Tag of the output file of the macro.
Definition: PndMasterRunAna.h:108
muminus
CandList muminus
Definition: softtrigger_kin5.C:144
PndMasterRunAna::SetParamAsciiFile
void SetParamAsciiFile(TString par)
Setter of the parameter ascii file.
Definition: PndMasterRunAna.h:118
void countDoubles ( RhoCandList &  l,
int &  n1,
int &  n2,
int &  n3 
)

Definition at line 31 of file production/dayone-2017/fullsim/ana_day1.C.

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

32 {
33  int n_smc = 0;
34  int n_strk = 0;
35  int n_both = 0;
36  double d = 0.00001;
37 
38  for (int i=0;i<l.GetLength()-1;++i)
39  {
40  for (int j=i+1;j<l.GetLength();++j)
41  {
42  TLorentzVector dl = l[i]->P4() - l[j]->P4();
43 
44  bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth());
45  bool chktrk = (fabs(dl.X())<d) && (fabs(dl.Y())<d) && (fabs(dl.Z())<d) && (fabs(dl.E())<d);
46  if (chkmc) n_smc++;
47  if (chktrk) n_strk++;
48  if (chktrk && chkmc) n_both++;
49  }
50  }
51  n1 = n_strk;
52  n2 = n_smc;
53  n3 = n_both;
54 }
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 24 of file production/dayone-2017/fullsim/ana_day1.C.

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

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

Definition at line 8 of file production/dayone-2017/fullsim/ana_day1.C.

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

9 {
10  int removed = 0;
11 
12  for (int ii=l.GetLength()-1;ii>=0;--ii)
13  {
14  if ( !(ana->McTruthMatch(l[ii])) )
15  {
16  l.Remove(l[ii]);
17  removed++;
18  }
19  }
20 
21  return removed;
22 }
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 Wed Apr 3 2019 10:03:10 for FairRoot/PandaRoot by   doxygen 1.8.5