FairRoot/PandaRoot
Functions
thailand2017/solution/tut_ana_fast.C File Reference
#include "auxtut.C"

Go to the source code of this file.

Functions

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

Function Documentation

int SelectTruePid ( PndAnalysis *  ana,
RhoCandList &  l 
)

Definition at line 18 of file thailand2017/solution/tut_ana_fast.C.

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

19 {
20  int removed = 0;
21 
22  for (int ii=l.GetLength()-1;ii>=0;--ii)
23  {
24  if ( !(ana->McTruthMatch(l[ii])) )
25  {
26  l.Remove(l[ii]);
27  removed++;
28  }
29  }
30 
31  return removed;
32 }
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_fast ( int  nevts = 0,
TString  prefix = "signal" 
)

Definition at line 35 of file thailand2017/solution/tut_ana_fast.C.

References RhoKinFitter::Add4MomConstraint(), RhoKinFitter::AddMassConstraint(), RhoCandList::Combine(), PndAnalysis::FillList(), RhoKinFitter::Fit(), RhoFitterBase::Fit(), RhoCalculationTools::ForceConstantBz(), fRun, RhoFitterBase::GetChi2(), PndAnalysis::GetEntries(), PndAnalysis::GetEvent(), RhoCandList::GetLength(), RhoFitterBase::GetProb(), hvpos, i, RhoCandidate::M(), PndAnalysis::McTruthMatch(), muminus, muplus, out, piminus, piplus, plotmyhistos(), PndAnalysis::PndAnalysis(), RhoCandidate::Pos(), RhoCandList::Select(), SelectTruePid(), RhoCandList::SetType(), TString, and writemyhistos().

