FairRoot/PandaRoot
Functions
anatut_psi2s.C File Reference

Go to the source code of this file.

Functions

int SelectPdgCode (TCandList &mct, TCandList &l)
 
void anatut_psi2s (int nevts=0)
 

Function Documentation

void anatut_psi2s ( int  nevts = 0)

Definition at line 28 of file anatut_psi2s.C.

References RhoKinFitter::AddMassConstraint(), all, PndAnalysis::FillList(), RhoKinFitter::Fit(), RhoFitterBase::Fit(), Rho4CFitter::FitConserveMasses(), fRun, RhoFitterBase::GetChi2(), PndAnalysis::GetEntries(), PndAnalysis::GetEvent(), hvpos, i, inRecoFile, inSimFile, muminus, muplus, out, piminus, piplus, rtdb, SelectPdgCode(), and TString.

29 {
30  // *** some variables
31  int i=0,j=0, k=0, l=0;
32 
33  TString OutFile="output.root";
34 
35  // *** the files coming from the simulation
36  TString inPidFile = "pid_sttcombi.root"; // this file contains the PndPidCandidates
37  TString inRecoFile = "reco_sttcombi.root";
38  TString inSimFile = "points_sttcombi.root"; // this file contains the MC truth
39  TString inParFile = "params_sttcombi.root";
40 
41  gStyle->SetOptFit(1011);
42  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
43 
44  FairLogger::GetLogger()->SetLogToFile(kFALSE);
45 
46  // *** initialization
47  FairRunAna* fRun = new FairRunAna();
48  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
49  fRun->SetInputFile(inSimFile);
50  fRun->AddFriend(inPidFile);
51  fRun->AddFriend(inRecoFile);
52 
53  FairParRootFileIo* parIO = new FairParRootFileIo();
54  parIO->open(inParFile);
55  rtdb->setFirstInput(parIO);
56  rtdb->setOutput(parIO);
57 
58  fRun->SetOutputFile(OutFile);
59  fRun->Init();
60 
61  // *** create an output file for all histograms
62  TFile *out = TFile::Open("output_psi2sana.root","RECREATE");
63 
64  // *** create some histograms
65  TH1F *hjpsim_nopid = new TH1F("hjpsim_nopid","J/#psi mass (no pid)",200,0,4);
66  TH1F *hpsim_nopid = new TH1F("hpsim_nopid","#psi(2S) mass (no pid)",200,0,5);
67  TH1F *hjpsim_lpid = new TH1F("hjpsim_lpid","J/#psi mass (loose pid)",200,0,4);
68  TH1F *hpsim_lpid = new TH1F("hpsim_lpid","#psi(2S) mass (loose pid)",200,0,5);
69  TH1F *hjpsim_tpid = new TH1F("hjpsim_tpid","J/#psi mass (true pid)",200,0,4);
70  TH1F *hpsim_tpid = new TH1F("hpsim_tpid","#psi(2S) mass (true pid)",200,0,5);
71 
72  TH1F *hjpsim_ftm = new TH1F("hjpsim_ftm","J/#psi mass (full truth match)",200,0,4);
73  TH1F *hpsim_ftm = new TH1F("hpsim_ftm","#psi(2S) mass ( truth match)",200,0,5);
74  TH1F *hjpsim_nm = new TH1F("hjpsim_nm","J/#psi mass (no truth match)",200,0,4);
75  TH1F *hpsim_nm = new TH1F("hpsim_nm","#psi(2S) mass (no truth match)",200,0,5);
76 
77  TH1F *hjpsim_diff = new TH1F("hjpsim_diff","J/#psi mass diff to truth",100,-2,2);
78  TH1F *hpsim_diff = new TH1F("hpsim_diff","#psi(2S) mass diff to truth",100,-2,2);
79 
80  TH1F *hjpsim_vf = new TH1F("hjpsim_vf","J/#psi mass vertex fit",200,0,4);
81  TH1F *hjpsim_4cf = new TH1F("hjpsim_4cf","J/#psi mass (4C fit)",200,0,4);
82  TH1F *hjpsim_mcf = new TH1F("hjpsim_mcf","J/#psi mass (4C fit)",200,0,4);
83 
84  TH1F *hjpsi_chi2_vf = new TH1F("hjpsi_chi2_vf","J/#psi, #chi^{2} vertex fit",100,0,10);
85  TH1F *hpsi_chi2_4c = new TH1F("hpsi_chi2_4c","#psi, #chi^{2} 4C fit",100,0,250);
86  TH1F *hjpsi_chi2_mf = new TH1F("hjpsi_chi2_mf","J/#psi, #chi^{2} mass fit",100,0,10);
87 
88  TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay vertex",100,-2,2,100,-2,2);
89 
90  // *** the data reader object
91  PndAnalysis* theAnalysis = new PndAnalysis();
92  if (nevts==0) nevts= theAnalysis->GetEntries();
93 
94  // *** TCandLists for the analysis
95  TCandList all, chrg, mctrk, muplus, muminus, piplus, piminus, jpsi, psi2s;
96 
97  // *** Mass selector for the jpsi cands
98  TPidMassSelector *jpsiMassSel=new TPidMassSelector("jpsi",3.096,1.0);
99 
100  // *** the MC truth matcher
101  PndMcTruthMatch mcm;
102 
103  // *** the lorentz vector of the initial psi(2S)
104  TLorentzVector ini(0, 0, 6.231552, 7.240065);
105 
106  // ***
107  // the event loop
108  // ***
109  while (theAnalysis->GetEvent() && i++<nevts)
110  {
111  if ((i%100)==0) cout<<"evt " << i << endl;
112 
113  // *** the MC Truth objects
114  theAnalysis->FillList(mctrk,"McTruth");
115 
116  // *** Select with no PID info ('All'); type and mass are set
117  theAnalysis->FillList(muplus, "MuonAllPlus");
118  theAnalysis->FillList(muminus, "MuonAllMinus");
119  theAnalysis->FillList(piplus, "PionAllPlus");
120  theAnalysis->FillList(piminus, "PionAllMinus");
121 
122  jpsi.Combine(muplus, muminus);
123 
124  // *** do the truth match for jpsi
125  mcm.SetType(jpsi,"J/psi");
126  for (j=0;j<jpsi.GetLength();++j)
127  {
128  hjpsim_nopid->Fill( jpsi[j].M() );
129  if (mcm.MctMatch(jpsi[j], mctrk)) hjpsim_ftm->Fill( jpsi[j].M() );
130  else hjpsim_nm->Fill( jpsi[j].M() );
131  }
132 
133  // *** do vertex fitting (J/psi)
134  for (j=0;j<jpsi.GetLength();++j)
135  {
136  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
137 
138  vtxfitter.Fit();
139  TCandidate *jfit = vtxfitter.FittedCand(jpsi[j]); // access the fitted cand
140  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
141 
142  double chi2_vtx=vtxfitter.GlobalChi2(); // access chi2 of fit
143  hjpsi_chi2_vf->Fill(chi2_vtx);
144 
145  if (chi2_vtx<2) // when good enough, fill some histos
146  {
147  hjpsim_vf->Fill(jfit->M());
148  hvpos->Fill(jVtx.X(),jVtx.Y());
149  }
150  }
151 
152  // *** some rough mass selection
153  //jpsi.Select(jpsiMassSel);
154 
155  psi2s.Combine(jpsi, piplus, piminus);
156 
157  // *** do the truth match for psi2s
158  mcm.SetType(psi2s,"psi(2S)");
159 
160  for (j=0;j<psi2s.GetLength();++j)
161  {
162  hpsim_nopid->Fill( psi2s[j].M() );
163  if (mcm.MctMatch(psi2s[j], mctrk))
164  {
165  hpsim_ftm->Fill( psi2s[j].M() );
166 
167  TCandidate truePsi = mctrk[psi2s[j].GetMcIdx()];
168  TCandidate trueJ = mctrk[psi2s[j].Daughter(0)->GetMcIdx()];
169 
170  hjpsim_diff->Fill(trueJ.M()-psi2s[j].Daughter(0)->M());
171  hpsim_diff->Fill(truePsi.M()-psi2s[j].M());
172  }
173  else hpsim_nm->Fill( psi2s[j].M() );
174  }
175 
176  // *** do 4c fit (initial psi(2S) system)
177  for (j=0;j<psi2s.GetLength();++j)
178  {
179  Rho4CFitter fitter(psi2s[j],ini);
180  fitter.FitConserveMasses();
181 
182  double chi2_4c=fitter.GetChi2();
183  hpsi_chi2_4c->Fill(chi2_4c);
184 
185  TCandidate *jfit = fitter.FittedCand(*(psi2s[j].Daughter(0)));
186 
187  TCandidate *mupfit = fitter.FittedCand(*(psi2s[j].Daughter(0)->Daughter(0)));
188  TCandidate *mumfit = fitter.FittedCand(*(psi2s[j].Daughter(0)->Daughter(1)));
189 
190  TLorentzVector tlvmupf = mupfit->P4();
191  TLorentzVector tlvmumf = mumfit->P4();
192 
193  hjpsim_4cf->Fill((tlvmupf+tlvmumf).M());
194  }
195 
196  // do mass constraint fit
197  for (j=0;j<jpsi.GetLength();++j)
198  {
199  RhoKinFitter mfitter(jpsi[j]);
200  mfitter.AddMassConstraint(3.0965);
201  mfitter.Fit();
202  double chi2_m = mfitter.GlobalChi2();
203  hjpsi_chi2_mf->Fill(chi2_m);
204 
205  if (chi2_m<2) hjpsim_mcf->Fill(jpsi[j].M());
206  }
207 
208  // *** do MC truth match for PID type
209  SelectPdgCode(mctrk, muplus);
210  SelectPdgCode(mctrk, muminus);
211  SelectPdgCode(mctrk, piplus);
212  SelectPdgCode(mctrk, piminus);
213 
214  // *** all combinatorics again with true PID
215  jpsi.Combine(muplus, muminus);
216  for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j].M() );
217  //jpsi.Select(jpsiMassSel);
218 
219  psi2s.Combine(jpsi, piplus, piminus);
220  for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j].M() );
221 
222  // *** and again with PidAlgoEmcBayes and loose selection
223  theAnalysis->FillList(muplus, "MuonLoosePlus","PidAlgoMvd;PidAlgoStt");
224  theAnalysis->FillList(muminus, "MuonLooseMinus","PidAlgoMvd;PidAlgoStt");
225  theAnalysis->FillList(piplus, "PionLoosePlus","PidAlgoMvd;PidAlgoStt");
226  theAnalysis->FillList(piminus, "PionLooseMinus","PidAlgoMvd;PidAlgoStt");
227 
228  jpsi.Combine(muplus, muminus);
229  for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j].M() );
230  //jpsi.Select(jpsiMassSel);
231 
232  psi2s.Combine(jpsi, piplus, piminus);
233  for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j].M() );
234 
235  }
236 
237  // *** write out all the histos
238  out->cd();
239 
240  hjpsim_nopid->Write();
241  hpsim_nopid->Write();
242  hjpsim_lpid->Write();
243  hpsim_lpid->Write();
244  hjpsim_tpid->Write();
245  hpsim_tpid->Write();
246  hjpsim_ftm->Write();
247  hpsim_ftm->Write();
248 
249  hpsim_diff->Write();
250  hjpsim_diff->Write();
251 
252  hjpsim_vf->Write();
253  hjpsim_4cf->Write();
254  hjpsim_mcf->Write();
255 
256  hjpsi_chi2_vf->Write();
257  hpsi_chi2_4c->Write();
258  hjpsi_chi2_mf->Write();
259 
260  hvpos->Write();
261 
262  out->Save();
263 
264 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
i
Int_t i
Definition: run_full.C:25
inRecoFile
TString inRecoFile
Definition: macro/examples/pid/run_pid_stt.C:15
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
SelectPdgCode
int SelectPdgCode(TCandList &mct, TCandList &l)
Definition: anatut_psi2s.C:9
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
muplus
CandList muplus
Definition: softtrigger_kin5.C:144
inSimFile
TString inSimFile
Definition: run_reco_llbar_evtgen.C:8
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
out
TFile * out
Definition: reco_muo.C:20
piminus
CandList piminus
Definition: softtrigger_kin5.C:144
PndAnalysis::GetEvent
Int_t GetEvent(Int_t n=-1)
Definition: PndAnalysis.cxx:257
muminus
CandList muminus
Definition: softtrigger_kin5.C:144
int SelectPdgCode ( TCandList &  mct,
TCandList &  l 
)

Definition at line 9 of file anatut_psi2s.C.

Referenced by anatut_psi2s().

10 {
11  int removed = 0;
12  int pdgcode=0;
13 
14  PndMcTruthMatch mcm;
15 
16  if (l.GetLength()>0) pdgcode = l[0].PdgCode();
17  for (int ii=l.GetLength()-1;ii>=0;--ii)
18  {
19  if (!mcm.MctMatch(l[ii],mct))
20  {
21  l.Remove(l[ii]);
22  removed++;
23  }
24  }
25  return removed;
26 }
mct
PndMCTrack * mct
Definition: hist-t7.C:147
Generated on Wed Apr 3 2019 10:03:10 for FairRoot/PandaRoot by   doxygen 1.8.5