FairRoot/PandaRoot
Functions
run_ana_mertens_evt7.C File Reference
#include "TFile.h"
#include "TH1D.h"
#include "TCanvas.h"
#include "TStopwatch.h"
#include "TROOT.h"
#include "TSystem.h"
#include "TTree.h"
#include "TString.h"

Go to the source code of this file.

Functions

void SaveAndUpdateHisto (TH1D *currenthisto, TFile &storagefile)
 
int run_ana_mertens_evt7 (TString fname, TString simfname, int nevts=0, TString outfilename="psiana.root")
 

Function Documentation

int run_ana_mertens_evt7 ( TString  fname,
TString  simfname,
int  nevts = 0,
TString  outfilename = "psiana.root" 
)

Definition at line 19 of file run_ana_mertens_evt7.C.

References ctime, Double_t, RhoFitterBase::Fit(), PndMCTrack::GetPdgCode(), RhoVtxPoca::GetPocaVtx(), PndMCTrack::GetStartVertex(), i, Mass, mctrack, printf(), rootlogon(), rtime, SaveAndUpdateHisto(), timer, x, y, and z.

20 {
21 
22 TStopwatch timer;
23 timer.Start();
24 
25 gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");rootlogon();
26 
27 PndEventReader evr(fname);
28 if (nevts==0) nevts=evr.GetEntries();
29 
30 // Event Statistics Histos
31 TH1D* hkpperevent = new TH1D("hkpperevent","K+ per Event", 20, 0, 20);
32 TH1D* hkmperevent = new TH1D("hkmperevent","K- per Event", 20, 0, 20);
33 TH1D* hpipperevent = new TH1D("hpipperevent","Pi+ per Event", 20, 0, 20);
34 TH1D* hpimperevent = new TH1D("hpimperevent","Pi- per Event", 20, 0, 20);
35 TH1D* hniceevents = new TH1D("hniceevents", "Nice Events", 10, 0, 10);
36 
37 // D Histos
38 TH1D* hdprawmass = new TH1D("hdprawmass","D+ Raw Mass",2000,0,4);
39 TH1D* hdmrawmass = new TH1D("hdmrawmass","D- Raw Mass",2000,0,4);
40 TH1D* hdpmselmass = new TH1D("hdpmselmass","D+ Selected Mass",2000,0,4);
41 TH1D* hdmmselmass = new TH1D("hdmmselmass","D- Selected Mass",2000,0,4);
42 TH1D* hdpvtxacceptmass = new TH1D("hdpvtxacceptmass","D+ Vertex Fit Accepted Mass",2000,0,4);
43 TH1D* hdpvtxrejectmass = new TH1D("hdpvtxrejectmass","D+ Vertex Fit Rejected Mass",2000,0,4);
44 hdpvtxrejectmass->SetLineColor(kRed);
45 TH1D* hdmvtxacceptmass = new TH1D("hdmvtxacceptmass","D- Vertex Fit Accepted Mass",2000,0,4);
46 TH1D* hdmvtxrejectmass = new TH1D("hdmvtxrejectmass","D- Vertex Fit Rejected Mass",2000,0,4);
47 hdmvtxrejectmass->SetLineColor(kRed);
48 TH1D* hdpvtxchi2 = new TH1D("hdpvtxchi2","D+ Vertex Fit Chi2",400,0,200);
49 TH1D* hdmvtxchi2 = new TH1D("hdmvtxchi2","D- Vertex Fit Chi2",400,0,200);
50 
51 // Raw Psi Histos
52 TH1D* hrawpsirawmass = new TH1D("hrawpsirawmass","Psi Raw Mass (Raw D)",4000,0,8);
53 TH1D* hrawpsirawpt = new TH1D("hrawpsirawpt","Psi Raw Pt (Raw D)",4000,0,8);
54 TH1D* hrawpsirawpz = new TH1D("hrawpsirawpz","Psi Raw Pz (Raw D)",6000,0,12);
55 TH1D* hrawpsirawdpt = new TH1D("hrawpsirawdpt","Psi Raw D Pt (Raw D)",4000,0,8);
56 TH1D* hrawpsirawdpz = new TH1D("hrawpsirawdpz","Psi Raw D Pz (Raw D)",6000,0,12);
57 
58 TH1D* hrawpsimselmass = new TH1D("hrawpsimselmass","Psi Raw Mass (Mass Selected D)",4000,0,8);
59 TH1D* hrawpsimselpt = new TH1D("hrawpsimselpt","Psi Raw Pt (Mass Selected D)",4000,0,8);
60 TH1D* hrawpsimselpz = new TH1D("hrawpsimselpz","Psi Raw Pz (Mass Selected D)",6000,0,12);
61 TH1D* hrawpsimseldpt = new TH1D("hrawpsimseldpt","Psi Raw D Pt (Mass Selected D)",4000,0,8);
62 TH1D* hrawpsimseldpz = new TH1D("hrawpsimseldpz","Psi Raw D Pz (Mass Selected D)",6000,0,12);
63 
64 TH1D* hrawpsivtxmass = new TH1D("hrawpsivtxmass","Psi Raw Mass (Vertex Fit D)",4000,0,8);
65 TH1D* hrawpsivtxpt = new TH1D("hrawpsivtxpt","Psi Raw Pt (Vertex Fit D)",4000,0,8);
66 TH1D* hrawpsivtxpz = new TH1D("hrawpsivtxpz","Psi Raw Pz (Vertex Fit D)",6000,0,12);
67 TH1D* hrawpsivtxdpt = new TH1D("hrawpsivtxdpt","Psi Raw D Pt (Vertex Fit D)",4000,0,8);
68 TH1D* hrawpsivtxdpz = new TH1D("hrawpsivtxdpz","Psi Raw D Pz (Vertex Fit D)",6000,0,12);
69 
70 TH1D* hdpvertex = new TH1D("hdpvertex","D+ Vertex Resolution",800,0,8);
71 TH1D* hdpvertexfit = new TH1D("hdpvertexfit","D+ Vertex Resolution (After Vertex Fit)",800,0,8);
72 TH1D* hdmvertex = new TH1D("hdmvertex","D- Vertex Resolution",800,0,8);
73 TH1D* hdmvertexfit = new TH1D("hdmvertexfit","D- Vertex Resolution (After Vertex Fit)",800,0,8);
74 TH1D* hfinalpsimass = new TH1D("hfinalpsimass","Psi Mass after all cuts (Vertex Fit D)",800,0,8);
75 
76 TH1D* hdpvertexx = new TH1D("hdpvertexx","D+ Vertex Resolution X",800,-4,4);
77 TH1D* hdpvertexxfit = new TH1D("hdpvertexxfit","D+ Vertex Resolution X (After Vertex Fit)",800,-4,4);
78 TH1D* hdmvertexx = new TH1D("hdmvertexx","D- Vertex Resolution X",800,-4,4);
79 TH1D* hdmvertexxfit = new TH1D("hdmvertexxfit","D- Vertex Resolution X (After Vertex Fit)",800,-4,4);
80 
81 TH1D* hdpvertexy = new TH1D("hdpvertexy","D+ Vertex Resolution Y",800,-4,4);
82 TH1D* hdpvertexyfit = new TH1D("hdpvertexyfit","D+ Vertex Resolution Y (After Vertex Fit)",800,-4,4);
83 TH1D* hdmvertexy = new TH1D("hdmvertexy","D- Vertex Resolution Y",800,-4,4);
84 TH1D* hdmvertexyfit = new TH1D("hdmvertexyfit","D- Vertex Resolution Y (After Vertex Fit)",800,-4,4);
85 
86 TH1D* hdpvertexz = new TH1D("hdpvertexz","D+ Vertex Resolution Z",800,-4,4);
87 TH1D* hdpvertexzfit = new TH1D("hdpvertexzfit","D+ Vertex Resolution Z (After Vertex Fit)",800,-4,4);
88 TH1D* hdmvertexz = new TH1D("hdmvertexz","D- Vertex Resolution Z",800,-4,4);
89 TH1D* hdmvertexzfit = new TH1D("hdmvertexzfit","D- Vertex Resolution Z (After Vertex Fit)",800,-4,4);
90 
91 TH1D* hdpvertexxmc = new TH1D("hdpvertexxmc","D+ MC Vertex X",800,-4,4);
92 TH1D* hdpvertexymc = new TH1D("hdpvertexymc","D+ MC Vertex Y",800,-4,4);
93 TH1D* hdpvertexzmc = new TH1D("hdpvertexzmc","D+ MC Vertex Z",800,-4,4);
94 TH1D* hdmvertexxmc = new TH1D("hdmvertexxmc","D- MC Vertex X",800,-4,4);
95 TH1D* hdmvertexymc = new TH1D("hdmvertexymc","D- MC Vertex Y",800,-4,4);
96 TH1D* hdmvertexzmc = new TH1D("hdmvertexzmc","D- MC Vertex Z",800,-4,4);
97 
98 TH1D* hdpvertexxpoca = new TH1D("hdpvertexxpoca","D+ Poca Vertex X",800,-4,4);
99 TH1D* hdpvertexypoca = new TH1D("hdpvertexypoca","D+ Poca Vertex Y",800,-4,4);
100 TH1D* hdpvertexzpoca = new TH1D("hdpvertexzpoca","D+ Poca Vertex Z",800,-4,4);
101 TH1D* hdmvertexxpoca = new TH1D("hdmvertexxpoca","D- Poca Vertex X",800,-4,4);
102 TH1D* hdmvertexypoca = new TH1D("hdmvertexypoca","D- Poca Vertex Y",800,-4,4);
103 TH1D* hdmvertexzpoca = new TH1D("hdmvertexzpoca","D- Poca Vertex Z",800,-4,4);
104 
105 TH1D* hdpvertexxpocareso = new TH1D("hdpvertexxpocareso","D+ Poca Vertex Resolution X",800,-4,4);
106 TH1D* hdpvertexypocareso = new TH1D("hdpvertexypocareso","D+ Poca Vertex Resolution Y",800,-4,4);
107 TH1D* hdpvertexzpocareso = new TH1D("hdpvertexzpocareso","D+ Poca Vertex Resolution Z",800,-4,4);
108 TH1D* hdmvertexxpocareso = new TH1D("hdmvertexxpocareso","D- Poca Vertex Resolution X",800,-4,4);
109 TH1D* hdmvertexypocareso = new TH1D("hdmvertexypocareso","D- Poca Vertex Resolution Y",800,-4,4);
110 TH1D* hdmvertexzpocareso = new TH1D("hdmvertexzpocareso","D- Poca Vertex Resolution Z",800,-4,4);
111 
112 TH1D* hdpvertexxreco = new TH1D("hdpvertexxreco","D+ Reco Vertex X",800,-4,4);
113 TH1D* hdpvertexyreco = new TH1D("hdpvertexyreco","D+ Reco Vertex Y",800,-4,4);
114 TH1D* hdpvertexzreco = new TH1D("hdpvertexzreco","D+ Reco Vertex Z",800,-4,4);
115 TH1D* hdmvertexxreco = new TH1D("hdmvertexxreco","D- Reco Vertex X",800,-4,4);
116 TH1D* hdmvertexyreco = new TH1D("hdmvertexyreco","D- Reco Vertex Y",800,-4,4);
117 TH1D* hdmvertexzreco = new TH1D("hdmvertexzreco","D- Reco Vertex Z",800,-4,4);
118 
119 TH1D* hdpdecaylength = new TH1D("hdpdecaylength","D+ Decay Length",8000,0,8);
120 TH1D* hdmdecaylength = new TH1D("hdmdecaylength","D- Decay Length",8000,0,8);
121 
122 // Mass Selected Psi Histos
123 TH1D* hmselpsirawmass = new TH1D("hmselpsirawmass","Psi Mass Selected Mass (Raw D)",4000,0,8);
124 TH1D* hmselpsirawpt = new TH1D("hmselpsirawpt","Psi Mass Selected Pt (Raw D)",4000,0,8);
125 TH1D* hmselpsirawpz = new TH1D("hmselpsirawpz","Psi Mass Selected Pz (Raw D)",6000,0,12);
126 TH1D* hmselpsirawdpt = new TH1D("hmselpsirawdpt","Psi Mass Selected D Pt (Raw D)",4000,0,8);
127 TH1D* hmselpsirawdpz = new TH1D("hmselpsirawdpz","Psi Mass Selected D Pz (Raw D)",6000,0,12);
128 hmselpsirawmass->SetLineColor(kBlue);
129 hmselpsirawpt->SetLineColor(kBlue);
130 hmselpsirawpz->SetLineColor(kBlue);
131 hmselpsirawdpt->SetLineColor(kBlue);
132 hmselpsirawdpz->SetLineColor(kBlue);
133 
134 TH1D* hmselpsimselmass = new TH1D("hmselpsimselmass","Psi Mass Selected Mass (Mass Selected D)",4000,0,8);
135 TH1D* hmselpsimselpt = new TH1D("hmselpsimselpt","Psi Mass Selected Pt (Mass Selected D)",4000,0,8);
136 TH1D* hmselpsimselpz = new TH1D("hmselpsimselpz","Psi Mass Selected Pz (Mass Selected D)",6000,0,12);
137 TH1D* hmselpsimseldpt = new TH1D("hmselpsimseldpt","Psi Mass Selected D Pt (Mass Selected D)",4000,0,8);
138 TH1D* hmselpsimseldpz = new TH1D("hmselpsimseldpz","Psi Mass Selected D Pz (Mass Selected D)",6000,0,12);
139 hmselpsimselmass->SetLineColor(kBlue);
140 hmselpsimselpt->SetLineColor(kBlue);
141 hmselpsimselpz->SetLineColor(kBlue);
142 hmselpsimseldpt->SetLineColor(kBlue);
143 hmselpsimseldpz->SetLineColor(kBlue);
144 
145 TH1D* hmselpsivtxmass = new TH1D("hmselpsivtxmass","Psi Mass Selected Mass (Vertex Fit D)",4000,0,8);
146 TH1D* hmselpsivtxpt = new TH1D("hmselpsivtxpt","Psi Mass Selected Pt (Vertex Fit D)",4000,0,8);
147 TH1D* hmselpsivtxpz = new TH1D("hmselpsivtxpz","Psi Mass Selected Pz (Vertex Fit D)",6000,0,12);
148 TH1D* hmselpsivtxdpt = new TH1D("hmselpsivtxdpt","Psi Mass Selected D Pt (Vertex Fit D)",4000,0,8);
149 TH1D* hmselpsivtxdpz = new TH1D("hmselpsivtxdpz","Psi Mass Selected D Pz (Vertex Fit D)",6000,0,12);
150 hmselpsivtxmass->SetLineColor(kBlue);
151 hmselpsivtxpt->SetLineColor(kBlue);
152 hmselpsivtxpz->SetLineColor(kBlue);
153 hmselpsivtxdpt->SetLineColor(kBlue);
154 hmselpsivtxdpz->SetLineColor(kBlue);
155 
156 // D pt Selected Psi Histos
157 TH1D* hdptselpsirawmass = new TH1D("hdptselpsirawmass","Psi pt (D) Selected Mass (Raw D)",4000,0,8);
158 TH1D* hdptselpsirawpt = new TH1D("hdptselpsirawpt","Psi pt (D) Pt (Raw D)",4000,0,8);
159 TH1D* hdptselpsirawpz = new TH1D("hdptselpsirawpz","Psi pt (D) Selected Pz (Raw D)",6000,0,12);
160 TH1D* hdptselpsirawdpt = new TH1D("hdptselpsirawdpt","Psi pt (D) Selected D Pt (Raw D)",4000,0,8);
161 TH1D* hdptselpsirawdpz = new TH1D("hdptselpsirawdpz","Psi pt (D) Selected D Pz (Raw D)",6000,0,12);
162 
163 TH1D* hdptselpsimselmass = new TH1D("hdptselpsimselmass","Psi pt (D) Selected Mass (Mass Selected D)",4000,0,8);
164 TH1D* hdptselpsimselpt = new TH1D("hdptselpsimselpt","Psi pt (D) Selected Pt (Mass Selected D)",4000,0,8);
165 TH1D* hdptselpsimselpz = new TH1D("hdptselpsimselpz","Psi pt (D) Selected Pz (Mass Selected D)",6000,0,12);
166 TH1D* hdptselpsimseldpt = new TH1D("hdptselpsimseldpt","Psi pt (D) Selected D Pt (Mass Selected D)",4000,0,8);
167 TH1D* hdptselpsimseldpz = new TH1D("hdptselpsimseldpz","Psi pt (D) Selected D Pz (Mass Selected D)",6000,0,12);
168 
169 TH1D* hdptselpsivtxmass = new TH1D("hdptselpsivtxmass","Psi pt (D) Selected Mass (Vertex Fit D)",4000,0,8);
170 TH1D* hdptselpsivtxpt = new TH1D("hdptselpsivtxpt","Psi pt (D) Selected Pt (Vertex Fit D)",4000,0,8);
171 TH1D* hdptselpsivtxpz = new TH1D("hdptselpsivtxpz","Psi pt (D) Selected Pz (Vertex Fit D)",6000,0,12);
172 TH1D* hdptselpsivtxdpt = new TH1D("hdptselpsivtxdpt","Psi pt (D) Selected D Pt (Vertex Fit D)",4000,0,8);
173 TH1D* hdptselpsivtxdpz = new TH1D("hdptselpsivtxdpz","Psi pt (D) Selected D Pz (Vertex Fit D)",6000,0,12);
174 
175 // Psi pt Selected Psi Histos
176 TH1D* hpptselpsirawmass = new TH1D("hpptselpsirawmass","Psi pt (Psi) Selected Mass (Raw D)",4000,0,8);
177 TH1D* hpptselpsirawpt = new TH1D("hpptselpsirawpt","Psi pt (Psi) Pt (Raw D)",4000,0,8);
178 TH1D* hpptselpsirawpz = new TH1D("hpptselpsirawpz","Psi pt (Psi) Selected Pz (Raw D)",6000,0,12);
179 TH1D* hpptselpsirawdpt = new TH1D("hpptselpsirawdpt","Psi pt (Psi) Selected D Pt (Raw D)",4000,0,8);
180 TH1D* hpptselpsirawdpz = new TH1D("hpptselpsirawdpz","Psi pt (Psi) Selected D Pz (Raw D)",6000,0,12);
181 
182 TH1D* hpptselpsimselmass = new TH1D("hpptselpsimselmass","Psi pt (Psi) Selected Mass (Mass Selected D)",4000,0,8);
183 TH1D* hpptselpsimselpt = new TH1D("hpptselpsimselpt","Psi pt (Psi) Selected Pt (Mass Selected D)",4000,0,8);
184 TH1D* hpptselpsimselpz = new TH1D("hpptselpsimselpz","Psi pt (Psi) Selected Pz (Mass Selected D)",6000,0,12);
185 TH1D* hpptselpsimseldpt = new TH1D("hpptselpsimseldpt","Psi pt (Psi) Selected D Pt (Mass Selected D)",4000,0,8);
186 TH1D* hpptselpsimseldpz = new TH1D("hpptselpsimseldpz","Psi pt (Psi) Selected D Pz (Mass Selected D)",6000,0,12);
187 
188 TH1D* hpptselpsivtxmass = new TH1D("hpptselpsivtxmass","Psi pt (Psi) Selected Mass (Vertex Fit D)",4000,0,8);
189 TH1D* hpptselpsivtxpt = new TH1D("hpptselpsivtxpt","Psi pt (Psi) Selected Pt (Vertex Fit D)",4000,0,8);
190 TH1D* hpptselpsivtxpz = new TH1D("hpptselpsivtxpz","Psi pt (Psi) Selected Pz (Vertex Fit D)",6000,0,12);
191 TH1D* hpptselpsivtxdpt = new TH1D("hpptselpsivtxdpt","Psi pt (Psi) Selected D Pt (Vertex Fit D)",4000,0,8);
192 TH1D* hpptselpsivtxdpz = new TH1D("hpptselpsivtxdpz","Psi pt (Psi) Selected D Pz (Vertex Fit D)",6000,0,12);
193 
194 // D and Psi pt Selected Psi Histos
195 TH1D* hptselpsirawmass = new TH1D("hptselpsirawmass","Psi pt (Both) Selected Mass (Raw D)",4000,0,8);
196 TH1D* hptselpsirawpt = new TH1D("hptselpsirawpt","Psi pt (Both) Pt (Raw D)",4000,0,8);
197 TH1D* hptselpsirawpz = new TH1D("hptselpsirawpz","Psi pt (Both) Selected Pz (Raw D)",6000,0,12);
198 TH1D* hptselpsirawdpt = new TH1D("hptselpsirawdpt","Psi pt (Both) Selected D Pt (Raw D)",4000,0,8);
199 TH1D* hptselpsirawdpz = new TH1D("hptselpsirawdpz","Psi pt (Both) Selected D Pz (Raw D)",6000,0,12);
200 
201 TH1D* hptselpsimselmass = new TH1D("hptselpsimselmass","Psi pt (Both) Selected Mass (Mass Selected D)",4000,0,8);
202 TH1D* hptselpsimselpt = new TH1D("hptselpsimselpt","Psi pt (Both) Selected Pt (Mass Selected D)",4000,0,8);
203 TH1D* hptselpsimselpz = new TH1D("hptselpsimselpz","Psi pt (Both) Selected Pz (Mass Selected D)",6000,0,12);
204 TH1D* hptselpsimseldpt = new TH1D("hptselpsimseldpt","Psi pt (Both) Selected D Pt (Mass Selected D)",4000,0,8);
205 TH1D* hptselpsimseldpz = new TH1D("hptselpsimseldpz","Psi pt (Both) Selected D Pz (Mass Selected D)",6000,0,12);
206 
207 TH1D* hptselpsivtxmass = new TH1D("hptselpsivtxmass","Psi pt (Both) Selected Mass (Vertex Fit D)",4000,0,8);
208 TH1D* hptselpsivtxpt = new TH1D("hptselpsivtxpt","Psi pt (Both) Selected Pt (Vertex Fit D)",4000,0,8);
209 TH1D* hptselpsivtxpz = new TH1D("hptselpsivtxpz","Psi pt (Both) Selected Pz (Vertex Fit D)",6000,0,12);
210 TH1D* hptselpsivtxdpt = new TH1D("hptselpsivtxdpt","Psi pt (Both) Selected D Pt (Vertex Fit D)",4000,0,8);
211 TH1D* hptselpsivtxdpz = new TH1D("hptselpsivtxdpz","Psi pt (Both) Selected D Pz (Vertex Fit D)",6000,0,12);
212 
213 // Vertex Fit Psi Histos (only taken for preselected Psi)
214 TH1D* hvtxpsichi2 = new TH1D("hvtxpsichi2","Psi Vertex Fit Chi2",400,0,200);
215 TH1D* hvtxpsiacceptmass = new TH1D("hvtxpsiacceptmass","Psi Vertex Fit Accepted Mass",4000,0,8);
216 TH1D* hvtxpsiacceptpt = new TH1D("hvtxpsiacceptpt","Psi Vertex Fit Accepted Pt",4000,0,8);
217 TH1D* hvtxpsiacceptpz = new TH1D("hvtxpsiacceptpz","Psi Vertex Fit Accepted Pz",6000,0,12);
218 TH1D* hvtxpsiacceptdpt = new TH1D("hvtxpsiacceptdpt","Psi Vertex Fit Accepted D Pt",4000,0,8);
219 TH1D* hvtxpsiacceptdpz = new TH1D("hvtxpsiacceptdpz","Psi Vertex Fit Accepted D Pz",6000,0,12);
220 TH1D* hvtxpsirejectmass = new TH1D("hvtxpsirejectmass","Psi Vertex Fit Rejected Mass",4000,0,8);
221 TH1D* hvtxpsirejectpt = new TH1D("hvtxpsirejectpt","Psi Vertex Fit Rejected Pt",4000,0,8);
222 TH1D* hvtxpsirejectpz = new TH1D("hvtxpsirejectpz","Psi Vertex Fit Rejected Pz",6000,0,12);
223 TH1D* hvtxpsirejectdpt = new TH1D("hvtxpsirejectdpt","Psi Vertex Fit Rejected D Pt",4000,0,8);
224 TH1D* hvtxpsirejectdpz = new TH1D("hvtxpsirejectdpz","Psi Vertex Fit Rejected D Pz",6000,0,12);
225 hvtxpsirejectmass->SetLineColor(kRed);
226 hvtxpsirejectpt->SetLineColor(kRed);
227 hvtxpsirejectpz->SetLineColor(kRed);
228 hvtxpsirejectdpt->SetLineColor(kRed);
229 hvtxpsirejectdpz->SetLineColor(kRed);
230 
231 TH1D* hptselvtxpsimass = new TH1D("hptselvtxpsimass","Psi Vertex Fit pt (Both) Selected Mass",4000,0,8);
232 TH1D* hptselvtxpsipt = new TH1D("hptselvtxpsipt","Psi Vertex Fit pt (Both) Selected Pt",4000,0,8);
233 TH1D* hptselvtxpsipz = new TH1D("hptselvtxpsipz","Psi Vertex Fit pt (Both) Selected Pz",6000,0,12);
234 TH1D* hptselvtxpsidpt = new TH1D("hptselvtxpsidpt","Psi Vertex Fit pt (Both) Selected D Pt",4000,0,8);
235 TH1D* hptselvtxpsidpz = new TH1D("hptselvtxpsidpz","Psi Vertex Fit pt (Both) Selected D Pz",6000,0,12);
236 
237 TH1D* dpmass = new TH1D("dpmass","D+ Raw",2000,0,4);
238 TH1D* dmmass = new TH1D("dmmass","D- Raw",2000,0,4);
239 TH1D* dp2mass = new TH1D("dp2mass","D+ after Selection",2000,0,4);
240 TH1D* dm2mass = new TH1D("dm2mass","D- after Selection",2000,0,4);
241 TH1D* dp3mass = new TH1D("dp3mass","D+ after Vertex Fit",2000,0,4);
242 TH1D* dm3mass = new TH1D("dm3mass","D- after Vertex Fit",2000,0,4);
243 TH1D* dp3massrej = new TH1D("dp3massrej","D+ after Vertex Fit",2000,0,4);
244 dp3massrej->SetLineColor(kRed);
245 TH1D* dm3massrej = new TH1D("dm3massrej","D- after Vertex Fit",2000,0,4);
246 dm3massrej->SetLineColor(kRed);
247 
248 TH1D* dptfit = new TH1D("dptfit","Pt Difference",2000,-2,2);
249 TH1D* dptsel = new TH1D("dptsel","Pt Difference",2000,-2,2);
250 TH1D* psipt = new TH1D("psipt","Psi Pt",2000,-2,2);
251 
252 TH1D* dpchi2 = new TH1D("dpchi2","D+ Chi2",400,0,200);
253 TH1D* dmchi2 = new TH1D("dmchi2","D- Chi2",400,0,200);
254 TH1D* psichi2 = new TH1D("psichi2","Psi Chi2",400,0,200);
255 
256 TH1D* psimass = new TH1D("psimass","Psi Raw",4000,0,8);
257 TH1D* psimass2 = new TH1D("psimass2","Psi with selected D",4000,0,8);
258 TH1D* psimass3 = new TH1D("psimass3","Psi with Vertex Fit D",4000,0,8);
259 TH1D* psiptsel = new TH1D("psiptsel","Psi with small D Pt",4000,0,8);
260 TH1D* psiptselfine = new TH1D("psiptselfine","Psi with small Pt",4000,0,8);
261 
262 TH1D* ppmass = new TH1D("ppmass","4C Fit pbarp input mass",8000,0,16);
263 TH1D* ppmass2 = new TH1D("ppmass2","4C Fit pbarp fit mass",8000,0,16);
264 
265 // particle lists
266 TCandList pp;
267 TCandList pipbase, pimbase, kpbase, kmbase;
268 TCandList pip, pim, kp, km;
269 TCandList dpraw, dmraw;
270 TCandList dpmsel, dmmsel;
271 TCandList dpvtx, dmvtx;
272 TCandList rawpsiraw, rawpsimsel, rawpsivtx;
273 TCandList mselpsiraw, mselpsimsel, mselpsivtx;
274 TCandList dptselpsiraw, dptselpsimsel, dptselpsivtx;
275 TCandList pptselpsiraw, pptselpsimsel, pptselpsivtx;
276 TCandList ptselpsiraw, ptselpsimsel, ptselpsivtx;
277 TCandList vtxpsi, ptselvtxpsi, finalpsi;
278 TCandList mct;
279 TCandList dp,dm,pp,psi3770,psiraw,psi3770fit;
280 
281 //TLorentzVector ini(0,0,6.23164,7.24015);
282 TPidMassSelector *jpsiMSel = new TPidMassSelector("jpsiMSel" , 3.096 , 0.3);
283 TPidMassSelector *dMSel = new TPidMassSelector("dMSel", TRho::Instance()->GetPDG()->GetParticle(411)->Mass(), 0.9);
284 TPidMassSelector *dMSelfine = new TPidMassSelector("dMSelfine", TRho::Instance()->GetPDG()->GetParticle(411)->Mass(), 0.5);
285 TPidMassSelector* psiMSel = new TPidMassSelector("psiMSel", TRho::Instance()->GetPDG()->GetParticle(40443)->Mass(), 0.9);
286 
287 Double_t dvtxfitchi2limit = 18;
288 Double_t ptlimit = 0.5;
289 Double_t inipz = 6.5788;
290 Double_t pzlimit = 0.5;
291 
292 int i=0,j=0;
293 
294 TFile outputstorage(outfilename, "UPDATE");
295 TFile mcfile(simfname, "READ");
296 TTree* mctree = (TTree*)mcfile.Get("pndsim");
297 TClonesArray* mcarray = new TClonesArray("PndMCTrack");
298 mctree->SetBranchAddress("MCTrack", &mcarray);
299 
300 TVector3 dpstartvertex;
301 TVector3 dmstartvertex;
302 
303 //PndAnalysis* theAnalysis = new PndAnalysis("BarrelGenTrack");
304 
305 // **** loop over all _events_
306 //
307 while (evr.GetEvent() && ++i<nevts) {
308 
309  // sorry, main loop too long. will not indent code
310 
311 if (!(i%100)) cout <<"evt "<<i<<endl;
312 
313 mctree->GetEntry(i);
314 
315 pp.Cleanup();
316 pip.Cleanup();
317 pim.Cleanup();
318 kp.Cleanup();
319 km.Cleanup();
320 dpraw.Cleanup();
321 dmraw.Cleanup();
322 dpmsel.Cleanup();
323 dmmsel.Cleanup();
324 dpvtx.Cleanup();
325 dmvtx.Cleanup();
326 rawpsiraw.Cleanup();
327 rawpsimsel.Cleanup();
328 rawpsivtx.Cleanup();
329 mselpsiraw.Cleanup();
330 mselpsimsel.Cleanup();
331 mselpsivtx.Cleanup();
332 dptselpsiraw.Cleanup();
333 dptselpsimsel.Cleanup();
334 dptselpsivtx.Cleanup();
335 pptselpsiraw.Cleanup();
336 pptselpsimsel.Cleanup();
337 pptselpsivtx.Cleanup();
338 ptselpsiraw.Cleanup();
339 ptselpsimsel.Cleanup();
340 ptselpsivtx.Cleanup();
341 ptselvtxpsi.Cleanup();
342 vtxpsi.Cleanup();
343 finalpsi.Cleanup();
344 
345 evr.FillList(pipbase,"PionVeryLoosePlus");
346 evr.FillList(pimbase,"PionVeryLooseMinus");
347 evr.FillList(kpbase, "KaonVeryLoosePlus");
348 evr.FillList(kmbase, "KaonVeryLooseMinus");
349 
350 // use Monte Carlo PID information to clean particle lists
351 for (j=0; j<kpbase.GetLength(); ++j) {
352  int mcindex = kpbase[j].GetMcIdx();
353  if (mcindex >= 0) PndMCTrack *mctrack = (PndMCTrack*)mcarray->At(mcindex);
354  if ((mctrack != 0) && (mctrack->GetPdgCode() == 321)) kp.Add(kpbase[j]);
355 }
356 for (j=0; j<kmbase.GetLength(); ++j) {
357  int mcindex = kmbase[j].GetMcIdx();
358  if (mcindex >= 0) PndMCTrack *mctrack = (PndMCTrack*)mcarray->At(mcindex);
359  if ((mctrack != 0) && (mctrack->GetPdgCode() == -321)) km.Add(kmbase[j]);
360 }
361 
362 // store Monte Carlo Start Vertex information for later comparison
363 for (j=0; j<pipbase.GetLength(); ++j) {
364  int mcindex = pipbase[j].GetMcIdx();
365  if (mcindex >= 0) PndMCTrack *mctrack = (PndMCTrack*)mcarray->At(mcindex);
366  if ((mctrack != 0) && (mctrack->GetPdgCode() == 211)) {
367  Double_t x;
368  Double_t y;
369  Double_t z;
370  std::stringstream strstrx;
371  std::stringstream strstry;
372  std::stringstream strstrz;
373  strstrx << mctrack->GetStartVertex().x();
374  strstrx >> x;
375  strstry << mctrack->GetStartVertex().y();
376  strstry >> y;
377  strstrz << mctrack->GetStartVertex().z();
378  strstrz >> z;
379  dpstartvertex.SetXYZ(x, y, z);
380  pip.Add(pipbase[j]);
381  //dpstartvertex = mctrack->GetStartVertex();
382  }
383 }
384 for (j=0; j<pimbase.GetLength(); ++j) {
385  int mcindex = pimbase[j].GetMcIdx();
386  if (mcindex >= 0) PndMCTrack *mctrack = (PndMCTrack*)mcarray->At(mcindex);
387  if ((mctrack != 0) && (mctrack->GetPdgCode() == -211)) {
388  Double_t x;
389  Double_t y;
390  Double_t z;
391  std::stringstream strstrx;
392  std::stringstream strstry;
393  std::stringstream strstrz;
394  strstrx << mctrack->GetStartVertex().x();
395  strstrx >> x;
396  strstry << mctrack->GetStartVertex().y();
397  strstry >> y;
398  strstrz << mctrack->GetStartVertex().z();
399  strstrz >> z;
400  dmstartvertex.SetXYZ(x, y, z);
401  pim.Add(pimbase[j]);
402  //dmstartvertex = mctrack->GetStartVertex();
403  }
404 }
405 
406 hkpperevent->Fill(kp.GetLength());
407 hkmperevent->Fill(km.GetLength());
408 hpipperevent->Fill(pip.GetLength());
409 hpimperevent->Fill(pim.GetLength());
410 if ((kp.GetLength() == 1) && (km.GetLength() == 1) && (pip.GetLength() == 2) && (pim.GetLength() == 2)) hniceevents->Fill(0);
411 if ((kp.GetLength() >= 1) && (km.GetLength() >= 1) && (pip.GetLength() >= 2) && (pim.GetLength() >= 2)) hniceevents->Fill(1);
412 if ((kp.GetLength() + pim.GetLength() == 3) && (km.GetLength() + pip.GetLength() == 3)) hniceevents->Fill(2);
413 if ((kp.GetLength() + pim.GetLength() == 3) || (km.GetLength() + pip.GetLength() == 3)) hniceevents->Fill(3);
414 
415 if ((kp.GetLength() > 0) && (km.GetLength() > 0) && (pip.GetLength() > 1) && (pim.GetLength() > 1)) {
416  dpraw.Combine(km,pip,pip);
417  //cout << "D+ TCandList entries: " << dpraw.GetLength() << endl;
418  dmraw.Combine(kp,pim,pim);
419  //cout << "D- TCandList entries: " << dmraw.GetLength() << endl;
420  dpmsel.Select(dpraw, dMSel);
421  dmmsel.Select(dmraw, dMSel);
422 }
423 
424 for (j=0;j<dpraw.GetLength();++j) {
425  hdprawmass->Fill(dpraw[j].M());
426 }
427 for (j=0;j<dmraw.GetLength();++j) {
428  hdmrawmass->Fill(dmraw[j].M());
429 }
430 for (j=0;j<dpmsel.GetLength();++j) {
431  hdpmselmass->Fill(dpmsel[j].M());
432 }
433 for (j=0;j<dmmsel.GetLength();++j) {
434  hdmmselmass->Fill(dmmsel[j].M());
435 }
436 
437 //POCA test
438 
439 for (j=0; j < dpmsel.GetLength(); ++j) {
440  RhoVtxPoca vtxfinder(dpmsel[j]);
441  TVector3 pocavertex(0,0,0);
442  double distval=-4444;
443  distval = vtxfinder.GetPocaVtx(pocavertex);
444  //reco vertex position
445  hdpvertexxpoca->Fill(pocavertex.x());
446  hdpvertexypoca->Fill(pocavertex.y());
447  hdpvertexzpoca->Fill(pocavertex.z());
448  //poca vertex resolution
449  TVector3 vertexdisp = dpstartvertex - pocavertex;
450  hdpvertexxpocareso->Fill(vertexdisp.x());
451  hdpvertexypocareso->Fill(vertexdisp.y());
452  hdpvertexzpocareso->Fill(vertexdisp.z());
453 }
454 for (j=0; j < dmmsel.GetLength(); ++j) {
455  RhoVtxPoca vtxfinder(dmmsel[j]);
456  TVector3 pocavertex(0,0,0);
457  double distval=-4444;
458  distval = vtxfinder.GetPocaVtx(pocavertex);
459  //reco vertex position
460  hdmvertexxpoca->Fill(pocavertex.x());
461  hdmvertexypoca->Fill(pocavertex.y());
462  hdmvertexzpoca->Fill(pocavertex.z());
463  //poca vertex resolution
464  TVector3 vertexdisp = dmstartvertex - pocavertex;
465  hdmvertexxpocareso->Fill(vertexdisp.x());
466  hdmvertexypocareso->Fill(vertexdisp.y());
467  hdmvertexzpocareso->Fill(vertexdisp.z());
468 }
469 
470 
471 //vertex fit for D+ and D-
472 
473 dpvtx.Cleanup();
474 dmvtx.Cleanup();
475 
476 int bestfitindex;
477 Double_t bestfitchi2 = dvtxfitchi2limit;
478 Double_t bestfitmass = 0;
479 for (j=0; j < dpmsel.GetLength(); ++j) {
480  RhoKinVtxFitter vtxfitter(dpmsel[j]);
481  //vtxfitter.AddMassConstraint(TRho::Instance()->GetPDG()->GetParticle(411)->Mass());
482  vtxfitter.Fit();
483  TCandidate fitcand=*(vtxfitter.FittedCand(dpmsel[j]));
484  TVector3 fitvertex=fitcand.Pos();
485  hdpvtxchi2->Fill(vtxfitter.GlobalChi2());
486  if (vtxfitter.GlobalChi2()<bestfitchi2) {
487  if (bestfitchi2 < dvtxfitchi2limit) {
488  hdpvtxrejectmass->Fill(bestfitmass);
489  }
490  bestfitchi2 = vtxfitter.GlobalChi2();
491  bestfitindex = j;
492  bestfitmass = fitcand.M();
493  } else {
494  hdpvtxrejectmass->Fill(fitcand.M());
495  }
496 }
497 //process best candidate
498 if (bestfitchi2 < dvtxfitchi2limit) {
499  RhoKinVtxFitter vtxfitter(dpmsel[bestfitindex]);
500  //vtxfitter.AddMassConstraint(TRho::Instance()->GetPDG()->GetParticle(411)->Mass());
501  vtxfitter.Fit();
502  TCandidate fitcand=*(vtxfitter.FittedCand(dpmsel[bestfitindex]));
503  TVector3 fitvertex=fitcand.Pos();
504  dpvtx.Add(fitcand);
505  //distance between fitted vertex and MC vertex
506  TVector3 vertexdisp = dpstartvertex - fitvertex;
507  hdpvertexfit->Fill(vertexdisp.Mag());
508  hdpvertexxfit->Fill(vertexdisp.x());
509  hdpvertexyfit->Fill(vertexdisp.y());
510  hdpvertexzfit->Fill(vertexdisp.z());
511  //distance between initial vertex and MC vertex (initial vertex always 0,0,0)
512  vertexdisp = dpstartvertex - dpmsel[bestfitindex].Pos();
513  hdpvertex->Fill(vertexdisp.Mag());
514  hdpvertexx->Fill(vertexdisp.x());
515  hdpvertexy->Fill(vertexdisp.y());
516  hdpvertexz->Fill(vertexdisp.z());
517  //MC vertex position
518  hdpvertexxmc->Fill(dpstartvertex.x());
519  hdpvertexymc->Fill(dpstartvertex.y());
520  hdpvertexzmc->Fill(dpstartvertex.z());
521  //reco vertex position
522  hdpvertexxreco->Fill(fitvertex.x());
523  hdpvertexyreco->Fill(fitvertex.y());
524  hdpvertexzreco->Fill(fitvertex.z());
525  //D+ decay length
526  hdpdecaylength->Fill(fitcand.M()*fitvertex.Mag()/(fitcand.P()*30));
527 }
528 
529 bestfitchi2 = dvtxfitchi2limit;
530 for (j=0; j<dmmsel.GetLength(); ++j) {
531  RhoKinVtxFitter vtxfitter(dmmsel[j]);
532  //vtxfitter.AddMassConstraint(TRho::Instance()->GetPDG()->GetParticle(411)->Mass());
533  vtxfitter.Fit();
534  TCandidate fitcand=*(vtxfitter.FittedCand(dmmsel[j]));
535  TVector3 fitvertex=fitcand.Pos();
536  hdmvtxchi2->Fill(vtxfitter.GlobalChi2());
537  if (vtxfitter.GlobalChi2()<bestfitchi2) {
538  if (bestfitchi2 < dvtxfitchi2limit) {
539  hdmvtxrejectmass->Fill(bestfitmass);
540  }
541  bestfitchi2 = vtxfitter.GlobalChi2();
542  bestfitindex = j;
543  bestfitmass = fitcand.M();
544  } else {
545  hdmvtxrejectmass->Fill(fitcand.M());
546  }
547 }
548 //process best candidate
549 if (bestfitchi2 < dvtxfitchi2limit) {
550  RhoKinVtxFitter vtxfitter(dmmsel[bestfitindex]);
551  //vtxfitter.AddMassConstraint(TRho::Instance()->GetPDG()->GetParticle(411)->Mass());
552  vtxfitter.Fit();
553  TCandidate fitcand=*(vtxfitter.FittedCand(dmmsel[bestfitindex]));
554  TVector3 fitvertex=fitcand.Pos();
555  dmvtx.Add(fitcand);
556  //distance between fitted vertex and MC vertex
557  TVector3 vertexdisp = dmstartvertex - fitvertex;
558  hdmvertexfit->Fill(vertexdisp.Mag());
559  hdmvertexxfit->Fill(vertexdisp.x());
560  hdmvertexyfit->Fill(vertexdisp.y());
561  hdmvertexzfit->Fill(vertexdisp.z());
562  vertexdisp = dmstartvertex - dmmsel[bestfitindex].Pos();
563  hdmvertex->Fill(vertexdisp.Mag());
564  hdmvertexx->Fill(vertexdisp.x());
565  hdmvertexy->Fill(vertexdisp.y());
566  hdmvertexz->Fill(vertexdisp.z());
567  hdmvertexxmc->Fill(dmstartvertex.x());
568  hdmvertexymc->Fill(dmstartvertex.y());
569  hdmvertexzmc->Fill(dmstartvertex.z());
570  hdmvertexxreco->Fill(fitvertex.x());
571  hdmvertexyreco->Fill(fitvertex.y());
572  hdmvertexzreco->Fill(fitvertex.z());
573  hdmdecaylength->Fill(fitcand.M()*fitvertex.Mag()/(fitcand.P()*30));
574 }
575 
576 dpvtx.Select(dMSelfine);
577 dmvtx.Select(dMSelfine);
578 
579 for (j=0;j<dpvtx.GetLength();++j) {
580  hdpvtxacceptmass->Fill(dpvtx[j].M());
581 }
582 for (j=0;j<dmvtx.GetLength();++j) {
583  hdmvtxacceptmass->Fill(dmvtx[j].M());
584 }
585 
586 // Make RAW Psi and prepare lists to hold pt selected Psi
587 
588 TLorentzVector sum;
589 
590 dptselpsiraw.Cleanup();
591 pptselpsiraw.Cleanup();
592 ptselpsiraw.Cleanup();
593 dptselpsimsel.Cleanup();
594 pptselpsimsel.Cleanup();
595 ptselpsimsel.Cleanup();
596 dptselpsivtx.Cleanup();
597 pptselpsivtx.Cleanup();
598 ptselpsivtx.Cleanup();
599 
600 // Psi from all D mesons
601 rawpsiraw.Combine(dpraw,dmraw);
602 for (j=0;j<rawpsiraw.GetLength();++j) {
603  hrawpsirawmass->Fill(rawpsiraw[j].M());
604  hrawpsirawpt->Fill(rawpsiraw[j].Pt());
605  hrawpsirawpz->Fill(rawpsiraw[j].Pz());
606  sum = rawpsiraw[j].Daughter(0)->P4() + rawpsiraw[j].Daughter(1)->P4();
607  hrawpsirawdpt->Fill(sum.Pt());
608  hrawpsirawdpz->Fill(sum.Pz());
609  if (sum.Pt() < ptlimit) {
610  dptselpsiraw.Add(rawpsiraw[j]);
611  }
612  if (rawpsiraw[j].Pt() < ptlimit) {
613  pptselpsiraw.Add(rawpsiraw[j]);
614  }
615  if ((sum.Pt() < ptlimit)&&(rawpsiraw[j].Pt() < ptlimit)) {
616  ptselpsiraw.Add(rawpsiraw[j]);
617  }
618 }
619 
620 // Psi from mass selected D mesons
621 rawpsimsel.Combine(dpmsel,dmmsel);
622 for (j=0;j<rawpsimsel.GetLength();++j) {
623  hrawpsimselmass->Fill(rawpsimsel[j].M());
624  hrawpsimselpt->Fill(rawpsimsel[j].Pt());
625  hrawpsimselpz->Fill(rawpsimsel[j].Pz());
626  sum = rawpsimsel[j].Daughter(0)->P4() + rawpsimsel[j].Daughter(1)->P4();
627  hrawpsimseldpt->Fill(sum.Pt());
628  hrawpsimseldpz->Fill(sum.Pz());
629  if (sum.Pt() < ptlimit) {
630  dptselpsimsel.Add(rawpsimsel[j]);
631  }
632  if (rawpsimsel[j].Pt() < ptlimit) {
633  pptselpsimsel.Add(rawpsimsel[j]);
634  }
635  if ((sum.Pt() < ptlimit)&&(rawpsimsel[j].Pt() < ptlimit)) {
636  ptselpsimsel.Add(rawpsimsel[j]);
637  }
638 }
639 
640 // Psi from vertex fitted D mesons
641 rawpsivtx.Combine(dpvtx,dmvtx);
642 for (j=0;j<rawpsivtx.GetLength();++j) {
643  hrawpsivtxmass->Fill(rawpsivtx[j].M());
644  hrawpsivtxpt->Fill(rawpsivtx[j].Pt());
645  hrawpsivtxpz->Fill(rawpsivtx[j].Pz());
646  sum = rawpsivtx[j].Daughter(0)->P4() + rawpsivtx[j].Daughter(1)->P4();
647  hrawpsivtxdpt->Fill(sum.Pt());
648  hrawpsivtxdpz->Fill(sum.Pz());
649  if (sum.Pt() < ptlimit) {
650  dptselpsivtx.Add(rawpsivtx[j]);
651  }
652  if (rawpsivtx[j].Pt() < ptlimit) {
653  pptselpsivtx.Add(rawpsivtx[j]);
654  }
655  if ((sum.Pt() < ptlimit)&&(rawpsivtx[j].Pt() < ptlimit)) {
656  ptselpsivtx.Add(rawpsivtx[j]);
657  }
658  if ((sum.Pt() < ptlimit)&&(rawpsivtx[j].Pt() < ptlimit) && (rawpsivtx[j].Pz() > (inipz - pzlimit)) && (rawpsivtx[j].Pz() < (inipz + pzlimit))) {
659  finalpsi.Add(rawpsivtx[j]);
660  }
661 }
662 
663 for (j=0;j<finalpsi.GetLength();++j) {
664  hfinalpsimass->Fill(finalpsi[j].M());
665 }
666 
667 // make mass selected Psi from Psi from all D mesons
668 mselpsiraw.Select(rawpsiraw, psiMSel);
669 for (j=0;j<mselpsiraw.GetLength();++j) {
670  hmselpsirawmass->Fill(mselpsiraw[j].M());
671  hmselpsirawpt->Fill(mselpsiraw[j].Pt());
672  hmselpsirawpz->Fill(mselpsiraw[j].Pz());
673  sum = mselpsiraw[j].Daughter(0)->P4() + mselpsiraw[j].Daughter(1)->P4();
674  hmselpsirawdpt->Fill(sum.Pt());
675  hmselpsirawdpz->Fill(sum.Pz());
676 }
677 
678 // make mass selected Psi from Psi from mass selected D mesons
679 mselpsimsel.Select(rawpsimsel, psiMSel);
680 for (j=0;j<mselpsimsel.GetLength();++j) {
681  hmselpsimselmass->Fill(mselpsimsel[j].M());
682  hmselpsimselpt->Fill(mselpsimsel[j].Pt());
683  hmselpsimselpz->Fill(mselpsimsel[j].Pz());
684  sum = mselpsimsel[j].Daughter(0)->P4() + mselpsimsel[j].Daughter(1)->P4();
685  hmselpsimseldpt->Fill(sum.Pt());
686  hmselpsimseldpz->Fill(sum.Pz());
687 }
688 
689 // make mass selected Psi from Psi from vertex fitted D mesons
690 mselpsivtx.Select(rawpsivtx, psiMSel);
691 for (j=0;j<mselpsivtx.GetLength();++j) {
692  hmselpsivtxmass->Fill(mselpsivtx[j].M());
693  hmselpsivtxpt->Fill(mselpsivtx[j].Pt());
694  hmselpsivtxpz->Fill(mselpsivtx[j].Pz());
695  sum = mselpsivtx[j].Daughter(0)->P4() + mselpsivtx[j].Daughter(1)->P4();
696  hmselpsivtxdpt->Fill(sum.Pt());
697  hmselpsivtxdpz->Fill(sum.Pz());
698 }
699 
700 
701 // plot pt (D) selected Psi
702 for (j=0;j<dptselpsiraw.GetLength();++j) {
703  hdptselpsirawmass->Fill(dptselpsiraw[j].M());
704  hdptselpsirawpt->Fill(dptselpsiraw[j].Pt());
705  hdptselpsirawpz->Fill(dptselpsiraw[j].Pz());
706  sum = dptselpsiraw[j].Daughter(0)->P4() + dptselpsiraw[j].Daughter(1)->P4();
707  hdptselpsirawdpt->Fill(sum.Pt());
708  hdptselpsirawdpz->Fill(sum.Pz());
709 }
710 for (j=0;j<dptselpsimsel.GetLength();++j) {
711  hdptselpsimselmass->Fill(dptselpsimsel[j].M());
712  hdptselpsimselpt->Fill(dptselpsimsel[j].Pt());
713  hdptselpsimselpz->Fill(dptselpsimsel[j].Pz());
714  sum = dptselpsimsel[j].Daughter(0)->P4() + dptselpsimsel[j].Daughter(1)->P4();
715  hdptselpsimseldpt->Fill(sum.Pt());
716  hdptselpsimseldpz->Fill(sum.Pz());
717 }
718 for (j=0;j<dptselpsivtx.GetLength();++j) {
719  hdptselpsivtxmass->Fill(dptselpsivtx[j].M());
720  hdptselpsivtxpt->Fill(dptselpsivtx[j].Pt());
721  hdptselpsivtxpz->Fill(dptselpsivtx[j].Pz());
722  sum = dptselpsivtx[j].Daughter(0)->P4() + dptselpsivtx[j].Daughter(1)->P4();
723  hdptselpsivtxdpt->Fill(sum.Pt());
724  hdptselpsivtxdpz->Fill(sum.Pz());
725 }
726 
727 // plot pt (Psi, momentum should be the same as D sum) selected Psi
728 for (j=0;j<pptselpsiraw.GetLength();++j) {
729  hpptselpsirawmass->Fill(pptselpsiraw[j].M());
730  hpptselpsirawpt->Fill(pptselpsiraw[j].Pt());
731  hpptselpsirawpz->Fill(pptselpsiraw[j].Pz());
732  sum = pptselpsiraw[j].Daughter(0)->P4() + pptselpsiraw[j].Daughter(1)->P4();
733  hpptselpsirawdpt->Fill(sum.Pt());
734  hpptselpsirawdpz->Fill(sum.Pz());
735 }
736 for (j=0;j<pptselpsimsel.GetLength();++j) {
737  hpptselpsimselmass->Fill(pptselpsimsel[j].M());
738  hpptselpsimselpt->Fill(pptselpsimsel[j].Pt());
739  hpptselpsimselpz->Fill(pptselpsimsel[j].Pz());
740  sum = pptselpsimsel[j].Daughter(0)->P4() + pptselpsimsel[j].Daughter(1)->P4();
741  hpptselpsimseldpt->Fill(sum.Pt());
742  hpptselpsimseldpz->Fill(sum.Pz());
743 }
744 for (j=0;j<pptselpsivtx.GetLength();++j) {
745  hpptselpsivtxmass->Fill(pptselpsivtx[j].M());
746  hpptselpsivtxpt->Fill(pptselpsivtx[j].Pt());
747  hpptselpsivtxpz->Fill(pptselpsivtx[j].Pz());
748  sum = pptselpsivtx[j].Daughter(0)->P4() + pptselpsivtx[j].Daughter(1)->P4();
749  hpptselpsivtxdpt->Fill(sum.Pt());
750  hpptselpsivtxdpz->Fill(sum.Pz());
751 }
752 
753 // plot pt (both D and Psi, should be redundant) selected Psi
754 for (j=0;j<ptselpsiraw.GetLength();++j) {
755  hptselpsirawmass->Fill(ptselpsiraw[j].M());
756  hptselpsirawpt->Fill(ptselpsiraw[j].Pt());
757  hptselpsirawpz->Fill(ptselpsiraw[j].Pz());
758  sum = ptselpsiraw[j].Daughter(0)->P4() + ptselpsiraw[j].Daughter(1)->P4();
759  hptselpsirawdpt->Fill(sum.Pt());
760  hptselpsirawdpz->Fill(sum.Pz());
761 }
762 for (j=0;j<ptselpsimsel.GetLength();++j) {
763  hptselpsimselmass->Fill(ptselpsimsel[j].M());
764  hptselpsimselpt->Fill(ptselpsimsel[j].Pt());
765  hptselpsimselpz->Fill(ptselpsimsel[j].Pz());
766  sum = ptselpsimsel[j].Daughter(0)->P4() + ptselpsimsel[j].Daughter(1)->P4();
767  hptselpsimseldpt->Fill(sum.Pt());
768  hptselpsimseldpz->Fill(sum.Pz());
769 }
770 for (j=0;j<ptselpsivtx.GetLength();++j) {
771  hptselpsivtxmass->Fill(ptselpsivtx[j].M());
772  hptselpsivtxpt->Fill(ptselpsivtx[j].Pt());
773  hptselpsivtxpz->Fill(ptselpsivtx[j].Pz());
774  sum = ptselpsivtx[j].Daughter(0)->P4() + ptselpsivtx[j].Daughter(1)->P4();
775  hptselpsivtxdpt->Fill(sum.Pt());
776  hptselpsivtxdpz->Fill(sum.Pz());
777 }
778 
779 // Make Psi Vertex Fit
780 
781 for (j=0;j<ptselpsivtx.GetLength();++j) {
782  RhoKinVtxFitter fitter(ptselpsivtx[j]);
783  fitter.Fit();
784  TCandidate fitcand=*(fitter.FittedCand(ptselpsivtx[j]));
785  TCandidate dpfit=*(fitter.FittedCand(*(ptselpsivtx[j].Daughter(0))));
786  TCandidate dmfit=*(fitter.FittedCand(*(ptselpsivtx[j].Daughter(1))));
787  //TCandidate* dpfit = fitcand->Daughter(0);
788  //TCandidate* dmfit = fitcand->Daughter(1);
789  TLorentzVector sum=dpfit.P4()+dmfit.P4();
790  hvtxpsichi2->Fill(fitter.GlobalChi2());
791  if (fitter.GlobalChi2()<50) {
792  hvtxpsiacceptmass->Fill(fitcand.M());
793  hvtxpsiacceptpt->Fill(fitcand.Pt());
794  hvtxpsiacceptpz->Fill(fitcand.Pz());
795  hvtxpsiacceptdpt->Fill(sum.Pt());
796  hvtxpsiacceptdpz->Fill(sum.Pz());
797  vtxpsi.Add(fitcand);
798  //pt selection here for now
799  if (fitcand.Pt() < 0.5) {
800  hptselvtxpsimass->Fill(fitcand.M());
801  hptselvtxpsipt->Fill(fitcand.Pt());
802  hptselvtxpsipz->Fill(fitcand.Pz());
803  hptselvtxpsidpt->Fill(sum.Pt());
804  hptselvtxpsidpz->Fill(sum.Pz());
805  }
806  } else {
807  hvtxpsirejectmass->Fill(fitcand.M());
808  hvtxpsirejectpt->Fill(fitcand.Pt());
809  hvtxpsirejectpz->Fill(fitcand.Pz());
810  hvtxpsirejectdpt->Fill(sum.Pt());
811  hvtxpsirejectdpz->Fill(sum.Pz());
812  }
813 }
814 
815 } // end main loop
816 
817 //TLorentzVector ini(0,0,6.5788,7.58364);
818 
819 //plot histograms
820 
821 TCanvas *eventstatisticscanvas = new TCanvas("eventstatisticscanvas","eventstatisticscanvas",600,600);
822 eventstatisticscanvas->Divide(2,3);
823  eventstatisticscanvas->cd(1); hniceevents->Draw();
824  eventstatisticscanvas->cd(3); hkpperevent->Draw();
825  eventstatisticscanvas->cd(4); hkmperevent->Draw();
826  eventstatisticscanvas->cd(5); hpipperevent->Draw();
827  eventstatisticscanvas->cd(6); hpimperevent->Draw();
828 
829 TCanvas *dcanvas=new TCanvas("dcanvas","dcanvas",600,600);
830 dcanvas->Divide(2,4);
831  dcanvas->cd(1); hdprawmass->Draw();
832  dcanvas->cd(2); hdmrawmass->Draw();
833  dcanvas->cd(3); hdpmselmass->Draw();
834  dcanvas->cd(4); hdmmselmass->Draw();
835  dcanvas->cd(5); hdpvtxacceptmass->Draw(); hdpvtxrejectmass->Draw("same");
836  dcanvas->cd(6); hdmvtxacceptmass->Draw(); hdmvtxrejectmass->Draw("same");
837  dcanvas->cd(7); hdpvtxchi2->Draw(); //psimass3->Draw(); //ppmass->Draw();
838  dcanvas->cd(8); hdmvtxchi2->Draw(); //ppmass->Draw();
839  dcanvas->cd();
840 
841 TCanvas *rawpsirawcanvas = new TCanvas("rawpsirawcanvas","rawpsirawcanvas",600,600);
842 rawpsirawcanvas->Divide(3,2);
843  rawpsirawcanvas->cd(1); hrawpsirawmass->Draw(); hmselpsirawmass->Draw("same");
844  rawpsirawcanvas->cd(2); hrawpsirawpt->Draw(); hmselpsirawpt->Draw("same");
845  rawpsirawcanvas->cd(3); hrawpsirawpz->Draw(); hmselpsirawpz->Draw("same");
846  rawpsirawcanvas->cd(5); hrawpsirawdpt->Draw(); hmselpsirawdpt->Draw("same");
847  rawpsirawcanvas->cd(6); hrawpsirawdpz->Draw(); hmselpsirawdpz->Draw("same");
848 
849 TCanvas *rawpsimselcanvas = new TCanvas("rawpsimselcanvas","rawpsimselcanvas",600,600);
850 rawpsimselcanvas->Divide(3,2);
851  rawpsimselcanvas->cd(1); hrawpsimselmass->Draw(); hmselpsimselmass->Draw("same");
852  rawpsimselcanvas->cd(2); hrawpsimselpt->Draw(); hmselpsimselpt->Draw("same");
853  rawpsimselcanvas->cd(3); hrawpsimselpz->Draw(); hmselpsimselpz->Draw("same");
854  rawpsimselcanvas->cd(5); hrawpsimseldpt->Draw(); hmselpsimseldpt->Draw("same");
855  rawpsimselcanvas->cd(6); hrawpsimseldpz->Draw(); hmselpsimseldpz->Draw("same");
856 
857 TCanvas *rawpsivtxcanvas = new TCanvas("rawpsivtxcanvas","rawpsivtxcanvas",600,600);
858 rawpsivtxcanvas->Divide(3,2);
859  rawpsivtxcanvas->cd(1); hrawpsivtxmass->Draw(); hmselpsivtxmass->Draw("same");
860  rawpsivtxcanvas->cd(2); hrawpsivtxpt->Draw(); hmselpsivtxpt->Draw("same");
861  rawpsivtxcanvas->cd(3); hrawpsivtxpz->Draw(); hmselpsivtxpz->Draw("same");
862  rawpsivtxcanvas->cd(4); hfinalpsimass->Draw();
863  rawpsivtxcanvas->cd(5); hrawpsivtxdpt->Draw(); hmselpsivtxdpt->Draw("same");
864  rawpsivtxcanvas->cd(6); hrawpsivtxdpz->Draw(); hmselpsivtxdpz->Draw("same");
865 
866 TCanvas *dptselpsirawcanvas = new TCanvas("dptselpsirawcanvas","dptselpsirawcanvas",600,600);
867 dptselpsirawcanvas->Divide(3,2);
868  dptselpsirawcanvas->cd(1); hdptselpsirawmass->Draw();
869  dptselpsirawcanvas->cd(2); hdptselpsirawpt->Draw();
870  dptselpsirawcanvas->cd(3); hdptselpsirawpz->Draw();
871  dptselpsirawcanvas->cd(5); hdptselpsirawdpt->Draw();
872  dptselpsirawcanvas->cd(6); hdptselpsirawdpz->Draw();
873 
874 TCanvas *dptselpsimselcanvas = new TCanvas("dptselpsimselcanvas","dptselpsimselcanvas",600,600);
875 dptselpsimselcanvas->Divide(3,2);
876  dptselpsimselcanvas->cd(1); hdptselpsimselmass->Draw();
877  dptselpsimselcanvas->cd(2); hdptselpsimselpt->Draw();
878  dptselpsimselcanvas->cd(3); hdptselpsimselpz->Draw();
879  dptselpsimselcanvas->cd(5); hdptselpsimseldpt->Draw();
880  dptselpsimselcanvas->cd(6); hdptselpsimseldpz->Draw();
881 
882 TCanvas *dptselpsivtxcanvas = new TCanvas("dptselpsivtxcanvas","dptselpsivtxcanvas",600,600);
883 dptselpsivtxcanvas->Divide(3,2);
884  dptselpsivtxcanvas->cd(1); hdptselpsivtxmass->Draw();
885  dptselpsivtxcanvas->cd(2); hdptselpsivtxpt->Draw();
886  dptselpsivtxcanvas->cd(3); hdptselpsivtxpz->Draw();
887  dptselpsivtxcanvas->cd(5); hdptselpsivtxdpt->Draw();
888  dptselpsivtxcanvas->cd(6); hdptselpsivtxdpz->Draw();
889 
890 TCanvas *pptselpsirawcanvas = new TCanvas("pptselpsirawcanvas","pptselpsirawcanvas",600,600);
891 pptselpsirawcanvas->Divide(3,2);
892  pptselpsirawcanvas->cd(1); hpptselpsirawmass->Draw();
893  pptselpsirawcanvas->cd(2); hpptselpsirawpt->Draw();
894  pptselpsirawcanvas->cd(3); hpptselpsirawpz->Draw();
895  pptselpsirawcanvas->cd(5); hpptselpsirawdpt->Draw();
896  pptselpsirawcanvas->cd(6); hpptselpsirawdpz->Draw();
897 
898 TCanvas *pptselpsimselcanvas = new TCanvas("pptselpsimselcanvas","pptselpsimselcanvas",600,600);
899 pptselpsimselcanvas->Divide(3,2);
900  pptselpsimselcanvas->cd(1); hpptselpsimselmass->Draw();
901  pptselpsimselcanvas->cd(2); hpptselpsimselpt->Draw();
902  pptselpsimselcanvas->cd(3); hpptselpsimselpz->Draw();
903  pptselpsimselcanvas->cd(5); hpptselpsimseldpt->Draw();
904  pptselpsimselcanvas->cd(6); hpptselpsimseldpz->Draw();
905 
906 TCanvas *pptselpsivtxcanvas = new TCanvas("pptselpsivtxcanvas","pptselpsivtxcanvas",600,600);
907 pptselpsivtxcanvas->Divide(3,2);
908  pptselpsivtxcanvas->cd(1); hpptselpsivtxmass->Draw();
909  pptselpsivtxcanvas->cd(2); hpptselpsivtxpt->Draw();
910  pptselpsivtxcanvas->cd(3); hpptselpsivtxpz->Draw();
911  pptselpsivtxcanvas->cd(5); hpptselpsivtxdpt->Draw();
912  pptselpsivtxcanvas->cd(6); hpptselpsivtxdpz->Draw();
913 
914 TCanvas *ptselpsirawcanvas = new TCanvas("ptselpsirawcanvas","ptselpsirawcanvas",600,600);
915 ptselpsirawcanvas->Divide(3,2);
916  ptselpsirawcanvas->cd(1); hptselpsirawmass->Draw();
917  ptselpsirawcanvas->cd(2); hptselpsirawpt->Draw();
918  ptselpsirawcanvas->cd(3); hptselpsirawpz->Draw();
919  ptselpsirawcanvas->cd(5); hptselpsirawdpt->Draw();
920  ptselpsirawcanvas->cd(6); hptselpsirawdpz->Draw();
921 
922 TCanvas *ptselpsimselcanvas = new TCanvas("ptselpsimselcanvas","ptselpsimselcanvas",600,600);
923 ptselpsimselcanvas->Divide(3,2);
924  ptselpsimselcanvas->cd(1); hptselpsimselmass->Draw();
925  ptselpsimselcanvas->cd(2); hptselpsimselpt->Draw();
926  ptselpsimselcanvas->cd(3); hptselpsimselpz->Draw();
927  ptselpsimselcanvas->cd(5); hptselpsimseldpt->Draw();
928  ptselpsimselcanvas->cd(6); hptselpsimseldpz->Draw();
929 
930 TCanvas *ptselpsivtxcanvas = new TCanvas("ptselpsivtxcanvas","ptselpsivtxcanvas",600,600);
931 ptselpsivtxcanvas->Divide(3,2);
932  ptselpsivtxcanvas->cd(1); hptselpsivtxmass->Draw();
933  ptselpsivtxcanvas->cd(2); hptselpsivtxpt->Draw();
934  ptselpsivtxcanvas->cd(3); hptselpsivtxpz->Draw();
935  ptselpsivtxcanvas->cd(5); hptselpsivtxdpt->Draw();
936  ptselpsivtxcanvas->cd(6); hptselpsivtxdpz->Draw();
937 
938 TCanvas* vtxpsicanvas = new TCanvas("vtxpsicanvas","vtxpsicanvas",600,600);
939 vtxpsicanvas->Divide(3,2);
940  vtxpsicanvas->cd(1); hvtxpsichi2->Draw();
941  vtxpsicanvas->cd(2); hvtxpsiacceptmass->Draw(); hvtxpsirejectmass->Draw("same");
942  vtxpsicanvas->cd(3); hvtxpsiacceptpt->Draw(); hvtxpsirejectpt->Draw("same");
943  vtxpsicanvas->cd(4); hvtxpsiacceptpz->Draw(); hvtxpsirejectpz->Draw("same");
944  vtxpsicanvas->cd(5); hvtxpsiacceptdpt->Draw(); hvtxpsirejectdpt->Draw("same");
945  vtxpsicanvas->cd(6); hvtxpsiacceptdpz->Draw(); hvtxpsirejectdpz->Draw("same");
946 
947 TCanvas* ptselvtxpsicanvas = new TCanvas("ptselvtxpsicanvas","ptselvtxpsicanvas",600,600);
948 ptselvtxpsicanvas->Divide(2,3);
949  ptselvtxpsicanvas->cd(1); hptselvtxpsimass->Draw();
950  ptselvtxpsicanvas->cd(3); hptselvtxpsipt->Draw();
951  ptselvtxpsicanvas->cd(4); hptselvtxpsipz->Draw();
952  ptselvtxpsicanvas->cd(5); hptselvtxpsidpt->Draw();
953  ptselvtxpsicanvas->cd(6); hptselvtxpsidpz->Draw();;
954 
955 TCanvas* vertexresolutioncanvas = new TCanvas("vertexresolutioncanvas","vertexresolutioncanvas",600,600);
956 vertexresolutioncanvas->Divide(2,2);
957  vertexresolutioncanvas->cd(1); hdpvertex->Draw();
958  vertexresolutioncanvas->cd(2); hdpvertexfit->Draw();
959  vertexresolutioncanvas->cd(3); hdmvertex->Draw();
960  vertexresolutioncanvas->cd(4); hdmvertexfit->Draw();
961 
962 TCanvas* vertexresolutioncanvas2 = new TCanvas("vertexresolutioncanvas2","vertexresolutioncanvas2",600,600);
963 vertexresolutioncanvas2->Divide(3,3);
964  vertexresolutioncanvas2->cd(1); hdmvertexx->Draw();
965  vertexresolutioncanvas2->cd(2); hdmvertexy->Draw();
966  vertexresolutioncanvas2->cd(3); hdmvertexz->Draw();
967  vertexresolutioncanvas2->cd(4); hdmvertexxfit->Draw();
968  vertexresolutioncanvas2->cd(5); hdmvertexyfit->Draw();
969  vertexresolutioncanvas2->cd(6); hdmvertexzfit->Draw();
970  vertexresolutioncanvas2->cd(7); hdmvertexxpocareso->Draw();
971  vertexresolutioncanvas2->cd(8); hdmvertexypocareso->Draw();
972  vertexresolutioncanvas2->cd(9); hdmvertexzpocareso->Draw();
973 
974 TCanvas* vertexresolutioncanvas3 = new TCanvas("vertexresolutioncanvas3","vertexresolutioncanvas3",600,600);
975 vertexresolutioncanvas3->Divide(3,3);
976  vertexresolutioncanvas3->cd(1); hdpvertexx->Draw();
977  vertexresolutioncanvas3->cd(2); hdpvertexy->Draw();
978  vertexresolutioncanvas3->cd(3); hdpvertexz->Draw();
979  vertexresolutioncanvas3->cd(4); hdpvertexxfit->Draw();
980  vertexresolutioncanvas3->cd(5); hdpvertexyfit->Draw();
981  vertexresolutioncanvas3->cd(6); hdpvertexzfit->Draw();
982  vertexresolutioncanvas3->cd(7); hdpvertexxpocareso->Draw();
983  vertexresolutioncanvas3->cd(8); hdpvertexypocareso->Draw();
984  vertexresolutioncanvas3->cd(9); hdpvertexzpocareso->Draw();
985 
986 TCanvas* vertexpositionpcanvas = new TCanvas("vertexpositionpcanvas","vertexpositionpcanvas",600,600);
987 vertexpositionpcanvas->Divide(3,3);
988  vertexpositionpcanvas->cd(1); hdpvertexxmc->Draw();
989  vertexpositionpcanvas->cd(2); hdpvertexymc->Draw();
990  vertexpositionpcanvas->cd(3); hdpvertexzmc->Draw();
991  vertexpositionpcanvas->cd(4); hdpvertexxreco->Draw();
992  vertexpositionpcanvas->cd(5); hdpvertexyreco->Draw();
993  vertexpositionpcanvas->cd(6); hdpvertexzreco->Draw();
994  vertexpositionpcanvas->cd(7); hdpvertexxpoca->Draw();
995  vertexpositionpcanvas->cd(8); hdpvertexypoca->Draw();
996  vertexpositionpcanvas->cd(9); hdpvertexzpoca->Draw();
997 
998 TCanvas* vertexpositionmcanvas = new TCanvas("vertexpositionmcanvas","vertexpositionmcanvas",600,600);
999 vertexpositionmcanvas->Divide(3,3);
1000  vertexpositionmcanvas->cd(1); hdmvertexxmc->Draw();
1001  vertexpositionmcanvas->cd(2); hdmvertexymc->Draw();
1002  vertexpositionmcanvas->cd(3); hdmvertexzmc->Draw();
1003  vertexpositionmcanvas->cd(4); hdmvertexxreco->Draw();
1004  vertexpositionmcanvas->cd(5); hdmvertexyreco->Draw();
1005  vertexpositionmcanvas->cd(6); hdmvertexzreco->Draw();
1006  vertexpositionmcanvas->cd(7); hdmvertexxpoca->Draw();
1007  vertexpositionmcanvas->cd(8); hdmvertexypoca->Draw();
1008  vertexpositionmcanvas->cd(9); hdmvertexzpoca->Draw();
1009 
1010 
1011 TCanvas* decaylengthcanvas = new TCanvas("decaylengthcanvas","decaylengthcanvas",600,600);
1012 decaylengthcanvas->Divide(2,1);
1013  decaylengthcanvas->cd(1); hdpdecaylength->Draw();
1014  decaylengthcanvas->cd(2); hdmdecaylength->Draw();
1015 
1016 //save and update the interesting histos
1017 outputstorage.cd();
1018 SaveAndUpdateHisto(hniceevents, outputstorage);
1019 SaveAndUpdateHisto(hkpperevent, outputstorage);
1020 SaveAndUpdateHisto(hkmperevent, outputstorage);
1021 SaveAndUpdateHisto(hpipperevent, outputstorage);
1022 SaveAndUpdateHisto(hpimperevent, outputstorage);
1023 
1024 SaveAndUpdateHisto(hdprawmass, outputstorage);
1025 SaveAndUpdateHisto(hdmrawmass, outputstorage);
1026 SaveAndUpdateHisto(hdpmselmass, outputstorage);
1027 SaveAndUpdateHisto(hdmmselmass, outputstorage);
1028 SaveAndUpdateHisto(hdpvtxacceptmass, outputstorage);
1029 SaveAndUpdateHisto(hdpvtxrejectmass, outputstorage);
1030 SaveAndUpdateHisto(hdmvtxacceptmass, outputstorage);
1031 SaveAndUpdateHisto(hdmvtxrejectmass, outputstorage);
1032 SaveAndUpdateHisto(hdpvtxchi2, outputstorage);
1033 SaveAndUpdateHisto(hdmvtxchi2, outputstorage);
1034 
1035 SaveAndUpdateHisto(hrawpsirawmass, outputstorage);
1036 SaveAndUpdateHisto(hmselpsirawmass, outputstorage);
1037 SaveAndUpdateHisto(hrawpsirawpt, outputstorage);
1038 SaveAndUpdateHisto(hmselpsirawpt, outputstorage);
1039 SaveAndUpdateHisto(hrawpsirawpz, outputstorage);
1040 SaveAndUpdateHisto(hmselpsirawpz, outputstorage);
1041 SaveAndUpdateHisto(hrawpsirawdpt, outputstorage);
1042 SaveAndUpdateHisto(hmselpsirawdpt, outputstorage);
1043 SaveAndUpdateHisto(hrawpsirawdpz, outputstorage);
1044 SaveAndUpdateHisto(hmselpsirawdpz, outputstorage);
1045 SaveAndUpdateHisto(hrawpsimselmass, outputstorage);
1046 SaveAndUpdateHisto(hmselpsimselmass, outputstorage);
1047 SaveAndUpdateHisto(hrawpsimselpt, outputstorage);
1048 SaveAndUpdateHisto(hmselpsimselpt, outputstorage);
1049 SaveAndUpdateHisto(hrawpsimselpz, outputstorage);
1050 SaveAndUpdateHisto(hmselpsimselpz, outputstorage);
1051 SaveAndUpdateHisto(hrawpsimseldpt, outputstorage);
1052 SaveAndUpdateHisto(hmselpsimseldpt, outputstorage);
1053 SaveAndUpdateHisto(hrawpsimseldpz, outputstorage);
1054 SaveAndUpdateHisto(hmselpsimseldpz, outputstorage);
1055 SaveAndUpdateHisto(hrawpsivtxmass, outputstorage);
1056 SaveAndUpdateHisto(hmselpsivtxmass, outputstorage);
1057 SaveAndUpdateHisto(hrawpsivtxpt, outputstorage);
1058 SaveAndUpdateHisto(hmselpsivtxpt, outputstorage);
1059 SaveAndUpdateHisto(hrawpsivtxpz, outputstorage);
1060 SaveAndUpdateHisto(hmselpsivtxpz, outputstorage);
1061 SaveAndUpdateHisto(hfinalpsimass, outputstorage);
1062 SaveAndUpdateHisto(hrawpsivtxdpt, outputstorage);
1063 SaveAndUpdateHisto(hmselpsivtxdpt, outputstorage);
1064 SaveAndUpdateHisto(hrawpsivtxdpz, outputstorage);
1065 SaveAndUpdateHisto(hmselpsivtxdpz, outputstorage);
1066 
1067 SaveAndUpdateHisto(hdptselpsirawmass, outputstorage);
1068 SaveAndUpdateHisto(hdptselpsirawpt, outputstorage);
1069 SaveAndUpdateHisto(hdptselpsirawpz, outputstorage);
1070 SaveAndUpdateHisto(hdptselpsirawdpt, outputstorage);
1071 SaveAndUpdateHisto(hdptselpsirawdpz, outputstorage);
1072 SaveAndUpdateHisto(hdptselpsimselmass, outputstorage);
1073 SaveAndUpdateHisto(hdptselpsimselpt, outputstorage);
1074 SaveAndUpdateHisto(hdptselpsimselpz, outputstorage);
1075 SaveAndUpdateHisto(hdptselpsimseldpt, outputstorage);
1076 SaveAndUpdateHisto(hdptselpsimseldpz, outputstorage);
1077 SaveAndUpdateHisto(hdptselpsivtxmass, outputstorage);
1078 SaveAndUpdateHisto(hdptselpsivtxpt, outputstorage);
1079 SaveAndUpdateHisto(hdptselpsivtxpz, outputstorage);
1080 SaveAndUpdateHisto(hdptselpsivtxdpt, outputstorage);
1081 SaveAndUpdateHisto(hdptselpsivtxdpz, outputstorage);
1082 
1083 SaveAndUpdateHisto(hpptselpsirawmass, outputstorage);
1084 SaveAndUpdateHisto(hpptselpsirawpt, outputstorage);
1085 SaveAndUpdateHisto(hpptselpsirawpz, outputstorage);
1086 SaveAndUpdateHisto(hpptselpsirawdpt, outputstorage);
1087 SaveAndUpdateHisto(hpptselpsirawdpz, outputstorage);
1088 SaveAndUpdateHisto(hpptselpsimselmass, outputstorage);
1089 SaveAndUpdateHisto(hpptselpsimselpt, outputstorage);
1090 SaveAndUpdateHisto(hpptselpsimselpz, outputstorage);
1091 SaveAndUpdateHisto(hpptselpsimseldpt, outputstorage);
1092 SaveAndUpdateHisto(hpptselpsimseldpz, outputstorage);
1093 SaveAndUpdateHisto(hpptselpsivtxmass, outputstorage);
1094 SaveAndUpdateHisto(hpptselpsivtxpt, outputstorage);
1095 SaveAndUpdateHisto(hpptselpsivtxpz, outputstorage);
1096 SaveAndUpdateHisto(hpptselpsivtxdpt, outputstorage);
1097 SaveAndUpdateHisto(hpptselpsivtxdpz, outputstorage);
1098 
1099 SaveAndUpdateHisto(hptselpsirawmass, outputstorage);
1100 SaveAndUpdateHisto(hptselpsirawpt, outputstorage);
1101 SaveAndUpdateHisto(hptselpsirawpz, outputstorage);
1102 SaveAndUpdateHisto(hptselpsirawdpt, outputstorage);
1103 SaveAndUpdateHisto(hptselpsirawdpz, outputstorage);
1104 SaveAndUpdateHisto(hptselpsimselmass, outputstorage);
1105 SaveAndUpdateHisto(hptselpsimselpt, outputstorage);
1106 SaveAndUpdateHisto(hptselpsimselpz, outputstorage);
1107 SaveAndUpdateHisto(hptselpsimseldpt, outputstorage);
1108 SaveAndUpdateHisto(hptselpsimseldpz, outputstorage);
1109 SaveAndUpdateHisto(hptselpsivtxmass, outputstorage);
1110 SaveAndUpdateHisto(hptselpsivtxpt, outputstorage);
1111 SaveAndUpdateHisto(hptselpsivtxpz, outputstorage);
1112 SaveAndUpdateHisto(hptselpsivtxdpt, outputstorage);
1113 SaveAndUpdateHisto(hptselpsivtxdpz, outputstorage);
1114 
1115 SaveAndUpdateHisto(hvtxpsichi2, outputstorage);
1116 SaveAndUpdateHisto(hvtxpsiacceptmass, outputstorage);
1117 SaveAndUpdateHisto(hvtxpsirejectmass, outputstorage);
1118 SaveAndUpdateHisto(hvtxpsiacceptpt, outputstorage);
1119 SaveAndUpdateHisto(hvtxpsirejectpt, outputstorage);
1120 SaveAndUpdateHisto(hvtxpsiacceptpz, outputstorage);
1121 SaveAndUpdateHisto(hvtxpsirejectpz, outputstorage);
1122 SaveAndUpdateHisto(hvtxpsiacceptdpt, outputstorage);
1123 SaveAndUpdateHisto(hvtxpsirejectdpt, outputstorage);
1124 SaveAndUpdateHisto(hvtxpsiacceptdpz, outputstorage);
1125 SaveAndUpdateHisto(hvtxpsirejectdpz, outputstorage);
1126 
1127 SaveAndUpdateHisto(hptselvtxpsimass, outputstorage);
1128 SaveAndUpdateHisto(hptselvtxpsipt, outputstorage);
1129 SaveAndUpdateHisto(hptselvtxpsipz, outputstorage);
1130 SaveAndUpdateHisto(hptselvtxpsidpt, outputstorage);
1131 SaveAndUpdateHisto(hptselvtxpsidpz, outputstorage);
1132 SaveAndUpdateHisto(hdpvertex, outputstorage);
1133 SaveAndUpdateHisto(hdpvertexfit, outputstorage);
1134 SaveAndUpdateHisto(hdmvertex, outputstorage);
1135 SaveAndUpdateHisto(hdmvertexfit, outputstorage);
1136 
1137 SaveAndUpdateHisto(hdmvertexx, outputstorage);
1138 SaveAndUpdateHisto(hdmvertexy, outputstorage);
1139 SaveAndUpdateHisto(hdmvertexz, outputstorage);
1140 SaveAndUpdateHisto(hdmvertexxfit, outputstorage);
1141 SaveAndUpdateHisto(hdmvertexyfit, outputstorage);
1142 SaveAndUpdateHisto(hdmvertexzfit, outputstorage);
1143 SaveAndUpdateHisto(hdpvertexx, outputstorage);
1144 SaveAndUpdateHisto(hdpvertexy, outputstorage);
1145 SaveAndUpdateHisto(hdpvertexz, outputstorage);
1146 SaveAndUpdateHisto(hdpvertexxfit, outputstorage);
1147 SaveAndUpdateHisto(hdpvertexyfit, outputstorage);
1148 SaveAndUpdateHisto(hdpvertexzfit, outputstorage);
1149 
1150 SaveAndUpdateHisto(hdpvertexxmc, outputstorage);
1151 SaveAndUpdateHisto(hdpvertexymc, outputstorage);
1152 SaveAndUpdateHisto(hdpvertexzmc, outputstorage);
1153 SaveAndUpdateHisto(hdpvertexxreco, outputstorage);
1154 SaveAndUpdateHisto(hdpvertexyreco, outputstorage);
1155 SaveAndUpdateHisto(hdpvertexzreco, outputstorage);
1156 SaveAndUpdateHisto(hdmvertexxmc, outputstorage);
1157 SaveAndUpdateHisto(hdmvertexymc, outputstorage);
1158 SaveAndUpdateHisto(hdmvertexzmc, outputstorage);
1159 SaveAndUpdateHisto(hdmvertexxreco, outputstorage);
1160 SaveAndUpdateHisto(hdmvertexyreco, outputstorage);
1161 SaveAndUpdateHisto(hdmvertexzreco, outputstorage);
1162 SaveAndUpdateHisto(hdpdecaylength, outputstorage);
1163 SaveAndUpdateHisto(hdmdecaylength, outputstorage);
1164 
1165 SaveAndUpdateHisto(hdpvertexxpoca, outputstorage);
1166 SaveAndUpdateHisto(hdpvertexypoca, outputstorage);
1167 SaveAndUpdateHisto(hdpvertexzpoca, outputstorage);
1168 SaveAndUpdateHisto(hdmvertexxpoca, outputstorage);
1169 SaveAndUpdateHisto(hdmvertexypoca, outputstorage);
1170 SaveAndUpdateHisto(hdmvertexzpoca, outputstorage);
1171 SaveAndUpdateHisto(hdpvertexxpocareso, outputstorage);
1172 SaveAndUpdateHisto(hdpvertexypocareso, outputstorage);
1173 SaveAndUpdateHisto(hdpvertexzpocareso, outputstorage);
1174 SaveAndUpdateHisto(hdmvertexxpocareso, outputstorage);
1175 SaveAndUpdateHisto(hdmvertexypocareso, outputstorage);
1176 SaveAndUpdateHisto(hdmvertexzpocareso, outputstorage);
1177 
1178 timer.Stop();
1179 Double_t rtime = timer.RealTime();
1180 Double_t ctime = timer.CpuTime();
1181 
1182 printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);
1183 
1184  return 0;
1185 } // end macro
PndMCTrack * mctrack
printf("RealTime=%f seconds, CpuTime=%f seconds\n", rtime, ctime)
Int_t i
Definition: run_full.C:25
Int_t GetPdgCode() const
Definition: PndMCTrack.h:73
void SaveAndUpdateHisto(TH1D *currenthisto, TFile &storagefile)
Double_t
TStopwatch timer
Definition: hit_dirc.C:51
Double_t z
Double_t ctime
Definition: hit_dirc.C:114
Double_t x
Double_t y
TVector3 GetStartVertex() const
Definition: PndMCTrack.h:76
Double_t rtime
Definition: hit_dirc.C:113
void SaveAndUpdateHisto ( TH1D *  currenthisto,
TFile &  storagefile 
)

Definition at line 10 of file run_ana_mertens_evt7.C.

Referenced by run_ana_mertens_evt7().

11 {
12  if (storagefile.Get(currenthisto->GetName()) != 0) {
13  currenthisto->Add((TH1D*)storagefile.Get(currenthisto->GetName()));
14  }
15  cout << currenthisto->GetName() << ": " << currenthisto->GetEntries() << endl;
16  currenthisto->Write();
17 }