ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/IPHCalignment2/analysis/main_old.cpp
Revision: 1.1
Committed: Fri Nov 25 16:33:02 2011 UTC (13 years, 5 months ago) by econte
Branch: MAIN
CVS Tags: TBD2011, HEAD
Log Message: 
new IPHC alignment

File Contents

# User Rev Content
1 econte 1.1 #include <fstream>
2     #include <iostream>
3     #include <math.h>
4     #include <sstream>
5     #include <stdlib.h>
6     #include "TreeProducer.h"
7     #include "TwoBodyDecaySpy.h"
8     #include "histograms.h"
9     #include <TLorentzVector.h>
10    
11     #define FITMODE 0
12    
13    
14     double cos(double a1, double b1, double c1, double a2, double b2, double c2)
15     {
16     double n1 = sqrt(a1*a1+b1*b1+c1*c1);
17     double n2 = sqrt(a2*a2+b2*b2+c2*c2);
18     double product = a1*a2+b1*b2+c1*c2;
19     return product/(n1*n2);
20     }
21    
22    
23     double stringTodouble(const std::string& word)
24     {
25     std::stringstream str;
26     str<<word;
27     double a;
28     str >> a;
29     return a;
30     }
31    
32     double massInvariant(const double& px1, const double& py1, const double& pz1,
33     const double& px2, const double& py2, const double& pz2)
34     {
35     double p1=sqrt(px1*px1+py1*py1+pz1*pz1);
36     double p2=sqrt(px2*px2+py2*py2+pz2*pz2);
37    
38    
39     double m_muon=105.658e-3;
40    
41     double e1 = sqrt(pow(m_muon,2)+pow(p1,2));
42     double e2 = sqrt(pow(m_muon,2)+pow(p2,2));
43     double E2 = pow(e1+e2,2);
44     double P2 = pow(px1+px2,2)
45     + pow(py1+py2,2)
46     + pow(pz1+pz2,2);
47     return sqrt(E2-P2);
48     }
49    
50     double pT(const double& px1, const double& py1)
51     {
52     double PT = sqrt( pow(px1,2)
53     + pow(py1,2));
54     return PT;
55     }
56    
57     double p(const double& px1, const double& py1, const double& pz1,
58     const double& px2, const double& py2, const double& pz2)
59     {
60     double P = sqrt( pow(px1+px2,2)
61     + pow(py1+py2,2)
62     + pow(pz1+pz2,2));
63     return P;
64     }
65    
66    
67     //################################ Eric and Jeremy ############################################
68    
69     /*
70     void LoopRef(TreeProducer* treeRef, HistoCollection& histos)
71     {
72     Long64_t nentries = treeRef->GetEntriesFast();
73    
74     Long64_t nbytes1 = 0, nb1 = 0;
75     Long64_t nbytes2 = 0, nb2 = 0;
76     for (Long64_t jentry=0; jentry<nentries;jentry++)
77     {
78     Long64_t ientry1 = treeRef->LoadTree(jentry);
79     if (ientry1 < 0) break;
80    
81     nb1 = treeRef->GetEntry(jentry); nbytes1 += nb1;
82    
83     double m=massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]);
84     double mgen=massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]);
85     histos.DeltaM_None->Fill((m-mgen)/mgen);
86    
87     histos.mref_histo->Fill( massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
88    
89    
90     if (treeRef->Track_charge[0]==treeRef->Track_charge[1]) std::cout << "aie" << std::endl;
91     if (treeRef->Track_charge[0]<0)
92     {
93     histos.mVSphi_m->Fill(treeRef->Track_phi[0],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
94     }else
95     {
96     histos.mVSphi_p->Fill(treeRef->Track_phi[0],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
97     }
98    
99     if (treeRef->Track_charge[1]<0)
100     {
101     histos.mVSphi_m->Fill(treeRef->Track_phi[1],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
102     }else
103     {
104     histos.mVSphi_p->Fill(treeRef->Track_phi[1],massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
105     }
106    
107    
108     histos.mgen_histo->Fill( massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]) );
109    
110    
111     for (unsigned int i=0;i<2;i++)
112     {
113     unsigned int j=0;
114     double theMax1 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[0],treeRef->TrackGen_py[0],treeRef->TrackGen_pz[0]);
115     double theMax2 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[1],treeRef->TrackGen_py[1],treeRef->TrackGen_pz[1]);
116     if (theMax1>theMax2) j=0; else j=1;
117    
118     double pt = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2));
119     double ptgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2));
120    
121     double p = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2) + pow(treeRef->Track_pz[i],2));
122     double pgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2) +pow(treeRef->TrackGen_pz[j],2));
123    
124     histos.rptref_histo->Fill((ptgen-pt)/ptgen);
125     histos.rpref_histo->Fill((pgen-p)/pgen);
126     }
127    
128    
129     }
130     }
131    
132    
133    
134    
135     void Loop(TreeProducer* tree, HistoCollection& histos)
136     {
137     Long64_t nentries = tree->GetEntriesFast();
138    
139     Long64_t nbytes1 = 0, nb1 = 0;
140     Long64_t nbytes2 = 0, nb2 = 0;
141     for (Long64_t jentry=0; jentry<nentries;jentry++)
142     {
143     Long64_t ientry2 = tree->LoadTree(jentry);
144     if (ientry2 < 0) break;
145    
146     nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
147    
148    
149    
150     double m=massInvariant(tree->Track_px[0], tree->Track_py[0], tree->Track_pz[0], tree->Track_px[1], tree->Track_py[1], tree->Track_pz[1]);
151     double mgen=massInvariant(tree->TrackGen_px[0], tree->TrackGen_py[0], tree->TrackGen_pz[0],tree->TrackGen_px[1], tree->TrackGen_py[1], tree->TrackGen_pz[1]);
152     histos.DeltaM_Refit->Fill((m-mgen)/mgen);
153    
154     histos.m_histo->Fill( massInvariant(tree->Track_px[0], tree->Track_py[0], tree->Track_pz[0],
155     tree->Track_px[1], tree->Track_py[1], tree->Track_pz[1]) );
156    
157    
158    
159     for (unsigned int i=0;i<2;i++)
160     {
161     unsigned int j=0;
162     double theMax1 = cos(tree->Track_px[i],tree->Track_py[i],tree->Track_pz[i],tree->TrackGen_px[0],tree->TrackGen_py[0],tree->TrackGen_pz[0]);
163     double theMax2 = cos(tree->Track_px[i],tree->Track_py[i],tree->Track_pz[i],tree->TrackGen_px[1],tree->TrackGen_py[1],tree->TrackGen_pz[1]);
164     if (theMax1>theMax2) j=0; else j=1;
165    
166     double pt = sqrt(pow(tree->Track_px[i],2)+ pow(tree->Track_py[i],2));
167     double ptgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2));
168     double p = sqrt(pow(tree->Track_px[i],2)+ pow(tree->Track_py[i],2) + pow(tree->Track_pz[i],2));
169     double pgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2) +pow(tree->TrackGen_pz[j],2));
170    
171     histos.rpt_histo->Fill((ptgen-pt)/ptgen);
172     histos.rp_histo->Fill((pgen-p)/pgen);
173     }
174    
175    
176     }
177    
178     }
179    
180    
181    
182    
183     void LoopTBD(TwoBodyDecaySpy* tree, HistoCollection& histos)
184     {
185     Long64_t nentries = tree->GetEntriesFast();
186    
187     Long64_t nbytes1 = 0, nb1 = 0;
188     Long64_t nbytes2 = 0, nb2 = 0;
189     for (Long64_t jentry=0; jentry<nentries;jentry++)
190     {
191     Long64_t ientry2 = tree->LoadTree(jentry);
192     if (ientry2 < 0) break;
193    
194     nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
195    
196     histos.pp_histo->Fill(tree->NewTrack_Eric_p[0]);
197     histos.pm_histo->Fill(tree->NewTrack_Eric_p[1]);
198     histos.pall_histo->Fill(tree->NewTrack_Eric_p[0]);
199     histos.pall_histo->Fill(tree->NewTrack_Eric_p[1]);
200    
201     histos.ep_histo->Fill(tree->NewTrack_Eric_pErr[0]/tree->NewTrack_Eric_p[0]);
202     histos.em_histo->Fill(tree->NewTrack_Eric_pErr[1]/tree->NewTrack_Eric_p[1]);
203     histos.eall_histo->Fill(tree->NewTrack_Eric_pErr[0]/tree->NewTrack_Eric_p[0]);
204     histos.eall_histo->Fill(tree->NewTrack_Eric_pErr[1]/tree->NewTrack_Eric_p[1]);
205     histos.mass_histo->Fill(tree->Zmass);
206    
207     histos.massrec_histo->Fill(massInvariant(tree->NewTrack_px[0],tree->NewTrack_py[0],tree->NewTrack_pz[0],tree->NewTrack_px[1],tree->NewTrack_py[1],tree->NewTrack_pz[1]));
208     histos.mbefore_histo->Fill(massInvariant(tree->OldTrack_px[0],tree->OldTrack_py[0],tree->OldTrack_pz[0],tree->OldTrack_px[1],tree->OldTrack_py[1],tree->OldTrack_pz[1]));
209    
210     }
211    
212     }
213    
214    
215     void LoopCommon(TreeProducer* tree, TwoBodyDecaySpy * treeTBD, HistoCollection& histos)
216     {
217     Long64_t nentries = tree->GetEntriesFast();
218    
219     Long64_t nbytes1 = 0, nb1 = 0;
220     Long64_t nbytes2 = 0, nb2 = 0;
221     for (Long64_t jentry=0; jentry<nentries;jentry++)
222     {
223     Long64_t ientry1 = treeTBD->LoadTree(jentry);
224     Long64_t ientry2 = tree->LoadTree(jentry);
225     if (ientry2 < 0) break;
226     if (ientry1 < 0) break;
227    
228     nb1 = treeTBD->GetEntry(jentry); nbytes1 += nb1;
229     nb2 = tree->GetEntry(jentry); nbytes2 += nb2;
230    
231    
232     double m=massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]);
233     double mgen=massInvariant(tree->TrackGen_px[0], tree->TrackGen_py[0], tree->TrackGen_pz[0],tree->TrackGen_px[1], tree->TrackGen_py[1], tree->TrackGen_pz[1]);
234     histos.DeltaM_TBD->Fill((m-mgen)/mgen);
235    
236    
237     std::vector<int> sign(2,0);
238     for (unsigned int i=0;i<2;i++)
239     {
240     double Max1 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->Track_px[0],tree->Track_py[0],tree->Track_pz[0]);
241     double Max2 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->Track_px[1],tree->Track_py[1],tree->Track_pz[1]);
242     if (Max1>Max2) sign[i]=tree->Track_charge[0]; else sign[i]=tree->Track_charge[1];
243     }
244    
245     TVector3 vec1(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0]);
246     TVector3 vec2(treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]);
247    
248     if (sign[0]<0)
249     {
250     histos.mVSphi_m_TBD->Fill(vec1.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
251     }
252     else
253     {
254     histos.mVSphi_p_TBD->Fill(vec1.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]));
255     }
256    
257     if (sign[1]<0)
258     {
259     histos.mVSphi_m_TBD->Fill(vec2.Phi(),massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
260     }else
261     {
262     histos.mVSphi_p_TBD->Fill(vec2.Phi(), massInvariant(treeTBD->NewTrack_px[0], treeTBD->NewTrack_py[0], treeTBD->NewTrack_pz[0], treeTBD->NewTrack_px[1], treeTBD->NewTrack_py[1], treeTBD->NewTrack_pz[1]) );
263     }
264    
265     for (unsigned int i=0; i<2; i++)
266     {
267     unsigned int j=0;
268     double theMax1 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->TrackGen_px[0],tree->TrackGen_py[0],tree->TrackGen_pz[0]);
269     double theMax2 = cos(treeTBD->NewTrack_px[i],treeTBD->NewTrack_py[i],treeTBD->NewTrack_pz[i],tree->TrackGen_px[1],tree->TrackGen_py[1],tree->TrackGen_pz[1]);
270     if (theMax1>theMax2) j=0; else j=1;
271    
272     double pt = sqrt(pow(treeTBD->NewTrack_px[i],2)+ pow(treeTBD->NewTrack_py[i],2));
273     double ptgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2));
274     double p = sqrt(pow(treeTBD->NewTrack_px[i],2)+ pow(treeTBD->NewTrack_py[i],2) + pow(treeTBD->NewTrack_pz[i],2));
275     double pgen = sqrt(pow(tree->TrackGen_px[j],2) +pow(tree->TrackGen_py[j],2) +pow(tree->TrackGen_pz[j],2));
276    
277     histos.rptTBD_histo->Fill((ptgen-pt)/ptgen);
278     histos.rpTBD_histo->Fill((pgen-p)/pgen);
279     }
280     }
281    
282     }
283     */
284    
285     //############################# Christophe #################################################
286    
287    
288     void LoopRef(TreeProducer* treeRef, HistoCollection& histos)
289     {
290     Long64_t nentries = treeRef->GetEntriesFast();
291    
292     Long64_t nbytes1 = 0, nb1 = 0;
293     Long64_t nbytes2 = 0, nb2 = 0;
294     for (Long64_t jentry=0; jentry<nentries;jentry++)
295     {
296     Long64_t ientry1 = treeRef->LoadTree(jentry);
297     if (ientry1 < 0) break;
298    
299     nb1 = treeRef->GetEntry(jentry); nbytes1 += nb1;
300    
301     histos.mRef_histo->Fill( massInvariant(treeRef->Track_px[0], treeRef->Track_py[0], treeRef->Track_pz[0], treeRef->Track_px[1], treeRef->Track_py[1], treeRef->Track_pz[1]) );
302    
303     histos.mGen_histo->Fill( massInvariant(treeRef->TrackGen_px[0], treeRef->TrackGen_py[0], treeRef->TrackGen_pz[0],treeRef->TrackGen_px[1], treeRef->TrackGen_py[1], treeRef->TrackGen_pz[1]) );
304    
305    
306     for (unsigned int i=0;i<2;i++)
307     {
308     unsigned int j=0;
309     double theMax1 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[0],treeRef->TrackGen_py[0],treeRef->TrackGen_pz[0]);
310     double theMax2 = cos(treeRef->Track_px[i],treeRef->Track_py[i],treeRef->Track_pz[i],treeRef->TrackGen_px[1],treeRef->TrackGen_py[1],treeRef->TrackGen_pz[1]);
311     if (theMax1>theMax2) j=0; else j=1;
312    
313     double pt = pT(treeRef->Track_px[i], treeRef->Track_py[i]);
314     double ptgen = pT(treeRef->TrackGen_px[j], treeRef->TrackGen_py[j]);
315    
316     double p = sqrt(pow(treeRef->Track_px[i],2)+ pow(treeRef->Track_py[i],2) + pow(treeRef->Track_pz[i],2));
317     double pgen = sqrt(pow(treeRef->TrackGen_px[j],2) +pow(treeRef->TrackGen_py[j],2) +pow(treeRef->TrackGen_pz[j],2));
318    
319     histos.rptRef_histo->Fill((ptgen-pt)/ptgen);
320     histos.rpRef_histo->Fill((pgen-p)/pgen);
321     }
322    
323    
324     }
325     }
326    
327    
328     void LoopVtx(TreeProducer* treeVtx, HistoCollection& histos)
329     {
330     Long64_t nentries = treeVtx->GetEntriesFast();
331    
332     Long64_t nbytes1 = 0, nb1 = 0;
333     Long64_t nbytes2 = 0, nb2 = 0;
334     for (Long64_t jentry=0; jentry<nentries;jentry++)
335     {
336     Long64_t ientry1 = treeVtx->LoadTree(jentry);
337     if (ientry1 < 0) break;
338    
339     nb1 = treeVtx->GetEntry(jentry); nbytes1 += nb1;
340    
341     histos.mVtx_histo->Fill( massInvariant(treeVtx->Track_px[0], treeVtx->Track_py[0], treeVtx->Track_pz[0], treeVtx->Track_px[1], treeVtx->Track_py[1], treeVtx->Track_pz[1]) );
342    
343     //histos.mGen_histo->Fill( massInvariant(treeVtx->TrackGen_px[0], treeVtx->TrackGen_py[0], treeVtx->TrackGen_pz[0],treeVtx->TrackGen_px[1], treeVtx->TrackGen_py[1], treeVtx->TrackGen_pz[1]) );
344    
345    
346     for (unsigned int i=0;i<2;i++)
347     {
348     unsigned int j=0;
349     double theMax1 = cos(treeVtx->Track_px[i],treeVtx->Track_py[i],treeVtx->Track_pz[i],treeVtx->TrackGen_px[0],treeVtx->TrackGen_py[0],treeVtx->TrackGen_pz[0]);
350     double theMax2 = cos(treeVtx->Track_px[i],treeVtx->Track_py[i],treeVtx->Track_pz[i],treeVtx->TrackGen_px[1],treeVtx->TrackGen_py[1],treeVtx->TrackGen_pz[1]);
351     if (theMax1>theMax2) j=0; else j=1;
352    
353     double pt = pT(treeVtx->Track_px[i], treeVtx->Track_py[i]);
354     double ptgen = pT(treeVtx->TrackGen_px[j], treeVtx->TrackGen_py[j]);
355    
356     double p = sqrt(pow(treeVtx->Track_px[i],2)+ pow(treeVtx->Track_py[i],2) + pow(treeVtx->Track_pz[i],2));
357     double pgen = sqrt(pow(treeVtx->TrackGen_px[j],2) +pow(treeVtx->TrackGen_py[j],2) +pow(treeVtx->TrackGen_pz[j],2));
358    
359     histos.rptVtx_histo->Fill((ptgen-pt)/ptgen);
360     histos.rpVtx_histo->Fill((pgen-p)/pgen);
361     }
362    
363    
364     }
365     }
366    
367    
368     void LoopMom(TreeProducer* treeMom, HistoCollection& histos)
369     {
370     Long64_t nentries = treeMom->GetEntriesFast();
371    
372     Long64_t nbytes1 = 0, nb1 = 0;
373     Long64_t nbytes2 = 0, nb2 = 0;
374     for (Long64_t jentry=0; jentry<nentries;jentry++)
375     {
376     Long64_t ientry1 = treeMom->LoadTree(jentry);
377     if (ientry1 < 0) break;
378    
379     nb1 = treeMom->GetEntry(jentry); nbytes1 += nb1;
380    
381     histos.mMom_histo->Fill( massInvariant(treeMom->Track_px[0], treeMom->Track_py[0], treeMom->Track_pz[0], treeMom->Track_px[1], treeMom->Track_py[1], treeMom->Track_pz[1]) );
382    
383     //histos.mGen_histo->Fill( massInvariant(treeMom->TrackGen_px[0], treeMom->TrackGen_py[0], treeMom->TrackGen_pz[0],treeMom->TrackGen_px[1], treeMom->TrackGen_py[1], treeMom->TrackGen_pz[1]) );
384    
385    
386     for (unsigned int i=0;i<2;i++)
387     {
388     unsigned int j=0;
389     double theMax1 = cos(treeMom->Track_px[i],treeMom->Track_py[i],treeMom->Track_pz[i],treeMom->TrackGen_px[0],treeMom->TrackGen_py[0],treeMom->TrackGen_pz[0]);
390     double theMax2 = cos(treeMom->Track_px[i],treeMom->Track_py[i],treeMom->Track_pz[i],treeMom->TrackGen_px[1],treeMom->TrackGen_py[1],treeMom->TrackGen_pz[1]);
391     if (theMax1>theMax2) j=0; else j=1;
392    
393     double pt = pT(treeMom->Track_px[i], treeMom->Track_py[i]);
394     double ptgen = pT(treeMom->TrackGen_px[j], treeMom->TrackGen_py[j]);
395    
396     double p = sqrt(pow(treeMom->Track_px[i],2)+ pow(treeMom->Track_py[i],2) + pow(treeMom->Track_pz[i],2));
397     double pgen = sqrt(pow(treeMom->TrackGen_px[j],2) +pow(treeMom->TrackGen_py[j],2) +pow(treeMom->TrackGen_pz[j],2));
398    
399     histos.rptMom_histo->Fill((ptgen-pt)/ptgen);
400     histos.rpMom_histo->Fill((pgen-p)/pgen);
401     }
402    
403    
404     }
405     }
406    
407    
408     void LoopFull(TreeProducer* treeFull, HistoCollection& histos)
409     {
410     Long64_t nentries = treeFull->GetEntriesFast();
411    
412     Long64_t nbytes1 = 0, nb1 = 0;
413     Long64_t nbytes2 = 0, nb2 = 0;
414     for (Long64_t jentry=0; jentry<nentries;jentry++)
415     {
416     Long64_t ientry1 = treeFull->LoadTree(jentry);
417     if (ientry1 < 0) break;
418    
419     nb1 = treeFull->GetEntry(jentry); nbytes1 += nb1;
420    
421     histos.mFull_histo->Fill( massInvariant(treeFull->Track_px[0], treeFull->Track_py[0], treeFull->Track_pz[0], treeFull->Track_px[1], treeFull->Track_py[1], treeFull->Track_pz[1]) );
422    
423     //histos.mGen_histo->Fill( massInvariant(treeFull->TrackGen_px[0], treeFull->TrackGen_py[0], treeFull->TrackGen_pz[0],treeFull->TrackGen_px[1], treeFull->TrackGen_py[1], treeFull->TrackGen_pz[1]) );
424    
425    
426     for (unsigned int i=0;i<2;i++)
427     {
428     unsigned int j=0;
429     double theMax1 = cos(treeFull->Track_px[i],treeFull->Track_py[i],treeFull->Track_pz[i],treeFull->TrackGen_px[0],treeFull->TrackGen_py[0],treeFull->TrackGen_pz[0]);
430     double theMax2 = cos(treeFull->Track_px[i],treeFull->Track_py[i],treeFull->Track_pz[i],treeFull->TrackGen_px[1],treeFull->TrackGen_py[1],treeFull->TrackGen_pz[1]);
431     if (theMax1>theMax2) j=0; else j=1;
432    
433     double pt = pT(treeFull->Track_px[i], treeFull->Track_py[i]);
434     double ptgen = pT(treeFull->TrackGen_px[j], treeFull->TrackGen_py[j]);
435    
436     double p = sqrt(pow(treeFull->Track_px[i],2)+ pow(treeFull->Track_py[i],2) + pow(treeFull->Track_pz[i],2));
437     double pgen = sqrt(pow(treeFull->TrackGen_px[j],2) +pow(treeFull->TrackGen_py[j],2) +pow(treeFull->TrackGen_pz[j],2));
438    
439     histos.rptFull_histo->Fill((ptgen-pt)/ptgen);
440     histos.rpFull_histo->Fill((pgen-p)/pgen);
441     }
442    
443    
444     }
445     }
446    
447    
448     void LoopKinFit(TreeProducer* treeKinFit, HistoCollection& histos)
449     {
450     Long64_t nentries = treeKinFit->GetEntriesFast();
451    
452     Long64_t nbytes1 = 0, nb1 = 0;
453     Long64_t nbytes2 = 0, nb2 = 0;
454     for (Long64_t jentry=0; jentry<nentries;jentry++)
455     {
456     Long64_t ientry1 = treeKinFit->LoadTree(jentry);
457     if (ientry1 < 0) break;
458    
459     nb1 = treeKinFit->GetEntry(jentry); nbytes1 += nb1;
460    
461     histos.mKinFit_histo->Fill( massInvariant(treeKinFit->Track_px[0], treeKinFit->Track_py[0], treeKinFit->Track_pz[0], treeKinFit->Track_px[1], treeKinFit->Track_py[1], treeKinFit->Track_pz[1]) );
462    
463     //histos.mGen_histo->Fill( massInvariant(treeKinFit->TrackGen_px[0], treeKinFit->TrackGen_py[0], treeKinFit->TrackGen_pz[0],treeKinFit->TrackGen_px[1], treeKinFit->TrackGen_py[1], treeKinFit->TrackGen_pz[1]) );
464    
465    
466     for (unsigned int i=0;i<2;i++)
467     {
468     unsigned int j=0;
469     double theMax1 = cos(treeKinFit->Track_px[i],treeKinFit->Track_py[i],treeKinFit->Track_pz[i],treeKinFit->TrackGen_px[0],treeKinFit->TrackGen_py[0],treeKinFit->TrackGen_pz[0]);
470     double theMax2 = cos(treeKinFit->Track_px[i],treeKinFit->Track_py[i],treeKinFit->Track_pz[i],treeKinFit->TrackGen_px[1],treeKinFit->TrackGen_py[1],treeKinFit->TrackGen_pz[1]);
471     if (theMax1>theMax2) j=0; else j=1;
472    
473     double pt = pT(treeKinFit->Track_px[i], treeKinFit->Track_py[i]);
474     double ptgen = pT(treeKinFit->TrackGen_px[j], treeKinFit->TrackGen_py[j]);
475    
476     double p = sqrt(pow(treeKinFit->Track_px[i],2)+ pow(treeKinFit->Track_py[i],2) + pow(treeKinFit->Track_pz[i],2));
477     double pgen = sqrt(pow(treeKinFit->TrackGen_px[j],2) +pow(treeKinFit->TrackGen_py[j],2) +pow(treeKinFit->TrackGen_pz[j],2));
478    
479     histos.rptKinFit_histo->Fill((ptgen-pt)/ptgen);
480     histos.rpKinFit_histo->Fill((pgen-p)/pgen);
481     }
482    
483    
484     }
485     }
486    
487    
488     int init()
489     {
490    
491     return 0;
492     }
493    
494    
495     int main(int argc, char *argv[])
496     {
497     std::cout << "Analysis begins ...." << std::endl;
498     std::cout << "Check argument ..." << std::endl;
499    
500    
501     if (argc!=2)
502     {
503     std::cout << "ERROR : wrong number of arguments" << std::endl;
504     std::cout << "CORRECT SYNTAX : ./analysis input.root" << std::endl;
505     return 1;
506     }
507    
508     std::string filename = argv[1];
509     std::cout << "Input file : " << filename << std::endl;
510    
511     /// First Fit (without constraint)
512     TreeProducer treeRef(filename,"treeproducer1/ttree_NoAlgo");
513     std::cout << "N_entries for treeRef :" << treeRef.GetEntriesFast() << std::endl;
514    
515    
516     /// For Momentum constraint
517     // TwoBodyDecaySpy treeTBDMom(filename,"TwoBodyDecayMomConstraint/TwoBodyDecaySpyMom");
518     // std::cout << "N_entries for treeTBDMom :" << treeTBDMom.GetEntriesFast() << std::endl;
519    
520     TreeProducer treeMom(filename,"treeproducer3/ttree_MomTBD");
521     std::cout << "N_entries for treeMom :" << treeMom.GetEntriesFast() << std::endl;
522    
523    
524    
525     /// For Vertex constraint
526    
527     // TwoBodyDecaySpy treeTBDVtx(filename,"TwoBodyDecayVertexConstraint/TwoBodyDecaySpyVtx");
528     // std::cout << "N_entries for treeTBDVtx :" << treeTBDVtx.GetEntriesFast() << std::endl;
529    
530     TreeProducer treeVtx(filename,"treeproducer5/ttree_VtxTBD");
531     std::cout << "N_entries for treeVtx :" << treeVtx.GetEntriesFast() << std::endl;
532    
533     /// For Full constraint
534     TreeProducer treeFull(filename,"treeproducer2/ttree_FullTBD");
535     std::cout << "N_entries for treeFull :" << treeFull.GetEntriesFast() << std::endl;
536    
537     /// For KinFitter (mass)constraint
538    
539     // TwoBodyDecaySpy treeKinFitSpy(filename,"doConstraintKinFit/KinFitterSpy");
540     // std::cout << "N_entries for treeKinFitSpy :" << treeTBDVtx.GetEntriesFast() << std::endl;
541    
542     TreeProducer treeKinFit(filename,"treeproducer4/ttree_KinFit");
543     std::cout << "N_entries for treeKinFit :" << treeKinFit.GetEntriesFast() << std::endl;
544    
545    
546     std::cout << "Initializing histogram" << std::endl;
547     HistoCollection histos;
548     histos.Initialize();
549    
550     std::cout << "Filling"<<std::endl;
551    
552     LoopRef(&treeRef,histos); // Fill histos related to measurements at First Fit and at Generator Level
553     LoopVtx(&treeVtx,histos);
554     LoopMom(&treeMom,histos);
555     LoopFull(&treeFull,histos);
556     LoopKinFit(&treeKinFit,histos);
557    
558     //LoopTBD(&treeTBD,histos);
559     //LoopCommon(&tree,&treeTBD,histos);
560     std::cout << std::endl;
561     std::cout << "Drawing plots ... "<< std::endl;
562     histos.DrawInit();
563     histos.Finalize();
564    
565     std::cout << "End"<<std::endl;
566     }
567