36 {
37  // *** some variables
38  int i=0,j=0, k=0, l=0;
39  gStyle->SetOptFit(1011);
40 
41  // *** the output file for FairRunAna
42  TString OutFile="out_dummy.root";
43 
44  // *** the files coming from the simulation
45  TString inPidFile = prefix+"_fast.root"; // this file contains the PndPidCandidates and McTruth
46 
47  // *** PID table with selection thresholds; can be modified by the user
48  TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";
49 
50  // *** initialization
51  FairRunAna* fRun = new FairRunAna();
52  fRun->SetSource(new FairFileSource(inPidFile));
53  fRun->SetOutputFile(OutFile);
54 
55  // *** take constant field; needed for PocaVtx
56  RhoCalculationTools::ForceConstantBz(20.0);
57 
58  fRun->Init();
59 
60  // *** create an output file for all histograms
61  TFile *out = TFile::Open(prefix+"_ana_fast.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  TString pidalgo = "PidChargedProbability";
121 
122  // ***
123  // the event loop
124  // ***
125  while (theAnalysis->GetEvent() && i++<nevts)
126  {
127  if ((i%100)==0) cout<<"evt " << i << endl;
128 
129  // *** Select with no PID info ('All'); type and mass are set
130  theAnalysis->FillList(muplus, "MuonAllPlus" , pidalgo );
131  theAnalysis->FillList(muminus, "MuonAllMinus" , pidalgo );
132  theAnalysis->FillList(piplus, "PionAllPlus" , pidalgo );
133  theAnalysis->FillList(piminus, "PionAllMinus" , pidalgo );
134 
135  // *** combinatorics for J/psi -> mu+ mu-
136  jpsi.Combine(muplus, muminus);
137 
138 
139  // ***
140  // *** do the TRUTH MATCH for jpsi
141  // ***
142  jpsi.SetType(443);
143 
144  for (j=0;j<jpsi.GetLength();++j)
145  {
146  hjpsim_all->Fill( jpsi[j]->M() );
147 
148  if (theAnalysis->McTruthMatch(jpsi[j]))
149  {
150  hjpsim_ftm->Fill( jpsi[j]->M() );
151  hjpsim_diff->Fill( jpsi[j]->GetMcTruth()->M() - jpsi[j]->M() );
152  }
153  else
154  hjpsim_nm->Fill( jpsi[j]->M() );
155  }
156 
157  // ***
158  // *** do VERTEX FIT (J/psi)
159  // ***
160  for (j=0;j<jpsi.GetLength();++j)
161  {
162  RhoKinVtxFitter vtxfitter(jpsi[j]); // instantiate a vertex fitter
163  vtxfitter.Fit();
164 
165  double chi2_vtx = vtxfitter.GetChi2(); // access chi2 of fit
166  double prob_vtx = vtxfitter.GetProb(); // access probability of fit
167  hjpsi_chi2_vf->Fill(chi2_vtx);
168  hjpsi_prob_vf->Fill(prob_vtx);
169 
170  if ( prob_vtx > 0.01 ) // when good enough, fill some histos
171  {
172  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
173  TVector3 jVtx=jfit->Pos(); // and the decay vertex position
174 
175  hjpsim_vf->Fill(jfit->M());
176  hvpos->Fill(jVtx.X(),jVtx.Y());
177  }
178  }
179 
180  // *** some rough mass selection
181  jpsi.Select(jpsiMassSel);
182 
183  // *** combinatorics for psi(2S) -> J/psi pi+ pi-
184  psi2s.Combine(jpsi, piplus, piminus);
185 
186 
187  // ***
188  // *** do the TRUTH MATCH for psi(2S)
189  // ***
190  psi2s.SetType(88888);
191 
192  for (j=0;j<psi2s.GetLength();++j)
193  {
194  hpsim_all->Fill( psi2s[j]->M() );
195 
196  if (theAnalysis->McTruthMatch(psi2s[j]))
197  {
198  hpsim_ftm->Fill( psi2s[j]->M() );
199  hpsim_diff->Fill( psi2s[j]->GetMcTruth()->M() - psi2s[j]->M() );
200  }
201  else
202  hpsim_nm->Fill( psi2s[j]->M() );
203  }
204 
205 
206  // ***
207  // *** do 4C FIT (initial psi(2S) system)
208  // ***
209  for (j=0;j<psi2s.GetLength();++j)
210  {
211  RhoKinFitter fitter(psi2s[j]); // instantiate the kin fitter in psi(2S)
212  fitter.Add4MomConstraint(ini); // set 4 constraint
213  fitter.Fit(); // do fit
214 
215  double chi2_4c = fitter.GetChi2(); // get chi2 of fit
216  double prob_4c = fitter.GetProb(); // access probability of fit
217  hpsi_chi2_4c->Fill(chi2_4c);
218  hpsi_prob_4c->Fill(prob_4c);
219 
220  if ( prob_4c > 0.01 ) // when good enough, fill some histo
221  {
222  RhoCandidate *jfit = psi2s[j]->Daughter(0)->GetFit(); // get fitted J/psi
223 
224  hjpsim_4cf->Fill(jfit->M());
225  }
226  }
227 
228 
229  // ***
230  // *** do MASS CONSTRAINT FIT (J/psi)
231  // ***
232  for (j=0;j<jpsi.GetLength();++j)
233  {
234  RhoKinFitter mfitter(jpsi[j]); // instantiate the RhoKinFitter in psi(2S)
235  mfitter.AddMassConstraint(m0_jpsi); // add the mass constraint
236  mfitter.Fit(); // do fit
237 
238  double chi2_m = mfitter.GetChi2(); // get chi2 of fit
239  double prob_m = mfitter.GetProb(); // access probability of fit
240  hjpsi_chi2_mf->Fill(chi2_m);
241  hjpsi_prob_mf->Fill(prob_m);
242 
243  if ( prob_m > 0.01 ) // when good enough, fill some histo
244  {
245  RhoCandidate *jfit = jpsi[j]->GetFit(); // access the fitted cand
246  hjpsim_mcf->Fill(jfit->M());
247  }
248  }
249 
250 
251  // ***
252  // *** TRUE PID combinatorics
253  // ***
254 
255  // *** do MC truth match for PID type
256  SelectTruePid(theAnalysis, muplus);
257  SelectTruePid(theAnalysis, muminus);
258  SelectTruePid(theAnalysis, piplus);
259  SelectTruePid(theAnalysis, piminus);
260 
261  // *** all combinatorics again with true PID
262  jpsi.Combine(muplus, muminus);
263  for (j=0;j<jpsi.GetLength();++j) hjpsim_trpid->Fill( jpsi[j]->M() );
264  jpsi.Select(jpsiMassSel);
265 
266  psi2s.Combine(jpsi, piplus, piminus);
267  for (j=0;j<psi2s.GetLength();++j) hpsim_trpid->Fill( psi2s[j]->M() );
268 
269 
270  // ***
271  // *** LOOSE PID combinatorics
272  // ***
273  pidalgo = "MvdPidProbability;SttPidProbability;DrcBarrelProbability";
274  // *** and again with Mvd, Stt, Drc and loose selection
275  theAnalysis->FillList(muplus, "MuonLoosePlus", pidalgo );
276  theAnalysis->FillList(muminus, "MuonLooseMinus", pidalgo );
277  theAnalysis->FillList(piplus, "PionLoosePlus", pidalgo );
278  theAnalysis->FillList(piminus, "PionLooseMinus", pidalgo);
279 
280  jpsi.Combine(muplus, muminus);
281  for (j=0;j<jpsi.GetLength();++j) hjpsim_lpid->Fill( jpsi[j]->M() );
282  jpsi.Select(jpsiMassSel);
283 
284  psi2s.Combine(jpsi, piplus, piminus);
285  for (j=0;j<psi2s.GetLength();++j) hpsim_lpid->Fill( psi2s[j]->M() );
286 
287 
288  // ***
289  // *** TIGHT PID combinatorics
290  // ***
291  TString pidalgomdt = "ScMdtPidBarrelProbability;ScMdtPidForwardProbability";
292  // *** and again with Mvd, Stt and tight selection
293  theAnalysis->FillList(muplus, "MuonTightPlus", pidalgomdt);
294  theAnalysis->FillList(muminus, "MuonTightMinus", pidalgomdt);
295  theAnalysis->FillList(piplus, "PionLoosePlus", pidalgo);
296  theAnalysis->FillList(piminus, "PionLooseMinus", pidalgo);
297 
298  jpsi.Combine(muplus, muminus);
299  for (j=0;j<jpsi.GetLength();++j) hjpsim_tpid->Fill( jpsi[j]->M() );
300  jpsi.Select(jpsiMassSel);
301 
302  psi2s.Combine(jpsi, piplus, piminus);
303  for (j=0;j<psi2s.GetLength();++j) hpsim_tpid->Fill( psi2s[j]->M() );
304 
305  }
306 
307  // *** change to directory where histograms are created
308  out->cd();
309 
310  // *** plot all histos
311  plotmyhistos();
312 
313  // *** write out all the histos to file
314  int nhist = writemyhistos();
315  cout<<"Writing "<<nhist<<" histograms to file"<<endl;
316  out->Save();
317 
318  // *** temporaty fix to avoid error on macro exit
319  //RemoveGeoManager();
320 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
SelectTruePid
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
Definition: master/ana_complete.C:8
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
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
RhoCalculationTools::ForceConstantBz
static void ForceConstantBz(Double_t bz=0.)
Force a constant B field value for all positions.
Definition: RhoCalculationTools.h:37
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
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
writemyhistos
int writemyhistos(int maxy=800, double asp=1.1)
Definition: QA/auxi.C:121
plotmyhistos
void plotmyhistos(std::vector< TH1 * > h, int maxy=700, double asp=1.1)
Definition: QA/auxi.C:62
Generated on Wed Apr 3 2019 08:44:25 for FairRoot/PandaRoot by   doxygen 1.8.5