FairRoot/PandaRoot
thailand2017/solution/tut_ana_ntp.C
Go to the documentation of this file.
1 class RhoCandList;
2 class RhoCandidate;
3 class PndAnaPidSelector;
4 class PndAnaPidCombiner;
5 class PndAnalysis;
6 class RhoTuple;
7 
8 // **** some auxilliary functions in auxtut.C ****
9 // - FairRunAna* initrun(TString prefix, TString outfile, int min=-1, int max=-1) --> Init FairRunAna
10 // - plotmyhistos() --> Plots all histograms in current TDirectory on a autosized canvas
11 // - writemyhistos() --> Writes all histos in current TFile
12 // - fillM(RhoCandList l, TH1* h) --> Fill mass histogram h with masses of candidates in l
13 // - RemoveGeoManager() --> Temporary fix for error on macro exit
14 // **** some auxilliary functions in auxtut.C ****
15 #include "auxtut.C"
16 
17 void tut_ana_ntp(int nevts = 0, TString prefix = "signal")
18 {
19  // *** some variables
20  int i=0,j=0, k=0, l=0;
21  gStyle->SetOptFit(1011);
22 
23  // *** Initialize FairRunAna with defaults
24  TString OutFile="out_dummy.root";
25  FairRunAna* fRun = initrun(prefix, OutFile);
26  fRun->Init();
27 
28  // *** create an output file for all histograms
29  TFile *out = TFile::Open(prefix+"_ana_ntp.root","RECREATE");
30 
31  // *** create ntuples for J/psi and psi(2S)
32  RhoTuple *njpsi = new RhoTuple("njpsi","J/psi Analysis");
33  RhoTuple *npsip = new RhoTuple("npsip","psi' Analysis");
34 
35  // *** create some columns which might not be filled sometimes
36  njpsi->Column("tjpsim", 0.0f, -999.9f);
37  njpsi->Column("tjpsip", 0.0f, -999.9f);
38  njpsi->Column("tjpsitht", 0.0f, -999.9f);
39 
40  npsip->Column("tpsim", 0.0f, -999.9f);
41  npsip->Column("tpsip", 0.0f, -999.9f);
42  npsip->Column("tpsitht", 0.0f, -999.9f);
43 
44  //
45  // Now the analysis stuff comes...
46  //
47 
48  // *** the data reader object
49  PndAnalysis* theAnalysis = new PndAnalysis();
50  if (nevts==0) nevts= theAnalysis->GetEntries();
51 
52  // *** RhoCandLists for the analysis
53  RhoCandList muplus, muminus, piplus, piminus, jpsi, psi2s, allpart;
54 
55  // *** Mass selector for the jpsi cands
56  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
57  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
58 
59  // *** Pid Selection Algorithms
60  TString pidSelection = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts";
61 
62  // *** the lorentz vector of the initial psi(2S)
63  TLorentzVector ini(0, 0, 6.231552, 7.240065);
64 
65  // *** the QA tool for easy ntuple output
66  PndRhoTupleQA qa(theAnalysis,ini.P());
67 
68 
69  // ***
70  // the event loop
71  // ***
72  while (theAnalysis->GetEvent() && i++<nevts)
73  {
74  if ((i%100)==0) cout<<"evt " << i << endl;
75 
76  // *** Select with no PID info ('All'); type and mass are set
77  theAnalysis->FillList(muplus, "MuonAllPlus", pidSelection);
78  theAnalysis->FillList(muminus, "MuonAllMinus", pidSelection);
79  theAnalysis->FillList(piplus, "PionAllPlus", pidSelection);
80  theAnalysis->FillList(piminus, "PionAllMinus", pidSelection);
81 
82  // *** prepare PndEventShape
83  theAnalysis->FillList(allpart, "All", pidSelection);
84  PndEventShape evtsh(allpart, ini, 0.05, 0.1);
85 
86  // *** combinatorics for J/psi -> mu+ mu-
87  jpsi.Combine(muplus, muminus);
88  jpsi.SetType(443);
89 
90  // *** some mass pre selection
91  //jpsi.Select(jpsiMassSel);
92 
93  // ***
94  // *** do all kind of analysis and store in N-tuple
95  // ***
96  for (j=0;j<jpsi.GetLength();++j)
97  {
98  // get daughters
99  RhoCandidate *mup = jpsi[j]->Daughter(0);
100  RhoCandidate *mum = jpsi[j]->Daughter(1);
101  PndPidCandidate *mup_rec = (PndPidCandidate*)mup->GetRecoCandidate();
102  PndPidCandidate *mum_rec = (PndPidCandidate*)mum->GetRecoCandidate();
103 
104  // get truth information
105  bool mct = theAnalysis->McTruthMatch(jpsi[j]);
106  RhoCandidate *true_jpsi = jpsi[j]->GetMcTruth();
107 
108  // perform vertex fitter
109  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
110  vtxfitter.Fit();
111 
112  RhoCandidate *fitvtx_jpsi = jpsi[j]->GetFit();
113  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
114  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
115  TVector3 vtxpos(-999.,-999.,-999.);
116  if (fitvtx_jpsi) vtxpos = fitvtx_jpsi->Daughter(0)->Pos();
117 
118  // perform mass fit
119  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
120  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
121  mfitter.Fit(); // do fit
122 
123  RhoCandidate *fitmass_jpsi = jpsi[j]->GetFit();
124  double chi2_mass = mfitter.GetChi2(); // get chi2 of fit
125  double prob_mass = mfitter.GetProb(); // access probability of fit
126 
127  // *** now write ntuple information
128 
129  // *** general event info
130  njpsi->Column("ev", (Float_t) i, -999.9f);
131  njpsi->Column("cand", (Float_t) j, -999.9f);
132 
133  // *** basic J/psi info
134  njpsi->Column("jpsim", (Float_t) jpsi[j]->M(), -999.9f);
135  njpsi->Column("jpsip", (Float_t) jpsi[j]->P(), -999.9f);
136  njpsi->Column("jpsipt", (Float_t) jpsi[j]->P3().Pt(), -999.9f);
137  njpsi->Column("jpsitht", (Float_t) jpsi[j]->P3().Theta(), -999.9f);
138  njpsi->Column("jpsimissm", (Float_t) (ini-(jpsi[j]->P4())).M(), -999.9f);
139 
140  // *** MC truth info
141  njpsi->Column("mct", (Float_t) mct, -999.9f);
142  if (true_jpsi)
143  {
144  njpsi->Column("tjpsim", (Float_t) true_jpsi->M(), -999.9f);
145  njpsi->Column("tjpsip", (Float_t) true_jpsi->M(), -999.9f);
146  njpsi->Column("tjpsitht",(Float_t) true_jpsi->P3().Theta(), -999.9f);
147  }
148 
149  // *** fitting info
150  njpsi->Column("jpsimvtx", (Float_t) fitvtx_jpsi->M(), -999.9f);
151  njpsi->Column("chi2vtx", (Float_t) chi2_vtx, -999.9f);
152  njpsi->Column("probvtx", (Float_t) prob_vtx, -999.9f);
153  njpsi->Column("vtxx", (Float_t) vtxpos.X(), -999.9f);
154  njpsi->Column("vtxy", (Float_t) vtxpos.Y(), -999.9f);
155  njpsi->Column("vtxz", (Float_t) vtxpos.Z(), -999.9f);
156 
157  njpsi->Column("jpsimmass", (Float_t) fitmass_jpsi->M(), -999.9f);
158  njpsi->Column("chi2mass", (Float_t) chi2_mass, -999.9f);
159  njpsi->Column("probmass", (Float_t) prob_mass, -999.9f);
160 
161  // *** kinematic info of daughters
162  njpsi->Column("mupp", (Float_t) mup->P(), -999.9f);
163  njpsi->Column("muppt", (Float_t) mup->P3().Pt(), -999.9f);
164  njpsi->Column("muptht", (Float_t) mup->P3().Theta(), -999.9f);
165 
166  njpsi->Column("mump", (Float_t) mum->P(), -999.9f);
167  njpsi->Column("mumpt", (Float_t) mum->P3().Pt(), -999.9f);
168  njpsi->Column("mumtht", (Float_t) mum->P3().Theta(), -999.9f);
169 
170  // *** PID info of daughters
171  njpsi->Column("muppid", (Float_t) mup->GetPidInfo(1), -999.9f);
172  njpsi->Column("mumpid", (Float_t) mum->GetPidInfo(1), -999.9f);
173 
174  // *** and finally FILL Ntuple
175  njpsi->DumpData();
176 
177  }
178 
179  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
180  psi2s.Combine(jpsi, piplus, piminus);
181  psi2s.SetType(88888);
182 
183  for (j=0;j<psi2s.GetLength();++j)
184  {
185  // get daughters
186  RhoCandidate *jp = psi2s[j]->Daughter(0);
187  RhoCandidate *pip = psi2s[j]->Daughter(1);
188  RhoCandidate *pim = psi2s[j]->Daughter(2);
189  RhoCandidate *mup = psi2s[j]->Daughter(0)->Daughter(0);
190  RhoCandidate *mum = psi2s[j]->Daughter(0)->Daughter(1);
191 
192  PndPidCandidate *pip_rec = (PndPidCandidate*)pip->GetRecoCandidate();
193  PndPidCandidate *pim_rec = (PndPidCandidate*)pim->GetRecoCandidate();
194 
195  // get truth information
196  bool mct = theAnalysis->McTruthMatch(psi2s[j]);
197  RhoCandidate *true_psi = psi2s[j]->GetMcTruth();
198 
199  // do 4C fit
200  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
201  fitter.Add4MomConstraint(ini); // set 4 constraint
202  fitter.Fit(); // do fit
203  RhoCandidate *fit4c_jpsi = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
204 
205  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
206  double prob_4c = fitter.GetProb(); // access probability of fit
207 
208  // *** general event info
209  npsip->Column("ev", (Float_t) i, -999.9f);
210  npsip->Column("cand", (Float_t) j, -999.9f);
211 
212  // *** basic psi(2s) info
213  qa.qaCand("psi", psi2s[j], npsip);
214 
215  // *** basic J/psi info
216  qa.qaCand("jpsi" , jp , npsip);
217  qa.qaCand("jpsi4cfit", fit4c_jpsi, npsip);
218 
219  // *** MC truth info
220  npsip->Column("mct", (Float_t) mct, -999.9f);
221  if (true_psi)
222  {
223  npsip->Column("tpsim", (Float_t) true_psi->M(), -999.9f);
224  npsip->Column("tpsip", (Float_t) true_psi->M(), -999.9f);
225  npsip->Column("tpsitht",(Float_t) true_psi->P3().Theta(), -999.9f);
226  }
227 
228  // *** kinematic info of daughters
229  qa.qaCand("pip", pip, npsip);
230  qa.qaCand("pim", pip, npsip);
231 
232  // *** PID info of daughters
233  qa.qaPid("pip", pip, npsip);
234  qa.qaPid("pim", pip, npsip);
235  qa.qaPid("mup", mup, npsip);
236  qa.qaPid("mum", mum, npsip);
237 
238  // *** Event shape info
239  qa.qaEventShapeShort("evtsh", &evtsh, npsip);
240 
241  // *** and finally FILL Ntuple
242  npsip->DumpData();
243  }
244 
245  }
246 
247  // *** write out all the histos
248  out->cd();
249 
250  njpsi->GetInternalTree()->Write();
251  npsip->GetInternalTree()->Write();
252 
253  out->Save();
254 
255  // *** temporaty fix to avoid error on macro exit
256  RemoveGeoManager();
257 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
initrun
FairRunAna * initrun(TString prefix, TString outfile, int min=-1, int max=-1)
Definition: QA/auxi.C:32
RhoKinFitter::AddMassConstraint
void AddMassConstraint(double mass)
Definition: RhoKinFitter.cxx:56
i
Int_t i
Definition: run_full.C:25
RemoveGeoManager
void RemoveGeoManager()
Definition: auxtut.C:11
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
RhoCandidate::GetRecoCandidate
PndPidCandidate * GetRecoCandidate() const
Definition: RhoCandidate.h:376
RhoTuple::RhoTuple
RhoTuple()
Definition: RhoTuple.cxx:27
RhoCandidate::Daughter
RhoCandidate * Daughter(Int_t n)
Definition: RhoCandidate.cxx:991
piplus
CandList piplus
Definition: softtrigger_kin5.C:144
RhoKinFitter::Add4MomConstraint
void Add4MomConstraint(TLorentzVector lv)
Definition: RhoKinFitter.cxx:34
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
f
TFile * f
Definition: bump_analys.C:12
RhoTuple::Column
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
Definition: RhoTuple.cxx:56
RhoFitterBase::GetProb
double GetProb() const
Definition: RhoFitterBase.h:50
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
RhoCandidate::P
Double_t P() const
Definition: RhoCandidate.cxx:1629
RhoTuple::DumpData
void DumpData()
Definition: RhoTuple.cxx:391
P
GeV c P
Definition: ReadMesDecCorrRes_PiT.C:42
tut_ana_ntp
void tut_ana_ntp(int nevts=0, TString prefix="signal")
Definition: rho/tut_ana_ntp.C:8
piminus
CandList piminus
Definition: softtrigger_kin5.C:144
RhoFitterBase::GetChi2
double GetChi2() const
Definition: RhoFitterBase.h:48
RhoTuple::GetInternalTree
TTree * GetInternalTree()
Definition: RhoTuple.h:207
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
RhoKinFitter::Fit
Bool_t Fit()
Definition: RhoKinFitter.cxx:62
muminus
CandList muminus
Definition: softtrigger_kin5.C:144
mct
PndMCTrack * mct
Definition: hist-t7.C:147
mup
static const double mup
Definition: mzparameters.h:21
RhoCandidate::GetPidInfo
double GetPidInfo(int hypo)
Definition: RhoCandidate.cxx:1752
RhoCandidate::P3
TVector3 P3() const
Definition: RhoCandidate.h:199
Generated on Wed Apr 3 2019 10:03:09 for FairRoot/PandaRoot by   doxygen 1.8.5