FairRoot/PandaRoot
Functions
test_check.C File Reference

Go to the source code of this file.

Functions

int test_check (Int_t nEntries=0)
 

Function Documentation

int test_check ( Int_t  nEntries = 0)

Definition at line 1 of file test_check.C.

References PndPidCandidate::GetCharge(), PndPidCandidate::GetEmcIndex(), PndPidCandidate::GetEmcRawEnergy(), PndPidCandidate::GetMcIndex(), PndPidCandidate::GetMomentum(), PndMCTrack::GetMotherID(), PndPidCandidate::GetMuoIndex(), PndPidCandidate::GetMuoIron(), PndPidCandidate::GetMuoModule(), PndPidCandidate::GetMuoMomentumIn(), PndPidProbability::GetMuonPidProb(), PndMCTrack::GetPdgCode(), hit, inFile, inSimFile, m2(), mc_array, mm, mup, out, rootlogon(), sqrt(), tree, and TString.

2 {
3  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
4  rootlogon();
5  TString inPidFile = "pid_sttcombi.root";
6  TString inSimFile = "points_sttcombi.root";
7 
8  TFile *inFile = TFile::Open(inPidFile,"READ");
9 
10  TTree *tree=(TTree *) inFile->Get("pndsim") ;
11  tree->AddFriend("pndsim",inSimFile);
12 
13  TClonesArray* cand_array=new TClonesArray("PndPidCandidate");
14  tree->SetBranchAddress("PidChargedCand", &cand_array);
15 
16  TClonesArray* mc_array=new TClonesArray("PndMCTrack");
17  tree->SetBranchAddress("MCTrack", &mc_array);
18 
19  TClonesArray* mdt_array=new TClonesArray("PndPidProbability");
20  tree->SetBranchAddress("PidAlgoMdtHardCuts", &mdt_array);
21 
22  TFile *out = TFile::Open("out_test2.root","RECREATE");
23  //TNtuple *nt = new TNtuple("nt","nt","evt:mc_p:mc_theta:mc_phi:mc_pid:p:theta:phi:mult");
24  TNtuple *ee = new TNtuple("ee","ee","p1:theta1:phi1:p2:theta2:phi2:pid1:emc1:emc_ene1:muo1:muo_mod1:muo_iron1:muo_mom1:muo_prob1:pid2:emc2:emc_ene2:muo2:muo_mod2:muo_iron2:muo_mom2:muo_prob2:mass_ee:mass_mu:mass_pi");
25  TNtuple *nt4 = new TNtuple("nt4","nt4","p1:theta1:phi1:p2:theta2:phi2:pid1:emc1:emc_ene1:muo1:muo_mod1:muo_iron1:muo_mom1:muo_prob1:pid2:emc2:emc_ene2:muo2:muo_mod2:muo_iron2:muo_mom2:muo_prob2:mass_ee:mass_mu:mass_pi:p3:theta3:phi3:p4:theta4:phi4:pid3:emc3:emc_ene3:muo3:muo_mod3:muo_iron3:muo_mom3:muo_prob3:pid4:emc4:emc_ene4:muo4:muo_mod4:muo_iron4:muo_mom4:muo_prob4:mass2_ee:mass2_mu:mass2_pi:X_e:X_mu");
26  //TNtuple *mm = new TNtuple("mm","mm","p3:theta3:phi1:p2:theta2:phi2:mass");
27  //TNtuple *pp = new TNtuple("pp","pp","p1:theta1:phi1:p2:theta2:phi2:mass:miss");
28  TH1F *m1=new TH1F("m1","m1",100,0,4);
29  TH1F *m2=new TH1F("m2","m2",100,0,4);
30  TH1F *m3=new TH1F("m3","m3",100,0,4);
31  TH1F *m4=new TH1F("m4","m4",100,0,6);
32  TH1F *m5=new TH1F("m5","m5",100,0,6);
33  TH1F *mc1=new TH1F("mc1","mc1",100,0,4);
34  TH1F *mc2=new TH1F("mc2","mc2",100,0,4);
35  TH1F *mc3=new TH1F("mc3","mc3",100,0,4);
36  TH1F *mc4=new TH1F("mc4","mc4",100,0,6);
37  TH1F *mc5=new TH1F("mc5","mc5",100,0,6);
38 
39  TH2F *ps1 = new TH2F("ps1","ps1",100,0,7,100,0,180);
40  TH2F *ps2 = new TH2F("ps2","ps2",100,0,7,100,0,1.5);
41  TH2F *ps3 = new TH2F("ps3","ps3",100,0,7,100,0,100);
42  TH2F *ps4 = new TH2F("ps4","ps4",100,0,7,100,0,100);
43  TNtuple *hit = new TNtuple("hit","hit","q:p:theta:phi:pid:emc:emc_ene:muo:muo_mod:muo_iron:muo_mom:muo_prob");
44 
45  TLorentzVector ini(0., 0., 6.98840096, 7.98937619);
46  if (nEntries==0) nEntries = tree->GetEntriesFast();
47  for (Int_t j=0; j< nEntries; j++){
48  tree->GetEntry(j);
49  //if (cand_array->GetEntriesFast()==0) continue;
50  cout << "processing event " << j << "\n";
51 
52  TLorentzVector ep[15], em[15], pip[15], pim[15], mup[15], mum[15], jpsi_e, jpsi_mu;
53  TLorentzVector mep[15], mem[15], mpip[15], mpim[15], mmup[15], mmum[15], mjpsi_e, mjpsi_mu;
54  Int_t ep_count = 0, em_count = 0, pip_count = 0, pim_count = 0, mup_count = 0, mum_count = 0;
55  Int_t ep_mcount = 0, em_mcount = 0, pip_mcount = 0, pim_mcount = 0, mup_mcount = 0, mum_mcount = 0;
56  for (Int_t pp=0; pp<cand_array->GetEntriesFast(); pp++)
57  {
58  PndPidCandidate *pidCand1 = (PndPidCandidate*)cand_array->At(pp);
59  if ((pidCand1->GetCharge()==1) && (pidCand1->GetMomentum().Mag()<8) )
60  {
61  PndPidProbability *mdtProb1 = (PndPidProbability*)mdt_array->At(pp);
62  if (pidCand1->GetMomentum().Mag()>8) continue;
63  Int_t mc_pid1 = 0;
64  if (pidCand1->GetMcIndex()>-1)
65  {
66  PndMCTrack* mcTrack = (PndMCTrack*)mc_array->At(pidCand1->GetMcIndex());
67  if (mcTrack->GetMotherID()==-1) mc_pid1 = mcTrack->GetPdgCode();
68  }
69  Float_t hit_ntuple1[] = {pidCand1->GetCharge(), pidCand1->GetMomentum().Mag(), pidCand1->GetMomentum().Theta()*TMath::RadToDeg(), pidCand1->GetMomentum().Phi()*TMath::RadToDeg(), mc_pid1,
70  pidCand1->GetEmcIndex(),pidCand1->GetEmcRawEnergy(),
71  pidCand1->GetMuoIndex(), pidCand1->GetMuoModule(), pidCand1->GetMuoIron(), pidCand1->GetMuoMomentumIn(), mdtProb1->GetMuonPidProb()
72  };
73  hit->Fill(hit_ntuple1);
74  for (Int_t mm=0; mm<cand_array->GetEntriesFast(); mm++)
75  {
76  PndPidCandidate *pidCand2 = (PndPidCandidate*)cand_array->At(mm);
77  if ( (pidCand2->GetCharge()==-1) && (pidCand2->GetMomentum().Mag()<8) )
78  {
79  PndPidProbability *mdtProb2 = (PndPidProbability*)mdt_array->At(mm);
80  if (pidCand2->GetMomentum().Mag()>8) continue;
81  Int_t mc_pid2 = 0;
82  if (pidCand2->GetMcIndex()>-1)
83  {
84  PndMCTrack* mcTrack = (PndMCTrack*)mc_array->At(pidCand2->GetMcIndex());
85  if (mcTrack->GetMotherID()==-1) mc_pid2 = mcTrack->GetPdgCode();
86  }
87  Float_t hit_ntuple2[] = {pidCand2->GetCharge(), pidCand2->GetMomentum().Mag(), pidCand2->GetMomentum().Theta()*TMath::RadToDeg(), pidCand2->GetMomentum().Phi()*TMath::RadToDeg(), mc_pid2,
88  pidCand2->GetEmcIndex(),pidCand2->GetEmcRawEnergy(),
89  pidCand2->GetMuoIndex(), pidCand2->GetMuoModule(), pidCand2->GetMuoIron(), pidCand2->GetMuoMomentumIn(), mdtProb2->GetMuonPidProb()
90  };
91  hit->Fill(hit_ntuple2);
92 
93  TLorentzVector lorep(pidCand1->GetMomentum(),sqrt(pidCand1->GetMomentum()*pidCand1->GetMomentum()+0.0005*0.0005));
94  TLorentzVector lorem(pidCand2->GetMomentum(),sqrt(pidCand2->GetMomentum()*pidCand2->GetMomentum()+0.0005*0.0005));
95  TLorentzVector lorpip(pidCand1->GetMomentum(),sqrt(pidCand1->GetMomentum()*pidCand1->GetMomentum()+0.140*0.140));
96  TLorentzVector lorpim(pidCand2->GetMomentum(),sqrt(pidCand2->GetMomentum()*pidCand2->GetMomentum()+0.140*0.140));
97  TLorentzVector lormup(pidCand1->GetMomentum(),sqrt(pidCand1->GetMomentum()*pidCand1->GetMomentum()+0.113*0.113));
98  TLorentzVector lormum(pidCand2->GetMomentum(),sqrt(pidCand2->GetMomentum()*pidCand2->GetMomentum()+0.113*0.113));
99  Float_t ee_ntuple[] = {
100  lorep.P(), lorep.Theta()*TMath::RadToDeg(), lorep.Phi()*TMath::RadToDeg(),
101  lorem.P(), lorem.Theta()*TMath::RadToDeg(), lorem.Phi()*TMath::RadToDeg(),
102  mc_pid1, pidCand1->GetEmcIndex(),pidCand1->GetEmcRawEnergy(),
103  pidCand1->GetMuoIndex(), pidCand1->GetMuoModule(), pidCand1->GetMuoIron(), pidCand1->GetMuoMomentumIn(), mdtProb1->GetMuonPidProb(),
104  mc_pid2, pidCand2->GetEmcIndex(),pidCand2->GetEmcRawEnergy(),
105  pidCand2->GetMuoIndex(), pidCand2->GetMuoModule(), pidCand2->GetMuoIron(), pidCand2->GetMuoMomentumIn(), mdtProb2->GetMuonPidProb(),
106  (lorep+lorem).M(), (lormup+lormum).M(), (ini-(lorpip+lorpim)).M()
107  };
108  ee->Fill(ee_ntuple);
109 
110  for (Int_t pp2=0; pp2<cand_array->GetEntriesFast(); pp2++)
111  {
112  PndPidCandidate *pidCand3 = (PndPidCandidate*)cand_array->At(pp2);
113  if ((pidCand3->GetCharge()==1) && (pidCand3->GetMomentum().Mag()<8) )
114  {
115  PndPidProbability *mdtProb3 = (PndPidProbability*)mdt_array->At(pp2);
116  if (pidCand3->GetMomentum().Mag()>8) continue;
117  Int_t mc_pid3 = 0;
118  if (pidCand3->GetMcIndex()>-1)
119  {
120  PndMCTrack* mcTrack = (PndMCTrack*)mc_array->At(pidCand3->GetMcIndex());
121  if (mcTrack->GetMotherID()==-1) mc_pid3 = mcTrack->GetPdgCode();
122  }
123 
124  for (Int_t mm2=0; mm2<cand_array->GetEntriesFast(); mm2++)
125  {
126  PndPidCandidate *pidCand4 = (PndPidCandidate*)cand_array->At(mm2);
127  if ( (pidCand4->GetCharge()==-1) && (pidCand4->GetMomentum().Mag()<8) )
128  {
129  PndPidProbability *mdtProb4 = (PndPidProbability*)mdt_array->At(mm2);
130  if (pidCand4->GetMomentum().Mag()>8) continue;
131  Int_t mc_pid4 = 0;
132  if (pidCand4->GetMcIndex()>-1)
133  {
134  PndMCTrack* mcTrack = (PndMCTrack*)mc_array->At(pidCand4->GetMcIndex());
135  if (mcTrack->GetMotherID()==-1) mc_pid4 = mcTrack->GetPdgCode();
136  }
137 
138 
139  TLorentzVector lor2ep(pidCand3->GetMomentum(),sqrt(pidCand3->GetMomentum()*pidCand3->GetMomentum()+0.0005*0.0005));
140  TLorentzVector lor2em(pidCand4->GetMomentum(),sqrt(pidCand4->GetMomentum()*pidCand4->GetMomentum()+0.0005*0.0005));
141  TLorentzVector lor2pip(pidCand3->GetMomentum(),sqrt(pidCand3->GetMomentum()*pidCand3->GetMomentum()+0.140*0.140));
142  TLorentzVector lor2pim(pidCand4->GetMomentum(),sqrt(pidCand4->GetMomentum()*pidCand4->GetMomentum()+0.140*0.140));
143  TLorentzVector lor2mup(pidCand3->GetMomentum(),sqrt(pidCand3->GetMomentum()*pidCand3->GetMomentum()+0.113*0.113));
144  TLorentzVector lor2mum(pidCand4->GetMomentum(),sqrt(pidCand4->GetMomentum()*pidCand4->GetMomentum()+0.113*0.113));
145  Float_t p4_ntuple[] = {
146  lorep.P(), lorep.Theta()*TMath::RadToDeg(), lorep.Phi()*TMath::RadToDeg(),
147  lorem.P(), lorem.Theta()*TMath::RadToDeg(), lorem.Phi()*TMath::RadToDeg(),
148  mc_pid1, pidCand1->GetEmcIndex(),pidCand1->GetEmcRawEnergy(),
149  pidCand1->GetMuoIndex(), pidCand1->GetMuoModule(), pidCand1->GetMuoIron(), pidCand1->GetMuoMomentumIn(), mdtProb1->GetMuonPidProb(),
150  mc_pid2, pidCand2->GetEmcIndex(),pidCand2->GetEmcRawEnergy(),
151  pidCand2->GetMuoIndex(), pidCand2->GetMuoModule(), pidCand2->GetMuoIron(), pidCand2->GetMuoMomentumIn(), mdtProb2->GetMuonPidProb(),
152  (lorep+lorem).M(), (lormup+lormum).M(), (ini-(lorpip+lorpim)).M() ,
153  lor2ep.P(), lor2ep.Theta()*TMath::RadToDeg(), lor2ep.Phi()*TMath::RadToDeg(),
154  lor2em.P(), lor2em.Theta()*TMath::RadToDeg(), lor2em.Phi()*TMath::RadToDeg(),
155  mc_pid3, pidCand3->GetEmcIndex(),pidCand3->GetEmcRawEnergy(),
156  pidCand3->GetMuoIndex(), pidCand3->GetMuoModule(), pidCand3->GetMuoIron(), pidCand3->GetMuoMomentumIn(), mdtProb3->GetMuonPidProb(),
157  mc_pid4, pidCand4->GetEmcIndex(),pidCand4->GetEmcRawEnergy(),
158  pidCand4->GetMuoIndex(), pidCand4->GetMuoModule(), pidCand4->GetMuoIron(), pidCand4->GetMuoMomentumIn(), mdtProb4->GetMuonPidProb(),
159  (lor2ep+lor2em).M(), (lor2mup+lor2mum).M(), (ini-(lor2pip+lor2pim)).M(),(lorpip+lorpim+lor2ep+lor2em).M(),(lorpip+lorpim+lor2mup+lor2mum).M()
160  };
161  nt4->Fill(p4_ntuple);
162  }
163  } // end of mm2
164  }
165  } // end of pp2
166  }
167  } // end of mm
168  }
169  } // end of pp
170 
171 
172 
173 
174 
175  } // end of event loop
176 
177 
178  out->cd();
179  ps1->Write();
180  ps2->Write();
181  ps3->Write();
182  ps4->Write();
183  m1->Write();
184  m2->Write();
185  m3->Write();
186  m4->Write();
187  m5->Write();
188  mc1->Write();
189  mc2->Write();
190  mc3->Write();
191  mc4->Write();
192  mc5->Write();
193  hit->Write();
194  ee->Write();
195  nt4->Write();
196  out->Save();
197 
198  return 0;
199 }
Int_t GetCharge() const
Float_t GetEmcRawEnergy() const
friend F32vec4 sqrt(const F32vec4 &a)
Definition: P4_F32vec4.h:29
TTree * tree
Definition: plot_dirc.C:12
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
Float_t GetMuoIron() const
TClonesArray * mc_array
Definition: anaLmdCluster.C:25
Float_t GetMuoMomentumIn() const
Int_t GetMuoIndex() const
Int_t GetMcIndex() const
Int_t GetEmcIndex() const
TString inFile
Definition: hit_dirc.C:8
TString m2(TString pts, TString exts="e px py pz")
Definition: invexp.C:117
TString inSimFile
Double_t GetMuonPidProb(PndPidProbability *flux=NULL) const
TFile * out
Definition: reco_muo.C:20
Double_t const mm
Int_t GetMuoModule() const
PndSdsMCPoint * hit
Definition: anasim.C:70
Int_t GetMotherID() const
Definition: PndMCTrack.h:74
TVector3 GetMomentum() const
static const double mup
Definition: mzparameters.h:21