FairRoot/PandaRoot
ana_multikalman.C
Go to the documentation of this file.
1 class RhoCandList;
2 class RhoCandidate;
3 class RhoTuple;
4 class PndRhoTupleQA;
5 class PndAnaPidSelector;
6 class PndAnaPidCombiner;
7 class PndAnalysis;
8 
9 // *** routine to only keep PID matched candidates in list
10 int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
11 {
12  int removed = 0;
13 
14  for (int ii=l.GetLength()-1;ii>=0;--ii)
15  {
16  if ( !(ana->McTruthMatch(l[ii])) )
17  {
18  l.Remove(l[ii]);
19  removed++;
20  }
21  }
22 
23  return removed;
24 }
25 
26 void printCand(RhoCandidate *c)
27 {
28  TLorentzVector lv=c->P4();
29 
30  cout <<c->PdgCode()<<" ("<<lv.X()<<"/"<<lv.Y()<<"/"<<lv.Z()<<"/"<<lv.E()<<")"<<endl;
31 }
32 
33 void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
34 {
35  int n_smc = 0;
36  int n_strk = 0;
37  int n_both = 0;
38  double d = 0.00001;
39 
40  for (int i=0;i<l.GetLength()-1;++i)
41  {
42  for (int j=i+1;j<l.GetLength();++j)
43  {
44  TLorentzVector dl = l[i]->P4() - l[j]->P4();
45 
46  bool chkmc = (l[i]->GetMcTruth()==l[j]->GetMcTruth());
47  bool chktrk = (fabs(dl.X())<d) && (fabs(dl.Y())<d) && (fabs(dl.Z())<d) && (fabs(dl.E())<d);
48  if (chkmc) n_smc++;
49  if (chktrk) n_strk++;
50  if (chktrk && chkmc) n_both++;
51  }
52  }
53  n1 = n_strk;
54  n2 = n_smc;
55  n3 = n_both;
56 }
57 
58 int ana_multikalman(int nevts=0)
59 {
60  //-----User Settings:------------------------------------------------------
61  TString parAsciiFile = "all.par";
62 
63  TString prefix = "jpsi";
64  TString input = "psi2s_Jpsi2pi_Jpsi_mumu.dec";
65  TString output = "ana";
66  TString friend1 = "reco";
67  TString friend2 = "pid";
68 
69 // TString prefix = "evtcomplete";
70 // TString input = "psi2s_Jpsi2pi_Jpsi_mumu.dec";
71 // TString output = "ana";
72 // TString friend1 = "reco_single";
73 // TString friend2 = "pid_single";
74 
75  TString friend3 = "";
76  TString friend4 = "";
77 
78  // // *** some variables
79  int i=0,j=0, k=0, l=0;
80  gStyle->SetOptFit(1011);
81 
82  // ----- Initial Settings --------------------------------------------
83 // PndMasterRunAna *fRun= new PndMasterRunAna();
84 // // fRun->SetInput(input);
85 // fRun->SetInput(friend2);
86 // fRun->SetOutput(output);
91 // fRun->SetParamAsciiFile(parAsciiFile);
92 // fRun->Setup(prefix);
93 //
94 //
95 // fRun->Init();
96 //
97 
98 
99 
100 
101 
102 
103 
104  // *** the output file for FairRunAna
105  TString OutFile="out_dummy.root";
106 
107  // *** the files coming from the simulation
108  TString inPidFile = prefix+"_pid.root"; // this file contains the PndPidCandidates and McTruth
109  TString inParFile = prefix+"_par.root";
110 
111  // *** PID table with selection thresholds; can be modified by the user
112  TString pidParFile = TString(gSystem->Getenv("VMCWORKDIR"))+"/macro/params/all.par";
113 
114  // *** initialization
115  FairLogger::GetLogger()->SetLogToFile(kFALSE);
116  FairRunAna* fRun = new FairRunAna();
117  FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
118  fRun->SetSource(new FairFileSource(inPidFile));
119  fRun->SetUseFairLinks(kTRUE);
120 
121  // *** setup parameter database
122  FairParRootFileIo* parIO = new FairParRootFileIo();
123  parIO->open(inParFile);
124  FairParAsciiFileIo* parIOPid = new FairParAsciiFileIo();
125  parIOPid->open(pidParFile.Data(),"in");
126 
127  rtdb->setFirstInput(parIO);
128  rtdb->setSecondInput(parIOPid);
129  rtdb->setOutput(parIO);
130 
131  fRun->SetOutputFile(OutFile);
132  fRun->Init();
133 
134 
135 
136 
137 
138 
139 
140 
141 
142 
143  // *** create rhotuples
144  RhoTuple * fEPlus = new RhoTuple("EPlus","EPlus Info");
145  RhoTuple * fEMinus = new RhoTuple("EMinus","EMinus Info");
146  RhoTuple * fMuPlus = new RhoTuple("MuPlus","MuPlus Info");
147  RhoTuple * fMuMinus = new RhoTuple("MuMinus","MuMinus Info");
148  RhoTuple * fPiPlus = new RhoTuple("PiPlus","PiPlus Info");
149  RhoTuple * fPiMinus = new RhoTuple("PiMinus","PiMinus Info");
150  RhoTuple * fKPlus = new RhoTuple("KPlus","KPlus Info");
151  RhoTuple * fKMinus = new RhoTuple("KMinus","KMinus Info");
152  RhoTuple * fPPlus = new RhoTuple("PPlus","PPlus Info");
153  RhoTuple * fPMinus = new RhoTuple("PMinus","PMinus Info");
154 
155  RhoTuple * fMcTruth = new RhoTuple("McTruth","McTruth Info");
156 
157  // *** create an output file for all histograms
158  TFile *out = TFile::Open(prefix+"_output_ana.root","RECREATE");
159 
160  //
161  // Now the analysis stuff comes...
162  //
163 
164 
165  // *** the data reader object
166  //Set true for multikalman, false (or not set at all) for standard single kalman
167  PndAnalysis* theAnalysis = new PndAnalysis();
168  if (nevts==0) nevts= theAnalysis->GetEntries();
169 
170  // *** RhoCandLists for the analysis
171  RhoCandList muplus, muminus, piplus, piminus,eplus,eminus,kplus,kminus,pplus,pminus, mctruth;
172 
173  // *** Mass selector for the jpsi cands
174  double m0_jpsi = TDatabasePDG::Instance()->GetParticle("J/psi")->Mass(); // Get nominal PDG mass of the J/psi
175  RhoMassParticleSelector *jpsiMassSel=new RhoMassParticleSelector("jpsi",m0_jpsi,1.0);
176 
177  // *** the lorentz vector of the initial psi(2S)
178  TLorentzVector ini(0, 0, 6.231552, 7.240065);
179 
180  // ***
181  // the event loop
182  // ***
183 
184  int cntdbltrk=0, cntdblmc=0, cntdblboth=0, cnttrk=0, cnt_dbl_jpsi=0, cnt_dbl_psip=0;
185 
186  while (theAnalysis->GetEvent() && i++<nevts)
187  {
188  if ((i%100)==0) cout<<"evt " << i << endl;
189  PndRhoTupleQA qa(theAnalysis, ini.Z());
190 
191  theAnalysis->FillList(mctruth, "McTruth");
192 
193  theAnalysis->FillList(eplus, "PionAllPlus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
194  theAnalysis->FillList(eminus, "PionAllMinus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
195  theAnalysis->FillList(muplus, "MuonAllPlus","PidAlgoMdtHardCuts");
196  theAnalysis->FillList(muminus, "MuonAllMinus","PidAlgoMdtHardCuts");
197  theAnalysis->FillList(piplus, "PionAllPlus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
198  theAnalysis->FillList(piminus, "PionAllMinus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
199  theAnalysis->FillList(kplus, "PionAllPlus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
200  theAnalysis->FillList(kminus, "PionAllMinus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
201  theAnalysis->FillList(pplus, "PionAllPlus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
202  theAnalysis->FillList(pminus, "PionAllMinus","PidAlgoStt;PidAlgoMvd;PidAlgoEmcBayes;PidAlgoDrc;PidAlgoDisc");
203 
204  for (j=0;j<mctruth.GetLength();++j)
205  {
206  qa.qaCand("",mctruth[j],fMcTruth);
207  fMcTruth->DumpData();
208  }
209 
210  // *** momentum and theta histograms
211  for (j=0;j<eplus.GetLength();++j)
212  {
213  fEPlus->Column("mct",theAnalysis->McTruthMatch(eplus[j]));
214  fEPlus->Column("recocandfitstatus",eplus[j]->GetRecoCandidate()->GetFitStatus());
215  qa.qaCand("",eplus[j],fEPlus);
216  qa.qaMc("",eplus[j],fEPlus,false);
217  qa.qaMcDiff("",eplus[j],fEPlus,false);
218  qa.qaPid("",eplus[j],fEPlus);
219  qa.qaTrk("",eplus[j],fEPlus);
220  fEPlus->DumpData();
221  }
222  for (j=0;j<eminus.GetLength();++j)
223  {
224  fEMinus->Column("mct",theAnalysis->McTruthMatch(eminus[j]));
225  fEMinus->Column("recocandfitstatus",eminus[j]->GetRecoCandidate()->GetFitStatus());
226  qa.qaCand("",eminus[j],fEMinus);
227  qa.qaMc("",eminus[j],fEMinus,false);
228  qa.qaMcDiff("",eminus[j],fEMinus,false);
229  qa.qaPid("",eminus[j],fEMinus);
230  qa.qaTrk("",eminus[j],fEMinus);
231  fEMinus->DumpData();
232  }
233  for (j=0;j<muplus.GetLength();++j)
234  {
235  fMuPlus->Column("mct",theAnalysis->McTruthMatch(muplus[j]));
236  fMuPlus->Column("recocandfitstatus",muplus[j]->GetRecoCandidate()->GetFitStatus());
237  qa.qaCand("",muplus[j],fMuPlus);
238  qa.qaMc("",muplus[j],fMuPlus,false);
239  qa.qaMcDiff("",muplus[j],fMuPlus,false);
240  qa.qaPid("",muplus[j],fMuPlus);
241  qa.qaTrk("",muplus[j],fMuPlus);
242  fMuPlus->DumpData();
243  }
244  for (j=0;j<muminus.GetLength();++j)
245  {
246  fMuMinus->Column("mct",theAnalysis->McTruthMatch(muminus[j]));
247  fMuMinus->Column("recocandfitstatus",muminus[j]->GetRecoCandidate()->GetFitStatus());
248  qa.qaCand("",muminus[j],fMuMinus);
249  qa.qaMc("",muminus[j],fMuMinus,false);
250  qa.qaMcDiff("",muminus[j],fMuMinus,false);
251  qa.qaPid("",muminus[j],fMuMinus);
252  qa.qaTrk("",muminus[j],fMuMinus);
253  fMuMinus->DumpData();
254  }
255  for (j=0;j<piplus.GetLength();++j)
256  {
257  fPiPlus->Column("mct",theAnalysis->McTruthMatch(piplus[j]));
258  fPiPlus->Column("recocandfitstatus",piplus[j]->GetRecoCandidate()->GetFitStatus());
259  qa.qaCand("",piplus[j],fPiPlus);
260  qa.qaMc("",piplus[j],fPiPlus,false);
261  qa.qaMcDiff("",piplus[j],fPiPlus,false);
262  qa.qaPid("",piplus[j],fPiPlus);
263  qa.qaTrk("",piplus[j],fPiPlus);
264  fPiPlus->DumpData();
265  }
266  for (j=0;j<piminus.GetLength();++j)
267  {
268  fPiMinus->Column("mct",theAnalysis->McTruthMatch(piminus[j]));
269  fPiMinus->Column("recocandfitstatus",piminus[j]->GetRecoCandidate()->GetFitStatus());
270  qa.qaCand("",piminus[j],fPiMinus);
271  qa.qaMc("",piminus[j],fPiMinus,false);
272  qa.qaMcDiff("",piminus[j],fPiMinus,false);
273  qa.qaPid("",piminus[j],fPiMinus);
274  qa.qaTrk("",piminus[j],fPiMinus);
275  fPiMinus->DumpData();
276  }
277  for (j=0;j<kplus.GetLength();++j)
278  {
279  fKPlus->Column("mct",theAnalysis->McTruthMatch(kplus[j]));
280  fKPlus->Column("recocandfitstatus",kplus[j]->GetRecoCandidate()->GetFitStatus());
281  qa.qaCand("",kplus[j],fKPlus);
282  qa.qaMc("",kplus[j],fKPlus,false);
283  qa.qaMcDiff("",kplus[j],fKPlus,false);
284  qa.qaPid("",kplus[j],fKPlus);
285  qa.qaTrk("",kplus[j],fKPlus);
286  fKPlus->DumpData();
287  }
288  for (j=0;j<kminus.GetLength();++j)
289  {
290  fKMinus->Column("mct",theAnalysis->McTruthMatch(kminus[j]));
291  fKMinus->Column("recocandfitstatus",kminus[j]->GetRecoCandidate()->GetFitStatus());
292  qa.qaCand("",kminus[j],fKMinus);
293  qa.qaMc("",kminus[j],fKMinus,false);
294  qa.qaMcDiff("",kminus[j],fKMinus,false);
295  qa.qaPid("",kminus[j],fKMinus);
296  qa.qaTrk("",kminus[j],fKMinus);
297  fKMinus->DumpData();
298  }
299  for (j=0;j<pplus.GetLength();++j)
300  {
301  fPPlus->Column("mct",theAnalysis->McTruthMatch(pplus[j]));
302  fPPlus->Column("recocandfitstatus",pplus[j]->GetRecoCandidate()->GetFitStatus());
303  qa.qaCand("",pplus[j],fPPlus);
304  qa.qaMc("",pplus[j],fPPlus,false);
305  qa.qaMcDiff("",pplus[j],fPPlus,false);
306  qa.qaPid("",pplus[j],fPPlus);
307  qa.qaTrk("",pplus[j],fPPlus);
308  fPPlus->DumpData();
309  }
310  for (j=0;j<pminus.GetLength();++j)
311  {
312  fPMinus->Column("mct",theAnalysis->McTruthMatch(pminus[j]));
313  fPMinus->Column("recocandfitstatus",pminus[j]->GetRecoCandidate()->GetFitStatus());
314  qa.qaCand("",pminus[j],fPMinus);
315  qa.qaMc("",pminus[j],fPMinus,false);
316  qa.qaMcDiff("",pminus[j],fPMinus,false);
317  qa.qaPid("",pminus[j],fPMinus);
318  qa.qaTrk("",pminus[j],fPMinus);
319  fPMinus->DumpData();
320  }
321 
322  }
323 
324  // *** write out all the histos
325  out->cd();
326 
327  fEPlus->GetInternalTree()->Write();
328  fEMinus->GetInternalTree()->Write();
329  fMuPlus->GetInternalTree()->Write();
330  fMuMinus->GetInternalTree()->Write();
331  fPiPlus->GetInternalTree()->Write();
332  fPiMinus->GetInternalTree()->Write();
333  fKPlus->GetInternalTree()->Write();
334  fKMinus->GetInternalTree()->Write();
335  fPPlus->GetInternalTree()->Write();
336  fPMinus->GetInternalTree()->Write();
337  fMcTruth->GetInternalTree()->Write();
338 
339  out->Save();
340 
341  if (gROOT->GetVersionInt() >= 60602) {
342  gGeoManager->GetListOfVolumes()->Delete();
343  gGeoManager->GetListOfShapes()->Delete();
344  delete gGeoManager;
345  }
346  return 0;
347 
348 }
PndAnalysis::GetEntries
Int_t GetEntries()
Definition: PndAnalysis.cxx:471
pminus
CandList pminus
Definition: softtrigger_kin5.C:144
printCand
void printCand(RhoCandidate *c)
Definition: ana_multikalman.C:26
ana_multikalman
int ana_multikalman(int nevts=0)
Definition: ana_multikalman.C:58
d
TObjArray * d
Definition: Pnd_Hc_etaee7G.C:78
i
Int_t i
Definition: run_full.C:25
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
kminus
CandList kminus
Definition: softtrigger_kin5.C:144
piplus
CandList piplus
Definition: softtrigger_kin5.C:144
SelectTruePid
int SelectTruePid(PndAnalysis *ana, RhoCandList &l)
Definition: ana_multikalman.C:10
eplus
CandList eplus
Definition: softtrigger_kin5.C:144
gGeoManager
TGeoManager * gGeoManager
Definition: createRootGeoFile.C:97
output
FairParRootFileIo * output
Definition: sim_emc_apd.C:120
fRun
FairRunAna * fRun
Definition: hit_dirc.C:58
muplus
CandList muplus
Definition: softtrigger_kin5.C:144
countDoubles
void countDoubles(RhoCandList &l, int &n1, int &n2, int &n3)
Definition: ana_multikalman.C:33
PndAnalysis::PndAnalysis
PndAnalysis(TString tname1="", TString tname2="", TString algnamec="PidAlgoIdealCharged", TString algnamen="PidAlgoIdealNeutral")
Definition: PndAnalysis.cxx:48
rtdb
FairRuntimeDb * rtdb
Definition: hit_dirc.C:66
pplus
CandList pplus
Definition: softtrigger_kin5.C:144
RhoCandidate::P4
TLorentzVector P4() const
Definition: RhoCandidate.h:195
RhoTuple::Column
void Column(const char *label, Bool_t value, Bool_t defval=0, const char *block=0)
Definition: RhoTuple.cxx:56
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
out
TFile * out
Definition: reco_muo.C:20
RhoCandList::Remove
Int_t Remove(RhoCandidate *)
Definition: RhoCandList.cxx:131
c
GeV c
Definition: ReadMesDecCorrRes_PiT.C:28
RhoTuple::DumpData
void DumpData()
Definition: RhoTuple.cxx:391
piminus
CandList piminus
Definition: softtrigger_kin5.C:144
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
muminus
CandList muminus
Definition: softtrigger_kin5.C:144
kplus
CandList kplus
Definition: softtrigger_kin5.C:144
eminus
CandList eminus
Definition: softtrigger_kin5.C:144
Generated on Sun Aug 18 2019 03:01:05 for FairRoot/PandaRoot by   doxygen 1.8.8