FairRoot/PandaRoot
run_ana_eta_c_stt.C
Go to the documentation of this file.
1 class TCandList;
2 class TCandidate;
3 class TFitParams;
4 
5 int run_ana_eta_c_stt(int nevts=0, bool usePID=true)
6 {
7  TString OutFile="etac_histo_stt.root";
8 
9  enum DetectorId {kDCH,kDRC,kDSK,kEMC,kGEM,kLUMI,kMDT,
10  kMVD,kRPC,kSTT,kTPC,kTOF,kFTS,kHYPG,kHYP};
11  gStyle->SetOptFit(1011);
12 
13  TStopwatch timer;
14  timer.Start();
15 
16  gROOT->Macro("$VMCWORKDIR/gconfig/rootlogon.C");
17 
18  TString inPidFile = "evt_pid_stt.root";
19  TString inSimFile = "evt_points_stt.root";
20 
21  TFile *inFile = TFile::Open(inSimFile,"READ");
22  TTree *tree=(TTree *) inFile->Get("pndsim") ;
23  tree->AddFriend("pndsim",inPidFile);
24 
25  TClonesArray* mc_array=new TClonesArray("PndMCTrack");
26  tree->SetBranchAddress("MCTrack",&mc_array);
27 
28  TClonesArray* cand_array=new TClonesArray("PndPidCandidate");
29  tree->SetBranchAddress("PidChargedCand", &cand_array);
30 
31 
32  FairMCEventHeader* evthead;
33  tree->SetBranchAddress("MCEventHeader.", &evthead);
34 
35  TFile *out = TFile::Open(OutFile,"RECREATE");
36 
37  // the PndEventReader takes care about file/event handling
38  PndEventReader evr(inPidFile);
39 
40  TH1F *h_etac_nocut=new TH1F("h_etac_nocut","#eta_{c}m(#phi,#phi), (no cuts);E, GeV",100,2.5,3.5);
41  TH1F *h_mphi_nocuts=new TH1F("h_mphi_nocuts","#phi: m(K+ K-) (no cuts);E, GeV",100,0.95,1.5);
42 
43  TH1F *h_etac_pid=new TH1F("h_etac_pid","#eta_{c}m(#phi,#phi), (MC PID);E, GeV",100,2.5,3.5);
44  TH1F *h_mphi_pid=new TH1F("h_mphi_pid","#phi: m(K+ K-) (MC PID);E, GeV",100,0.95,1.5);
45 
46  TH1F *h_etac_4c=new TH1F("h_etac_4c","#eta_{c}m(#phi,#phi), 4C-fit;E, GeV",100,2.5,3.5);
47  TH1F *h_etac_4c_refit=new TH1F("h_etac_4c_refit","#eta_{c}m(#phi,#phi), 4C-fit;E, GeV",100,2.5,3.5);
48  TH1F *h_mphi_4c=new TH1F("h_mphi_4c","#phi: m(K+ K-) (4C-fit);E, GeV",100,0.95,1.5);
49 
50  TH1F *h_etac_vtx=new TH1F("h_etac_vtx","#eta_{c}m(#phi,#phi), Vertex fit;E, GeV",100,2.5,3.5);
51  TH1F *h_mphi_vtx=new TH1F("h_mphi_vtx","#phi: m(K+ K-) (Vertex fit);E, GeV",100,0.95,1.5);
52 
53  TH1F *h_etac_vtx_2=new TH1F("h_etac_vtx_2","#eta_{c}m(#phi,#phi), Vertex fit;E, GeV",100,2.5,3.5);
54  TH1F *h_mphi_vtx_2=new TH1F("h_mphi_vtx_2","#phi: m(K+ K-) (Vertex fit);E, GeV",100,0.95,1.5);
55 
56  TH1F *h_etac_phimass_4c=new TH1F("h_etac_phimass_4c","#eta_{c}m(#phi,#phi), (cut on #phi mass);E,
57  GeV",100,2.8,3.2);
58  TH1F *h_mphi_final_4c=new TH1F("h_mphi_final_4c","#phi: m(K+ K-);E, GeV",100,0.95,1.1);
59 
60  TH1F *h_etac_phimass_vtx=new TH1F("h_etac_phimass_vtx","#eta_{c}m(#phi,#phi), (cut on #phi mass);E,
61  GeV",100,2.8,3.2);
62  TH1F *h_mphi_final_vtx=new TH1F("h_mphi_final_vtx","#phi: m(K+ K-);E, GeV",100,0.95,1.1);
63 
64  TH1F *h_etac_phimass_vtx_2=new TH1F("h_etac_phimass_vtx_2","#eta_{c}m(#phi,#phi), (cut on #phi mass);E,
65  GeV",100,2.8,3.2);
66  TH1F *h_mphi_final_vtx_2=new TH1F("h_mphi_final_vtx_2","#phi: m(K+ K-);E, GeV",100,0.95,1.1);
67 
68  TH1F *h_etac_phimassfit=new TH1F("h_etac_phimassfit","#eta_{c}m(#phi,#phi);E, GeV",100,2.8,3.2);
69  TH1F *h_mphi_final_massfit=new TH1F("h_mphi_final_massfit","#phi: m(K+ K-);E, GeV",100,0.95,1.1);
70 
71  TH1F *nc=new TH1F("nc","n charged",20,0,20);
72 
73  TH1F *h_chi2_4c=new TH1F("h_chi2_4c","#chi^{2} 4C-fit;#chi^{2}/N_{df}",100,0,100);
74  TH1F *h_chi2b_4c=new TH1F("h_chi2b_4c","#chi^{2} 4C-fit;#chi^{2}/N_{df}",100,0,100);
75  TH1F *h_chi2_vtx=new TH1F("h_chi2_vtx","#chi^{2} vertex;#chi^{2}/N_{df}",100,0,100);
76  TH1F *h_chi2b_vtx=new TH1F("h_chi2b_vtx","#chi^{2} vertex;#chi^{2}/N_{df}",100,0,100);
77 
78  TH1F *h_chi2_prefit=new TH1F("h_chi2_prefit","#chi^{2} prefit;#chi^{2}",100,0,100);
79 
80  TH1F *h_chi2_mass=new TH1F("h_chi2_mass","#chi^{2} Mass constraint fit;#chi^{2}",100,0,100);
81 
82  TH2F *hvpos = new TH2F("hvpos","(x,y) projection of fitted decay
83  vertex",100,-5,5,100,-5,5);
84  TH1F *hvzpos = new TH1F("hvzpos","z position of fitted decay vertex",100,-10,10);
85  TH1F *hvtxresX = new TH1F("hvtxresX","X resolution of fitted decay
86  vertex",100,-0.1,0.1);
87  TH1F *hvtxresY = new TH1F("hvtxresY","Y resolution of fitted decay
88  vertex",100,-0.1,0.1);
89  TH1F *hvtxresZ = new TH1F("hvtxresZ","Z resolution of fitted decay
90  vertex",100,-0.1,0.1);
91 
92  TH2F *h_theta_p = new TH2F("h_theta_p","Theta vs p",100,0,180,100,0.,3);
93  TH2F *h_theta_p_nr = new TH2F("h_theta_p_nr","Theta vs p",100,0,180,100,0.,3); // not reconstructed
94 
95  TH2F *h_dp_p = new TH2F("h_dp_p","Delta p vs p",100,-3.,3,100,0,3);
96  TH2F *h_dp_theta = new TH2F("h_dp_theta","Delta p vs theta",100,-3.,3,100,0,180);
97  TH1F *h_dp=new TH1F("h_dp","Delta p",100,-3.,3.);
98  TH1F *h_dp_low=new TH1F("h_dp_low","Delta p",100,-3.,3);
99  TH1F *h_dp_high=new TH1F("h_dp_high","Delta p",100,-3.,3);
100 
101  TPidMassSelector *phiMassSel=new TPidMassSelector("phi",1.02,0.2);
102 
103  TPidPlusSelector *kplusSel=new TPidPlusSelector("kplus");
104  TPidMinusSelector *kminusSel=new TPidMinusSelector("kminus");
105 
106  // the candidates lists we need
107  TCandList p1, p2, phi1, phi1_pid, phi1_massfit, phi1_massfit_cut, etac, etac_pid, etac_nocut, etac_massfit, etac_massfit_cut;
108  TCandList etac_vtx;
109  int n_reco_4c=0, n_reco_vtx=0, n_reco_vtx_2=0;
110  // Number of events in file and number of reconstructed eta_c to store in root file
111  TH1F *n_events=new TH1F("n_events","total number of events",1,0,1);
112  TH1F *n_etac_4c=new TH1F("n_etac_4c","number of reconstructed eta_c (4C-fit)",1,0,1);
113  TH1F *n_etac_vtx=new TH1F("n_etac_vtx","number of reconstructed eta_c (Vertex fit)",1,0,1);
114  TH1F *n_etac_vtx_2=new TH1F("n_etac_vtx_2","number of reconstructed eta_c (Vertex fit)",1,0,1);
115 
116  TLorentzVector ini(0,0,3.6772,4.7333);
117 
118  if (nevts==0) nevts=evr.GetEntries();
119 
120  n_events->SetBinContent(1,nevts);
121  // cout << "nevts " << nevts << "\n";
122  int i=0,j=0, k=0, l=0;
123 
124 
125  TH1F *h_acc_stt=new TH1F("h_acc_stt","STT acceptance",6,0,6);
126 
127  Double_t mc_mom, mc_theta, rec_mom, rec_theta, delta_p;
128  int nEntries = tree->GetEntriesFast();
129  for (Int_t j=0; j< nEntries; j++)
130  {
131  tree->GetEntry(j);
132  int mccount=0;
133  for (Int_t mc = 0; mc < mc_array->GetEntriesFast(); mc++)
134  {
135  PndMCTrack *mctrack = (PndMCTrack*)mc_array->At(mc);
136  if (mctrack->GetMotherID()!=-1) continue;
137  if (mctrack->GetNPoints(kSTT)) mccount++;
138  mc_mom = mctrack->GetMomentum().Mag();
139  mc_theta = mctrack->GetMomentum().Theta()*TMath::RadToDeg();
140  Int_t cand_mult = 0;
141  for (Int_t pp=0; pp<cand_array->GetEntriesFast(); pp++)
142  {
143  PndPidCandidate *pidCand = (PndPidCandidate*)cand_array->At(pp);
144  if (pidCand->GetMcIndex()!=mc) continue;
145  if ( (cand_mult==0) || ((cand_mult>0) && (fabs(rec_mom-mc_mom)> fabs(pidCand->GetMomentum().Mag()-mc_mom))) )
146  {
147  rec_mom = pidCand->GetMomentum().Mag();
148  rec_theta = pidCand->GetMomentum().Theta();
149  }
150  cand_mult++;
151  }
152 
153  if (cand_mult==0)
154  {
155  h_theta_p_nr->Fill(mc_theta, mc_mom);
156  }
157  delta_p=rec_mom-mc_mom;
158  h_theta_p->Fill(mc_theta, mc_mom);
159  h_dp_p->Fill(delta_p,mc_mom);
160  h_dp_theta->Fill(delta_p,mc_theta);
161  h_dp->Fill(delta_p);
162  if (mc_mom<0.5)
163  h_dp_low->Fill(delta_p);
164  else
165  h_dp_high->Fill(delta_p);
166  }
167  h_acc_stt->Fill(mccount);
168  }
169 
170  // *************
171  // this is the loop through the events ... as simple as this...
172  // ****************
173  while (evr.GetEvent() && i++<nevts)
174  {
175  etac_vtx.Cleanup();
176  phi1_massfit.Cleanup();
177  phi1_massfit_cut.Cleanup();
178  //cout << "evt: " << i << endl;
179  if ((i%100)==0) cout<<"evt " << i << "\n";
180  //if (!((i+1)%100)) cout<<"evt " << i << "\n";
181  evr.FillList(p1,"Charged");
182  evr.FillList(p2,"Charged");
183 
184  p1.Select(kplusSel);
185  p2.Select(kminusSel);
186 
187  int n_removed=0;
188  int ii=0;
189  for (Int_t l=0;l<p1.GetLength();l++){
190  ii=l-n_removed;
191  if((p1[ii].GetMicroCandidate().GetSttHits())==0){
192  p1.Remove(p1[ii]);
193  n_removed++;
194  }
195  }
196 
197  int n_removed=0;
198  int ii=0;
199  for (Int_t l=0;l<p2.GetLength();l++){
200  ii=l-n_removed;
201  if((p2[ii].GetMicroCandidate().GetSttHits())==0){
202  p2.Remove(p2[ii]);
203  n_removed++;
204  }
205  }
206 
207  int nchrg=p1.GetLength()+p2.GetLength();
208  nc->Fill(nchrg);
209 
210  for (j=0;j<p1.GetLength();++j) {
211  p1[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(321)->Mass());
212  }
213  for (j=0;j<p2.GetLength();++j) {
214  p2[j].SetMass(TRho::Instance()->GetPDG()->GetParticle(321)->Mass());
215  }
216 
217  phi1.Combine(p1,p2);
218 
219  for (j=0;j<phi1.GetLength();++j) h_mphi_nocuts->Fill(phi1[j].M());
220 
221  phi1.Select(phiMassSel);
222  etac_nocut.Combine(phi1,phi1);
223 
224  for (l=0;l<etac_nocut.GetLength();++l) {
225  h_etac_nocut->Fill(etac_nocut[l].M());
226  }
227 
228  tree->GetEntry(i-1);
229  TVector3 mcVertex, mcD1Vertex, mcD2vertex;
230  evthead->GetVertex(mcVertex);
231 
232  // MC PID
233  // Leave only kaons in particle lists
234  if (usePID)
235  {
236  int n_removed=0;
237  int ii=0;
238  for (l=0;l<p1.GetLength();++l) {
239  ii=l-n_removed;
240  if (p1[ii].GetMicroCandidate().GetMcIndex()>-1){
241  PndMCTrack *mcTrack =
242  (PndMCTrack*)mc_array->At(p1[ii].GetMicroCandidate().GetMcIndex());
243  if (mcTrack!=0)
244  {
245  if ((mcTrack->GetPdgCode()!=321))
246  {
247  p1.Remove(p1[ii]);
248  n_removed++;
249  }
250  if (mcTrack==0) mcD1Vertex = mcTrack->GetStartVertex();
251  if (mcTrack==2) mcD2Vertex = mcTrack->GetStartVertex();
252  }
253  else
254  {
255  std::cout<<"stt h: " << p1[ii].GetMicroCandidate().GetSttHits() << std::endl;
256  std::cout<<"Kaon list 1, element "<<l<<" has no assosiated mcTRack"<<std::endl;
257  p1.Remove(p1[ii]);
258  n_removed++;
259  }
260  }
261  else {
262  std::cout<<"MC index =-1"<<std::endl;
263  p1.Remove(p1[ii]);
264  n_removed++;
265  }
266  }
267 
268  n_removed=0;
269  ii=0;
270  for (l=0;l<p2.GetLength();++l) {
271  ii=l-n_removed;
272  if (p2[ii].GetMicroCandidate().GetMcIndex()>-1){
273  PndMCTrack *mcTrack =
274  (PndMCTrack*)mc_array->At(p2[ii].GetMicroCandidate().GetMcIndex());
275  if (mcTrack!=0)
276  {
277  if ((mcTrack->GetPdgCode()!=-321))
278  {
279  p2.Remove(p2[ii]);
280  n_removed++;
281  }
282  if (mcTrack==1) mcD1Vertex = mcTrack->GetStartVertex();
283  if (mcTrack==3) mcD2Vertex = mcTrack->GetStartVertex();
284  }
285  else
286  {
287  std::cout<<"stt h: " << p2[ii].GetMicroCandidate().GetSttHits() << std::endl;
288  std::cout<<"Kaon list 2, element "<<l<<" has no assosiated mcTRack"<<std::endl;
289  p2.Remove(p2[ii]);
290  n_removed++;
291  }
292  }
293  else {
294  std::cout<<"MC index =-1"<<std::endl;
295  p2.Remove(p2[ii]);
296  n_removed++;
297  }
298  }
299  }
300 
301  phi1_pid.Combine(p1,p2);
302 
303  for (j=0;j<phi1_pid.GetLength();++j) h_mphi_pid->Fill(phi1_pid[j].M());
304  phi1_pid.Select(phiMassSel);
305  etac.Combine(phi1_pid,phi1_pid);
306 
307  for (l=0;l<etac.GetLength();++l) {
308  h_etac_pid->Fill(etac[l].M());
309  }
310 
312  int best_i=0;
313  double best_chi2=10000;
314  TCandidate *ccfit = new TCandidate();
315  double m_phi1, m_phi2;
316  for (l=0;l<etac.GetLength();++l) {
317  Rho4CFitter fitter(etac[l],ini);
318  fitter.FitConserveMasses();
319  double chi2=fitter.GetChi2();
320  if (chi2<best_chi2)
321  {
322  best_chi2=chi2;
323  best_i = l;
324  ccfit = (TCandidate*)fitter.FittedCand(etac[l]);
325  TCandidate *phi1best = fitter.FittedCand(*(etac[l].Daughter(0)));
326  TCandidate *phi2best = fitter.FittedCand(*(etac[l].Daughter(1)));
327  TCandidate *k1best = fitter.FittedCand(*(etac[l].Daughter(0)->Daughter(0)));
328  TCandidate *k2best = fitter.FittedCand(*(etac[l].Daughter(0)->Daughter(1)));
329  TCandidate *k3best = fitter.FittedCand(*(etac[l].Daughter(1)->Daughter(0)));
330  TCandidate *k4best = fitter.FittedCand(*(etac[l].Daughter(1)->Daughter(1)));
331  TLorentzVector tlvk1 = k1best->P4();
332  TLorentzVector tlvk2 = k2best->P4();
333  TLorentzVector tlvk3 = k3best->P4();
334  TLorentzVector tlvk4 = k4best->P4();
335  m_phi1= (tlvk1+tlvk2).M();
336  m_phi2= (tlvk3+tlvk4).M();
337  }
338  h_chi2_4c->Fill(chi2/9); // Ndf=3N-3=9
339  }
340 
341  if(/*(best_chi2<270)&&*/(etac.GetLength()!=0))
342  {
343  h_chi2b_4c->Fill(best_chi2/9); // Ndf=3N-3=9
344  h_etac_4c->Fill(etac[best_i].M());
345  h_etac_4c_refit->Fill(ccfit->M());
346  h_mphi_4c->Fill(m_phi1);
347  h_mphi_4c->Fill(m_phi2);
348  h_mphi_final_4c->Fill(m_phi1);
349  h_mphi_final_4c->Fill(m_phi2);
350  if (((m_phi1>1.02-0.02)&&(m_phi1<1.02+0.02))
351  && ((m_phi2>1.02-0.02)&&(m_phi2<1.02+0.02)))
352  {
353  h_etac_phimass_4c->Fill(etac[best_i].M());
354  if ((etac[best_i].M()>2.9)&&(etac[best_i].M()<3.06))
355  n_reco_4c++;
356  }
357  }
358 
360  TCandidate *k1, *k2, *k3, *k4, *phi1tmp, *phi2tmp, *etac_tmp;
361 
362  //Combine 4 kaons directly to candidates
363  for (j=0;j<etac.GetLength();++j)
364  {
365  phi1tmp=etac[j].Daughter(0);
366  phi2tmp=etac[j].Daughter(1);
367 
368  k1=phi1tmp->Daughter(0);
369  k2=phi1tmp->Daughter(1);
370  k3=phi2tmp->Daughter(0);
371  k4=phi2tmp->Daughter(1);
372 
373  etac_tmp=k1->Combine(*k2,*k3,*k4);
374  etac_vtx.Add(*etac_tmp);
375  }
376 
377  for (l=0;l<etac_vtx.GetLength();++l) {
378  h_etac_pid->Fill(etac_vtx[l].M());
379  }
380 
381  int best_i=0;
382  double best_chi2=10000;
383  TCandidate *etacfit_best=0;
384  TCandidate *k1fit_best=0, *k2fit_best=0, *k3fit_best=0, *k4fit_best=0;
385  TCandidate *phi1fit_best, *phi2fit_best;
386  TVector3 bestPos;
387  Float_t etacvtx_mass;
388  for (j=0;j<etac_vtx.GetLength();++j)
389  {
390  RhoKinVtxFitter vtxfitter(etac_vtx[j]); // instantiate a vertex fitter
391  vtxfitter.Fit(); // do the vertex fit
392 
393  TCandidate *etacfit=vtxfitter.FittedCand(etac_vtx[j]); // request the fitted EtaC candidate
394  TVector3 etacVtx=etacfit->Pos(); // and the decay vertex position
395  double chi2_vtx=vtxfitter.GlobalChi2();
396  h_chi2_vtx->Fill(chi2_vtx/5); // Number degree of freedom 2N-3=5
397  // plot mass and vtx x,y projection after fit
398  hvpos->Fill(etacVtx.X(),etacVtx.Y());
399  hvzpos->Fill(etacVtx.Z());
400  if(chi2_vtx<best_chi2)
401  {
402  best_chi2=chi2_vtx;
403  best_i=l;
404  etacfit_best=etacfit;
405  k1fit_best=vtxfitter.FittedCand(*(etacfit_best->Daughter(0)));
406  k2fit_best=vtxfitter.FittedCand(*(etacfit_best->Daughter(1)));
407  k3fit_best=vtxfitter.FittedCand(*(etacfit_best->Daughter(2)));
408  k4fit_best=vtxfitter.FittedCand(*(etacfit_best->Daughter(3)));
409  etacvtx_mass = etacfit_best->M();
410  bestPos = etacfit->Pos();
411  }
412  }
413 
414  if(/*(best_chi2<150)&&*/(etac.GetLength()!=0)&&(k1fit_best!=0)&&(k3fit_best!=0))
415  {
416  //h_etac_vtx->Fill(etacfit_best->M());
417  h_chi2b_vtx->Fill(best_chi2/5);
418  h_etac_vtx->Fill(etacvtx_mass);
419  phi1fit_best=k1fit_best->Combine(*k2fit_best);
420  phi2fit_best=k3fit_best->Combine(*k4fit_best);
421  double m_phi1=phi1fit_best->M();
422  double m_phi2=phi2fit_best->M();
423  h_mphi_vtx->Fill(m_phi1);
424  h_mphi_vtx->Fill(m_phi2);
425  h_mphi_final_vtx->Fill(m_phi1);
426  h_mphi_final_vtx->Fill(m_phi2);
427  hvtxresX->Fill(mcVertex.X()-bestPos.X());
428  hvtxresY->Fill(mcVertex.Y()-bestPos.Y());
429  hvtxresZ->Fill(mcVertex.Z()-bestPos.Z());
430 
431  if (((m_phi1>1.02-0.02)&&(m_phi1<1.02+0.02))&&
432  ((m_phi2>1.02-0.02)&&(m_phi2<1.02+0.02)))
433  {
434  h_etac_phimass_vtx->Fill(etacvtx_mass);
435  if ((etacvtx_mass>2.9)&&(etacvtx_mass<3.06))
436  n_reco_vtx++;
437  }
438 
439  }
443  int best_i=0;
444  double best_chi2_prefit=1e6;
445  TCandidate *etacprefit_best=0;
446  TCandidate *k1prefit_best=0, *k2prefit_best=0, *k3prefit_best=0, *k4prefit_best=0;
447  TCandidate *phi1prefit_best, *phi2prefit_best;
448  Double_t mphi1, mphi2, metac;
449  Double_t sigma_etac=0.0331, sigma_phi=0.00392;
450  Double_t m_etac_pdg=2.98, m_phi_pdg=1.02;
451 
452  for (j=0;j<etac_vtx.GetLength();++j)
453  {
454  metac=etac_vtx[j].M();
455  k1=etac_vtx[j].Daughter(0);
456  k2=etac_vtx[j].Daughter(1);
457  k3=etac_vtx[j].Daughter(2);
458  k4=etac_vtx[j].Daughter(3);
459  phi1prefit_best=k1->Combine(*k2);
460  phi2prefit_best=k3->Combine(*k4);
461  m_phi1=phi1prefit_best->M();
462  m_phi2=phi2prefit_best->M();
463 
464  double chi2_prefit=pow(metac-m_etac_pdg,2)/pow(sigma_etac,2)+
465  pow(m_phi1-m_phi_pdg,2)/pow(sigma_phi,2)+pow(m_phi1-m_phi_pdg,2)/pow(sigma_phi,2);
466  h_chi2_prefit->Fill(chi2_prefit);
467 
468  if(chi2_prefit<best_chi2_prefit)
469  {
470  best_chi2_prefit=chi2_prefit;
471  best_i=j;
472  etacprefit_best=etac_vtx[j];
473  etacvtx_mass=etacprefit_best->M();
474  }
475  }
476 
477  if (etacprefit_best!=0)
478  {
479  RhoKinVtxFitter vtxfitter(*etacprefit_best); // instantiate a vertex fitter
480  vtxfitter.Fit(); // do the vertex fit
481  TCandidate *etacfit=vtxfitter.FittedCand(*etacprefit_best);
482  k1prefit_best=vtxfitter.FittedCand(*(etacfit->Daughter(0)));
483  k2prefit_best=vtxfitter.FittedCand(*(etacfit->Daughter(1)));
484  k3prefit_best=vtxfitter.FittedCand(*(etacfit->Daughter(2)));
485  k4prefit_best=vtxfitter.FittedCand(*(etacfit->Daughter(3)));
486  phi1prefit_best=k1prefit_best->Combine(*k2prefit_best);
487  phi2prefit_best=k3prefit_best->Combine(*k4prefit_best);
488  double m_phi1=phi1prefit_best->M();
489  double m_phi2=phi2prefit_best->M();
490  h_mphi_vtx_2->Fill(m_phi1);
491  h_mphi_vtx_2->Fill(m_phi2);
492  h_mphi_final_vtx_2->Fill(m_phi1);
493  h_mphi_final_vtx_2->Fill(m_phi2);
494  if (((m_phi1>1.02-0.02)&&(m_phi1<1.02+0.02))&&
495  ((m_phi2>1.02-0.02)&&(m_phi2<1.02+0.02)))
496  {
497  h_etac_phimass_vtx_2->Fill(etacvtx_mass);
498  if ((etacvtx_mass>2.9)&&(etacvtx_mass<3.06))
499  n_reco_vtx_2++;
500  }
501  }
502 
503 
506 // for (l=0;l<phi1_pid.GetLength();++l) {
507 // RhoKinFitter kinfitter(phi1_pid[l]);
508 // // set it's mass constraint to phi mass
509 // kinfitter.AddMassConstraint(1.02);
510 // kinfitter.Fit();
511 // TCandidate phifit=kinfitter.GetFitted(phi1_pid[l]);
512 //
513 // phi1_massfit.Add(phifit);
514 // h_mphi_final_massfit->Fill(phifit.M());
515 // }
516 //
517 // etac_massfit.Combine(phi1_massfit,phi1_massfit);
518 //
519 // Double_t sigma_etac=0.0331, m_etac_pdg=2.98;
520 // Double_t chi2_m, best_chi2_m=1000, etac_m_best_mass;
521 // TCandidate *etac_m_best=0;
522 // for (j=0;j<etac_massfit.GetLength();++j)
523 // {
524 // metac=etac_massfit[j].M();
525 // chi2_m=pow(metac-m_etac_pdg,2)/pow(sigma_etac,2);
526 // if(chi2_m<best_chi2_m)
527 // {
528 // best_chi2_m=chi2_m;
529 // best_i=j;
530 // etac_m_best=&etac_massfit[j];
531 // etac_m_best_mass=etac_m_best->M();
532 // }
533 // }
534 //
535 // if (etac_m_best!=0)
536 // {
537 // double m_phi1=etac_m_best->Daughter(0)->M();
538 // double m_phi2=etac_m_best->Daughter(1)->M();
539 // if (((m_phi1>1.02-0.02)&&(m_phi1<1.02+0.02))&&
540 // ((m_phi2>1.02-0.02)&&(m_phi2<1.02+0.02)))
541 // {
542 // h_etac_phimassfit->Fill(etac_m_best_mass);
543 // }
544 // }
545 
546 
547  }
548 
549  std::cout<<"Number of reconstructed eta_c (4C) = "<<n_reco_4c<<std::endl;
550  std::cout<<"Number of reconstructed eta_c (Vertex fit)= "<<n_reco_vtx<<std::endl;
551  std::cout<<"Number of reconstructed eta_c (Vertex fit, preselection)= "<<n_reco_vtx_2<<std::endl;
552  n_etac_4c->SetBinContent(1,n_reco_4c);
553  n_etac_vtx->SetBinContent(1,n_reco_vtx);
554  n_etac_vtx_2->SetBinContent(1,n_reco_vtx_2);
555 
556  out->cd();
557  n_etac_4c->Write();
558  n_etac_vtx->Write();
559  n_etac_vtx_2->Write();
560  n_events->Write();
561  h_etac_nocut->Write();
562  h_etac_pid->Write();
563  h_etac_vtx->Write();
564  h_etac_vtx_2->Write();
565  h_etac_4c->Write();
566  h_etac_4c_refit->Write();
567  h_etac_phimass_4c->Write();
568  h_etac_phimass_vtx->Write();
569  h_etac_phimass_vtx_2->Write();
570  h_etac_phimassfit->Write();
571 
572  h_mphi_nocuts->Write();
573  h_mphi_pid->Write();
574  h_mphi_vtx->Write();
575  h_mphi_vtx_2->Write();
576  h_mphi_4c->Write();
577  h_mphi_final_4c->Write();
578  h_mphi_final_vtx->Write();
579  h_mphi_final_vtx_2->Write();
580  h_mphi_final_massfit->Write();
581 
582  nc->Write();
583 
584  h_chi2_4c->Write();
585  h_chi2b_4c->Write();
586  h_chi2_vtx->Write();
587  h_chi2b_vtx->Write();
588  h_chi2_mass->Write();
589  h_chi2_prefit->Write();
590  hvzpos->Write();
591  hvpos->Write();
592 
593  hvtxresX->Write();
594  hvtxresY->Write();
595  hvtxresZ->Write();
596 
597  h_theta_p->Write();
598  h_theta_p_nr->Write();
599  h_dp_p->Write();
600  h_dp_theta->Write();
601  h_dp->Write();
602  h_dp_low->Write();
603  h_dp_high->Write();
604 
605  h_acc_stt->Write();
606 
607  out->Save();
608 
609  timer.Stop();
610  Double_t rtime = timer.RealTime();
611  Double_t ctime = timer.CpuTime();
612  printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
613 
614  return 0;
615 }
mctrack
PndMCTrack * mctrack
Definition: analysis_point_hit_fwendcap.C:250
h_chi2_vtx
TH1F * h_chi2_vtx
Definition: plot_eta_c_stt.C:33
h_mphi_final_vtx_2
TH1F * h_mphi_final_vtx_2
Definition: plot_eta_c_stt.C:23
kMVD
Definition: production/tdrct/eta_c/track_check.C:4
h_mphi_final_massfit
TH1F * h_mphi_final_massfit
Definition: plot_eta_c_tpc.C:21
printf
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
h_etac_nocut
TH1F * h_etac_nocut
Definition: plot_eta_c_stt.C:7
hvtxresX
TH1F * hvtxresX
Definition: plot_eta_c_stt.C:38
kTOF
Definition: production/tdrct/eta_c/track_check.C:4
i
Int_t i
Definition: run_full.C:25
kHYP
Definition: production/tdrct/eta_c/track_check.C:4
tree
TTree * tree
Definition: plot_dirc.C:12
h_etac_4c
TH1F * h_etac_4c
Definition: plot_eta_c_stt.C:10
PndMCTrack::GetNPoints
Int_t GetNPoints(DetectorId detId) const
Definition: PndMCTrack.cxx:120
Rho4CFitter::FitConserveMasses
Bool_t FitConserveMasses()
Definition: Rho4CFitter.cxx:64
PndMCTrack::GetPdgCode
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
hvzpos
TH1F * hvzpos
Definition: plot_eta_c_stt.C:35
h_chi2_prefit
TH1F * h_chi2_prefit
Definition: plot_eta_c_stt.C:26
kLUMI
Definition: production/tdrct/eta_c/track_check.C:4
kTPC
Definition: production/tdrct/eta_c/track_check.C:4
PndMCTrack::GetMomentum
TVector3 GetMomentum() const
Definition: PndMCTrack.h:78
mc_array
TClonesArray * mc_array
Definition: anaLmdCluster.C:25
h_etac_phimass_4c
TH1F * h_etac_phimass_4c
Definition: plot_eta_c_stt.C:11
inFile
TString inFile
Definition: hit_dirc.C:8
h_mphi_nocuts
TH1F * h_mphi_nocuts
Definition: plot_eta_c_stt.C:14
h_mphi_final_4c
TH1F * h_mphi_final_4c
Definition: plot_eta_c_stt.C:18
n_etac_vtx
TH1F * n_etac_vtx
Definition: plot_eta_c_stt.C:43
PndPidCandidate::GetMcIndex
Int_t GetMcIndex() const
Definition: PndPidCandidate.h:83
h_mphi_4c
TH1F * h_mphi_4c
Definition: plot_eta_c_stt.C:17
kSTT
Definition: production/tdrct/eta_c/track_check.C:4
h_mphi_final_vtx
TH1F * h_mphi_final_vtx
Definition: plot_eta_c_stt.C:19
kDSK
Definition: production/tdrct/eta_c/track_check.C:4
h_etac_vtx
TH1F * h_etac_vtx
Definition: plot_eta_c_stt.C:9
inSimFile
TString inSimFile
Definition: run_reco_llbar_evtgen.C:8
h_etac_phimassfit
TH1F * h_etac_phimassfit
Definition: plot_eta_c_tpc.C:13
Double_t
Double_t
Definition: SimCompleteLinkDef.h:6
kDRC
Definition: production/tdrct/eta_c/track_check.C:4
h_chi2b_4c
TH1F * h_chi2b_4c
Definition: plot_eta_c_stt.C:31
DetectorId
DetectorId
Definition: production/tdrct/eta_c/track_check.C:2
h_chi2_4c
TH1F * h_chi2_4c
Definition: plot_eta_c_stt.C:30
h_mphi_vtx
TH1F * h_mphi_vtx
Definition: plot_eta_c_stt.C:16
timer
TStopwatch timer
Definition: hit_dirc.C:51
kMDT
Definition: production/tdrct/eta_c/track_check.C:4
h_etac_phimass_vtx
TH1F * h_etac_phimass_vtx
Definition: plot_eta_c_stt.C:12
kRPC
Definition: production/tdrct/eta_c/track_check.C:4
kHYPG
Definition: production/tdrct/eta_c/track_check.C:4
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: P4_F32vec4.h:47
out
TFile * out
Definition: reco_muo.C:20
n_events
TH1F * n_events
Definition: plot_eta_c_stt.C:44
h_etac_phimass_vtx_2
TH1F * h_etac_phimass_vtx_2
Definition: plot_eta_c_stt.C:21
run_ana_eta_c_stt
int run_ana_eta_c_stt(int nevts=0, bool usePID=true)
Definition: run_ana_eta_c_stt.C:5
ctime
Double_t ctime
Definition: hit_dirc.C:114
kEMC
Definition: production/tdrct/eta_c/track_check.C:4
h_etac_pid
TH1F * h_etac_pid
Definition: plot_eta_c_stt.C:8
p2
TPad * p2
Definition: hist-t7.C:117
h_mphi_vtx_2
TH1F * h_mphi_vtx_2
Definition: plot_eta_c_stt.C:22
n_etac_vtx_2
TH1F * n_etac_vtx_2
Definition: plot_eta_c_stt.C:25
hvtxresZ
TH1F * hvtxresZ
Definition: plot_eta_c_stt.C:40
p1
TPad * p1
Definition: hist-t7.C:116
RhoFitterBase::GetChi2
double GetChi2() const
Definition: RhoFitterBase.h:48
kGEM
Definition: production/tdrct/eta_c/track_check.C:4
h_chi2b_vtx
TH1F * h_chi2b_vtx
Definition: plot_eta_c_stt.C:34
PndMCTrack::GetStartVertex
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
PndMCTrack::GetMotherID
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
rtime
Double_t rtime
Definition: hit_dirc.C:113
PndPidCandidate::GetMomentum
TVector3 GetMomentum() const
Definition: PndPidCandidate.h:55
h_mphi_pid
TH1F * h_mphi_pid
Definition: plot_eta_c_stt.C:15
kFTS
Definition: production/tdrct/eta_c/track_check.C:4
n_etac_4c
TH1F * n_etac_4c
Definition: plot_eta_c_stt.C:42
hvtxresY
TH1F * hvtxresY
Definition: plot_eta_c_stt.C:39
hvpos
TH2F * hvpos
Definition: plot_eta_c_stt.C:36
nc
TH1F * nc
Definition: plot_eta_c_stt.C:28
kDCH
Definition: production/tdrct/eta_c/track_check.C:4
Generated on Sun Aug 18 2019 03:01:58 for FairRoot/PandaRoot by   doxygen 1.8.8