FairRoot/PandaRoot
Functions
rho/tut_ana.C File Reference

Go to the source code of this file.

Functions

int SelectTruePid (PndAnalysis *ana, RhoCandList &l)
 
void tut_ana (int nevts=0, TString prefix="signal")
 

Function Documentation

int SelectTruePid ( PndAnalysis *  ana,
RhoCandList &  l 
)

Definition at line 8 of file rho/tut_ana.C.

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

Referenced by tut_ana().

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
void tut_ana ( int  nevts = 0,
TString  prefix = "signal" 
)

Definition at line 25 of file rho/tut_ana.C.

References RhoKinFitter::Add4MomConstraint(), RhoKinFitter::AddMassConstraint(), RhoCandList::Combine(), 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, piminus, piplus, PndAnalysis::PndAnalysis(), RhoCandidate::Pos(), rtdb, RhoCandList::Select(), SelectTruePid(), RhoCandList::SetType(), and TString.

26 {
27  // *** some variables
28  int i=0,j=0, k=0, l=0;
29  gStyle->SetOptFit(1011);
30 
31  // *** the output file for FairRunAna
32  TString OutFile="out_dummy.root";
33 
34  // *** the files coming from the simulation
35  TString inPidFile = prefix+"_pid.root"; // this file contains the PndPidCandidates and McTruth
36  TString inParFile = prefix+"_par.root";
37 
38  // *** PID table with selection thresholds; can be modified by the user
39  TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";
40 
41  // *** initialization
42  //FairLogger::GetLogger()->SetLogToFile(kFALSE);
43  FairRunAna* fRun = new FairRunAna();
44  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
45  fRun->SetSource(new FairFileSource(inPidFile));
46 
47  // *** setup parameter database
48  FairParRootFileIo* parIO = new FairParRootFileIo();
49  parIO->open(inParFile);
50  FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo();
51  parIOPid->open(pidParFile.Data(),"in");
52 
53  rtdb->setFirstInput(parIO);
54  rtdb->setSecondInput(parIOPid);
55  rtdb->setOutput(parIO);
56 
57  fRun->SetOutputFile(OutFile);
58  fRun->Init();
59 
60  // *** create an output file for all histograms
61  TFile *out = TFile::Open(prefix+"_ana.root","RECREATE");
62 
63  // *** create some histograms
64  TH1F *hjpsim_all = new TH1F("hjpsim_all","J/#psi mass (all)",200,0,4.5);
65  TH1F *hpsim_all = new TH1F("hpsim_all","#psi(2S) mass (all)",200,0,5);
66 
67  TH1F *hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4.5);
68  TH1F *hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5);
69 
70  TH1F *hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (tight pid)",200,0,4.5);
71  TH1F *hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (tight pid)",200,0,5);
72 
73  TH1F *hjpsim_trpid = new TH1F("hjpsim_trpid","J/#psi mass (true pid)",200,0,4.5);
74  TH1F *hpsim_trpid = new TH1F("hpsim_trpid","#psi(2S) mass (true pid)",200,0,5);
75 
76 
77  TH1F *hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4.5);
78  TH1F *hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass (full truth match)",200,0,5);
79 
80  TH1F *hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4.5);
81  TH1F *hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5);
82 
83  TH1F *hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2);
84  TH1F *hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2);
85 
86 
87  TH1F *hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass (vertex fit)",200,0,4.5);
88  TH1F *hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4.5);
89  TH1F *hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (mass constraint fit)",200,0,4.5);
90 
91  TH1F *hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf", "J/#psi: #chi^{2} vertex fit",100,0,10);
92  TH1F *hpsi_chi2_4c = new TH1F("hpsi_chi2_4c", "#psi(2S): #chi^{2} 4C fit",100,0,250);
93  TH1F *hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf", "J/#psi: #chi^{2} mass fit",100,0,10);
94 
95  TH1F *hjpsi_prob_vf = new TH1F("hjpsi_prob_vf", "J/#psi: Prob vertex fit",100,0,1);
96  TH1F *hpsi_prob_4c = new TH1F("hpsi_prob_4c", "#psi(2S): Prob 4C fit",100,0,1);
97  TH1F *hjpsi_prob_mf = new TH1F("hjpsi_prob_mf", "J/#psi: Prob mass fit",100,0,1);
98 
99  TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2);
100 
101  //
102  // Now the analysis stuff comes...
103  //
104 
105 
106  // *** the data reader object
107  PndAnalysis* theAnalysis = new PndAnalysis();
108  if (nevts==0) nevts= theAnalysis->GetEntries();
109 
110  // *** RhoCandLists for the analysis
111  RhoCandList muplus, muminus, piplus, piminus, jpsi, psi2s;
112 
113  // *** Mass selector for the jpsi cands
114  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
115  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
116 
117  // *** the lorentz vector of the initial psi(2S)
118  TLorentzVector ini(0, 0, 6.231552, 7.240065);
119 
120  // ***
121  // the event loop
122  // ***
123  while (theAnalysis->GetEvent() && i++<nevts)
124  {
125  if ((i%100)==0) cout<<"evt " << i << endl;
126 
127  // *** Select with no PID info ('All'); type and mass are set
128  theAnalysis->FillList(muplus, "MuonAllPlus");
129  theAnalysis->FillList(muminus, "MuonAllMinus");
130  theAnalysis->FillList(piplus, "PionAllPlus");
131  theAnalysis->FillList(piminus, "PionAllMinus");
132 
133  // *** combinatorics for J/psi -> mu+ mu-
134  jpsi.Combine(muplus, muminus);
135 
136 
137  // ***
138  // *** do the TRUTH MATCH for jpsi
139  // ***
140  jpsi.SetType(443);
141 
142  for (j=0;j<jpsi.GetLength();++j)
143  {
144  hjpsim_all->Fill( jpsi[j]->M() );
145 
146  if (theAnalysis->McTruthMatch(jpsi[j]))
147  {
148  hjpsim_ftm->Fill( jpsi[j]->M() );
149  hjpsim_diff->Fill( jpsi[j]->GetMcTruth()->M() - jpsi[j]->M() );
150  }
151  else
152  hjpsim_nm->Fill( jpsi[j]->M() );
153  }
154 
155  // ***
156  // *** do VERTEX FIT (J/psi)
157  // ***
158  for (j=0;j<jpsi.GetLength();++j)
159  {
160  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
161  vtxfitter.Fit();
162 
163  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
164  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
165  hjpsi_chi2_vf->Fill(chi2_vtx);
166  hjpsi_prob_vf->Fill(prob_vtx);
167 
168  if ( prob_vtx > 0.01 ) // when good enough, fill some histos
169  {
170  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
171  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
172 
173  hjpsim_vf->Fill(jfit->M());
174  hvpos->Fill(jVtx.X(),jVtx.Y());
175  }
176  }
177 
178  // *** some rough mass selection
179  jpsi.Select(jpsiMassSel);
180 
181  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
182  psi2s.Combine(jpsi, piplus, piminus);
183 
184 
185  // ***
186  // *** do the TRUTH MATCH for psi(2S)
187  // ***
188  psi2s.SetType(88888);
189 
190  for (j=0;j<psi2s.GetLength();++j)
191  {
192  hpsim_all->Fill( psi2s[j]->M() );
193 
194  if (theAnalysis->McTruthMatch(psi2s[j]))
195  {
196  hpsim_ftm->Fill( psi2s[j]->M() );
197  hpsim_diff->Fill( psi2s[j]->GetMcTruth()->M() - psi2s[j]->M() );
198  }
199  else
200  hpsim_nm->Fill( psi2s[j]->M() );
201  }
202 
203 
204  // ***
205  // *** do 4C FIT (initial psi(2S) system)
206  // ***
207  for (j=0;j<psi2s.GetLength();++j)
208  {
209  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
210  fitter.Add4MomConstraint(ini); // set 4 constraint
211  fitter.Fit(); // do fit
212 
213  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
214  double prob_4c = fitter.GetProb(); // access probability of fit
215  hpsi_chi2_4c->Fill(chi2_4c);
216  hpsi_prob_4c->Fill(prob_4c);
217 
218  if ( prob_4c > 0.01 ) // when good enough, fill some histo
219  {
220  RhoCandidate *jfit = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
221 
222  hjpsim_4cf->Fill(jfit->M());
223  }
224  }
225 
226 
227  // ***
228  // *** do MASS CONSTRAINT FIT (J/psi)
229  // ***
230  for (j=0;j<jpsi.GetLength();++j)
231  {
232  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
233  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
234  mfitter.Fit(); // do fit
235 
236  double chi2_m = mfitter.GetChi2(); // get chi2 of fit
237  double prob_m = mfitter.GetProb(); // access probability of fit
238  hjpsi_chi2_mf->Fill(chi2_m);
239  hjpsi_prob_mf->Fill(prob_m);
240 
241  if ( prob_m > 0.01 ) // when good enough, fill some histo
242  {
243  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
244  hjpsim_mcf->Fill(jfit->M());
245  }
246  }
247 
248 
249  // ***
250  // *** TRUE PID combinatorics
251  // ***
252 
253  // *** do MC truth match for PID type
254  SelectTruePid(theAnalysis, muplus);
255  SelectTruePid(theAnalysis, muminus);
256  SelectTruePid(theAnalysis, piplus);
257  SelectTruePid(theAnalysis, piminus);
258 
259  // *** all combinatorics again with true PID
260  jpsi.Combine(muplus, muminus);
261  for (j=0;j<jpsi.GetLength();++j) hjpsim_trpid->Fill( jpsi[j]->M() );
262  jpsi.Select(jpsiMassSel);
263 
264  psi2s.Combine(jpsi, piplus, piminus);
265  for (j=0;j<psi2s.GetLength();++j) hpsim_trpid->Fill( psi2s[j]->M() );
266 
267 
268  // ***
269  // *** LOOSE PID combinatorics
270  // ***
271 
272  // *** and again with PidAlgoMvd;PidAlgoStt;PidAlgoDrc and loose selection
273  theAnalysis->FillList(muplus, "MuonLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
274  theAnalysis->FillList(muminus, "MuonLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
275  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
276  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
277 
278  jpsi.Combine(muplus, muminus);
279  for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j]->M() );
280  jpsi.Select(jpsiMassSel);
281 
282  psi2s.Combine(jpsi, piplus, piminus);
283  for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j]->M() );
284 
285 
286  // ***
287  // *** TIGHT PID combinatorics
288  // ***
289 
290  // *** and again with PidAlgoMvd;PidAlgoStt and tight selection
291  theAnalysis->FillList(muplus, "MuonTightPlus", "PidAlgoMdtHardCuts");
292  theAnalysis->FillList(muminus, "MuonTightMinus", "PidAlgoMdtHardCuts");
293  theAnalysis->FillList(piplus, "PionLoosePlus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
294  theAnalysis->FillList(piminus, "PionLooseMinus", "PidAlgoMvd;PidAlgoStt;PidAlgoDrc");
295 
296  jpsi.Combine(muplus, muminus);
297  for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j]->M() );
298  jpsi.Select(jpsiMassSel);
299 
300  psi2s.Combine(jpsi, piplus, piminus);
301  for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j]->M() );
302 
303  }
304 
305  // *** write out all the histos
306  out->cd();
307 
308  hjpsim_all->Write();
309  hpsim_all->Write();
310  hjpsim_lpid->Write();
311  hpsim_lpid->Write();
312  hjpsim_tpid->Write();
313  hpsim_tpid->Write();
314  hjpsim_trpid->Write();
315  hpsim_trpid->Write();
316 
317  hjpsim_ftm->Write();
318  hpsim_ftm->Write();
319  hjpsim_nm->Write();
320  hpsim_nm->Write();
321 
322  hpsim_diff->Write();
323  hjpsim_diff->Write();
324 
325  hjpsim_vf->Write();
326  hjpsim_4cf->Write();
327  hjpsim_mcf->Write();
328 
329  hjpsi_chi2_vf->Write();
330  hpsi_chi2_4c->Write();
331  hjpsi_chi2_mf->Write();
332 
333  hjpsi_prob_vf->Write();
334  hpsi_prob_4c->Write();
335  hjpsi_prob_mf->Write();
336 
337  hvpos->Write();
338 
339  out->Save();
340 
341 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
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
SelectTruePid
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
Definition: rho/tut_ana.C:8
piplus
CandList piplus
Definition: softtrigger_kin5.C:144
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
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
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
hvpos
TH2F * hvpos
Definition: plot_eta_c_stt.C:36
Generated on Sun Aug 18 2019 03:02:58 for FairRoot/PandaRoot by   doxygen 1.8.8