FairRoot/PandaRoot
rho/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 void tut_ana_ntp(int nevts = 0, TString prefix = "signal")
9 {
10  // *** some variables
11  int i=0,j=0, k=0, l=0;
12  gStyle->SetOptFit(1011);
13 
14  // *** the output file for FairRunAna
15  TString OutFile="out_dummy.root";
16 
17  // *** the files coming from the simulation
18  TString inPidFile = prefix+"_pid.root"; // this file contains the PndPidCandidates and McTruth
19  TString inParFile = prefix+"_par.root";
20 
21  // *** PID table with selection thresholds; can be modified by the user
22  TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";
23 
24  // *** initialization
25  FairLogger::GetLogger()->SetLogToFile(kFALSE);
26  FairRunAna* fRun = new FairRunAna();
27  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
28  fRun->SetSource(new FairFileSource(inPidFile));
29 
30  // *** setup parameter database
31  FairParRootFileIo* parIO = new FairParRootFileIo();
32  parIO->open(inParFile);
33  FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo();
34  parIOPid->open(pidParFile.Data(),"in");
35 
36  rtdb->setFirstInput(parIO);
37  rtdb->setSecondInput(parIOPid);
38  rtdb->setOutput(parIO);
39 
40  fRun->SetOutputFile(OutFile);
41  fRun->Init();
42 
43  // *** create an output file for all histograms
44  TFile *out = TFile::Open(prefix+"_ana_ntp.root","RECREATE");
45 
46  // *** create ntuples for J/psi and psi(2S)
47  RhoTuple *njpsi = new RhoTuple("njpsi","J/psi Analysis");
48  RhoTuple *npsip = new RhoTuple("npsip","psi' Analysis");
49 
50  // *** create some columns which might not be filled sometimes
51  njpsi->Column("tjpsim", 0.0f, -999.9f);
52  njpsi->Column("tjpsip", 0.0f, -999.9f);
53  njpsi->Column("tjpsitht", 0.0f, -999.9f);
54 
55  npsip->Column("tpsim", 0.0f, -999.9f);
56  npsip->Column("tpsip", 0.0f, -999.9f);
57  npsip->Column("tpsitht", 0.0f, -999.9f);
58 
59  //
60  // Now the analysis stuff comes...
61  //
62 
63 
64  // *** the data reader object
65  PndAnalysis* theAnalysis = new PndAnalysis();
66  if (nevts==0) nevts= theAnalysis->GetEntries();
67 
68  // *** RhoCandLists for the analysis
69  RhoCandList muplus, muminus, piplus, piminus, jpsi, psi2s;
70 
71  // *** Mass selector for the jpsi cands
72  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
73  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
74 
75  // *** Pid Selection Algorithms
76  TString pidSelection = "PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc;PidAlgoStt;PidAlgoMdtHardCuts";
77 
78  // *** the lorentz vector of the initial psi(2S)
79  TLorentzVector ini(0, 0, 6.231552, 7.240065);
80 
81  // ***
82  // the event loop
83  // ***
84  while (theAnalysis->GetEvent() && i++<nevts)
85  {
86  if ((i%100)==0) cout<<"evt " << i << endl;
87 
88  // *** Select with no PID info ('All'); type and mass are set
89  theAnalysis->FillList(muplus, "MuonAllPlus", pidSelection);
90  theAnalysis->FillList(muminus, "MuonAllMinus", pidSelection);
91  theAnalysis->FillList(piplus, "PionAllPlus", pidSelection);
92  theAnalysis->FillList(piminus, "PionAllMinus", pidSelection);
93 
94  // *** combinatorics for J/psi -> mu+ mu-
95  jpsi.Combine(muplus, muminus);
96  jpsi.SetType(443);
97 
98  // *** some mass pre selection
99  //jpsi.Select(jpsiMassSel);
100 
101  // ***
102  // *** do all kind of analysis and store in N-tuple
103  // ***
104  for (j=0;j<jpsi.GetLength();++j)
105  {
106  // get daughters
107  RhoCandidate *mup = jpsi[j]->Daughter(0);
108  RhoCandidate *mum = jpsi[j]->Daughter(1);
109  PndPidCandidate *mup_rec = (PndPidCandidate*)mup->GetRecoCandidate();
110  PndPidCandidate *mum_rec = (PndPidCandidate*)mum->GetRecoCandidate();
111 
112  // get truth information
113  bool mct = theAnalysis->McTruthMatch(jpsi[j]);
114  RhoCandidate *true_jpsi = jpsi[j]->GetMcTruth();
115 
116  // perform vertex fitter
117  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
118  vtxfitter.Fit();
119 
120  RhoCandidate *fitvtx_jpsi = jpsi[j]->GetFit();
121  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
122  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
123  TVector3 vtxpos(-999.,-999.,-999.);
124  if (fitvtx_jpsi) vtxpos = fitvtx_jpsi->Daughter(0)->Pos();
125 
126  // perform mass fit
127  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
128  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
129  mfitter.Fit(); // do fit
130 
131  RhoCandidate *fitmass_jpsi = jpsi[j]->GetFit();
132  double chi2_mass = mfitter.GetChi2(); // get chi2 of fit
133  double prob_mass = mfitter.GetProb(); // access probability of fit
134 
135  // *** now write ntuple information
136 
137  // *** general event info
138  njpsi->Column("ev", (Float_t) i, -999.9f);
139  njpsi->Column("cand", (Float_t) j, -999.9f);
140 
141  // *** basic J/psi info
142  njpsi->Column("jpsim", (Float_t) jpsi[j]->M(), -999.9f);
143  njpsi->Column("jpsip", (Float_t) jpsi[j]->P(), -999.9f);
144  njpsi->Column("jpsipt", (Float_t) jpsi[j]->P3().Pt(), -999.9f);
145  njpsi->Column("jpsitht", (Float_t) jpsi[j]->P3().Theta(), -999.9f);
146  njpsi->Column("jpsimissm", (Float_t) (ini-(jpsi[j]->P4())).M(), -999.9f);
147 
148  // *** MC truth info
149  njpsi->Column("mct", (Float_t) mct, -999.9f);
150  if (true_jpsi)
151  {
152  njpsi->Column("tjpsim", (Float_t) true_jpsi->M(), -999.9f);
153  njpsi->Column("tjpsip", (Float_t) true_jpsi->M(), -999.9f);
154  njpsi->Column("tjpsitht",(Float_t) true_jpsi->P3().Theta(), -999.9f);
155  }
156 
157  // *** fitting info
158  njpsi->Column("jpsimvtx", (Float_t) fitvtx_jpsi->M(), -999.9f);
159  njpsi->Column("chi2vtx", (Float_t) chi2_vtx, -999.9f);
160  njpsi->Column("probvtx", (Float_t) prob_vtx, -999.9f);
161  njpsi->Column("vtxx", (Float_t) vtxpos.X(), -999.9f);
162  njpsi->Column("vtxy", (Float_t) vtxpos.Y(), -999.9f);
163  njpsi->Column("vtxz", (Float_t) vtxpos.Z(), -999.9f);
164 
165  njpsi->Column("jpsimmass", (Float_t) fitmass_jpsi->M(), -999.9f);
166  njpsi->Column("chi2mass", (Float_t) chi2_mass, -999.9f);
167  njpsi->Column("probmass", (Float_t) prob_mass, -999.9f);
168 
169  // *** kinematic info of daughters
170  njpsi->Column("mupp", (Float_t) mup->P(), -999.9f);
171  njpsi->Column("muppt", (Float_t) mup->P3().Pt(), -999.9f);
172  njpsi->Column("muptht", (Float_t) mup->P3().Theta(), -999.9f);
173 
174  njpsi->Column("mump", (Float_t) mum->P(), -999.9f);
175  njpsi->Column("mumpt", (Float_t) mum->P3().Pt(), -999.9f);
176  njpsi->Column("mumtht", (Float_t) mum->P3().Theta(), -999.9f);
177 
178  // *** PID info of daughters
179  njpsi->Column("muppid", (Float_t) mup->GetPidInfo(1), -999.9f);
180  njpsi->Column("mumpid", (Float_t) mum->GetPidInfo(1), -999.9f);
181 
182  // *** and finally FILL Ntuple
183  njpsi->DumpData();
184 
185  }
186 
187  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
188  psi2s.Combine(jpsi, piplus, piminus);
189  psi2s.SetType(88888);
190 
191  for (j=0;j<psi2s.GetLength();++j)
192  {
193  // get daughters
194  RhoCandidate *jp = psi2s[j]->Daughter(0);
195  RhoCandidate *pip = psi2s[j]->Daughter(1);
196  RhoCandidate *pim = psi2s[j]->Daughter(2);
197 
198  PndPidCandidate *pip_rec = (PndPidCandidate*)pip->GetRecoCandidate();
199  PndPidCandidate *pim_rec = (PndPidCandidate*)pim->GetRecoCandidate();
200 
201  // get truth information
202  bool mct = theAnalysis->McTruthMatch(psi2s[j]);
203  RhoCandidate *true_psi = psi2s[j]->GetMcTruth();
204 
205  // do 4C fit
206  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
207  fitter.Add4MomConstraint(ini); // set 4 constraint
208  fitter.Fit(); // do fit
209  RhoCandidate *fit4c_jpsi = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
210 
211  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
212  double prob_4c = fitter.GetProb(); // access probability of fit
213 
214  // *** general event info
215  npsip->Column("ev", (Float_t) i, -999.9f);
216  npsip->Column("cand", (Float_t) j, -999.9f);
217 
218  // *** basic psi(2s) info
219  npsip->Column("psim", (Float_t) psi2s[j]->M(), -999.9f);
220  npsip->Column("psip", (Float_t) psi2s[j]->P(), -999.9f);
221  npsip->Column("psipt", (Float_t) psi2s[j]->P3().Pt(), -999.9f);
222  npsip->Column("psitht", (Float_t) psi2s[j]->P3().Theta(), -999.9f);
223 
224  // *** basic J/psi info
225  npsip->Column("jpsim", (Float_t) jp->M(), -999.9f);
226  npsip->Column("jpsip", (Float_t) jp->P(), -999.9f);
227  npsip->Column("jpsipt", (Float_t) jp->P3().Pt(), -999.9f);
228  npsip->Column("jpsitht",(Float_t) jp->P3().Theta(), -999.9f);
229 
230  npsip->Column("jpsim4c",(Float_t) fit4c_jpsi->M(), -999.9f);
231 
232  // *** MC truth info
233  npsip->Column("mct", (Float_t) mct, -999.9f);
234  if (true_psi)
235  {
236  npsip->Column("tpsim", (Float_t) true_psi->M(), -999.9f);
237  npsip->Column("tpsip", (Float_t) true_psi->M(), -999.9f);
238  npsip->Column("tpsitht",(Float_t) true_psi->P3().Theta(), -999.9f);
239  }
240 
241  // *** kinematic info of daughters
242  npsip->Column("pipp", (Float_t) pip->P(), -999.9f);
243  npsip->Column("pippt", (Float_t) pip->P3().Pt(), -999.9f);
244  npsip->Column("piptht", (Float_t) pip->P3().Theta(), -999.9f);
245 
246  npsip->Column("pimp", (Float_t) pim->P(), -999.9f);
247  npsip->Column("pimpt", (Float_t) pim->P3().Pt(), -999.9f);
248  npsip->Column("pimtht", (Float_t) pim->P3().Theta(), -999.9f);
249 
250  // *** PID info of daughters
251  npsip->Column("pippid", (Float_t) pip->GetPidInfo(2), -999.9f);
252  npsip->Column("pimpid", (Float_t) pim->GetPidInfo(2), -999.9f);
253 
254  // *** and finally FILL Ntuple
255  npsip->DumpData();
256  }
257 
258  }
259 
260  // *** write out all the histos
261  out->cd();
262 
263  njpsi->GetInternalTree()->Write();
264  npsip->GetInternalTree()->Write();
265 
266  out->Save();
267 
268 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
RhoKinFitter::AddMassConstraint
void AddMassConstraint(double mass)
Definition: RhoKinFitter.cxx:56
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
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
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
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 08:44:14 for FairRoot/PandaRoot by   doxygen 1.8.5