ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/PixelOnlyTracking/MyPixelAnalysisRatios.cc
Revision: 1.1
Committed: Sat Jun 5 17:50:09 2010 UTC (14 years, 10 months ago) by dinardo
Content type: text/plain
Branch: MAIN
Log Message:
*** empty log message ***

File Contents

# User Rev Content
1 dinardo 1.1 // #############################################################
2     // # Program to compute efficiencies, purities and resolutions #
3     // # for pixel stand alone tracking and vertexing #
4     // #############################################################
5     // # Author: Mauro Dinardo #
6     // #############################################################
7    
8     #include <TROOT.h>
9     #include <TF1.h>
10     #include <TH1.h>
11     #include <TH2.h>
12     #include <TFile.h>
13    
14     #include <math.h>
15     #include <map>
16     #include <stdlib.h>
17     #include <sstream>
18     #include <string>
19    
20    
21     using namespace std;
22    
23    
24     // ####################
25     // # Histogram TRACKS #
26     // ####################
27    
28     // @@@@@@ Pixel tracks @@@@@@
29     TH1F* H_pix_trk_pt;
30     TH1F* H_pix_trk_eta;
31    
32     TH1F* H_pix_trk_doubleCounter_pt;
33     TH1F* H_pix_trk_doubleCounter_eta;
34    
35     TH2F* H_pix_trk_normchi2_pt_eta[2];
36     TH2F* H_pix_trk_eta_phi;
37     TH2F* H_pix_trk_eta_pt;
38    
39     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
40     // @ Histograms containing ratios @
41     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
42     TH1F* H_pix_trk_doubleCounter_pt_rel;
43     TH1F* H_pix_trk_doubleCounter_eta_rel;
44    
45     // @@@@@@ Tracker tracks @@@@@@
46     TH1F* H_trk_trk_pt;
47     TH1F* H_trk_trk_eta;
48    
49     TH1F* H_trk_trk_doubleCounter_pt;
50     TH1F* H_trk_trk_doubleCounter_eta;
51    
52     TH2F* H_trk_trk_normchi2_pt_eta[2];
53     TH2F* H_trk_trk_eta_phi;
54     TH2F* H_trk_trk_eta_pt;
55    
56     TH2F* H_trk_trk_eta_phi_algo;
57     TH2F* H_trk_trk_eta_pt_algo;
58    
59     TH1F* H_trk_eff_pt_num;
60     TH1F* H_trk_eff_pt_den;
61     TH1F* H_trk_eff_eta_num;
62     TH1F* H_trk_eff_eta_den;
63    
64     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
65     // @ Histograms containing ratios @
66     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
67     TH1F* H_trk_trk_doubleCounter_pt_rel;
68     TH1F* H_trk_trk_doubleCounter_eta_rel;
69     TH1F* H_trk_eff_pt;
70     TH1F* H_trk_eff_eta;
71    
72     // @@@@@@ Common pixel&tracker tracks @@@@@@
73     typedef map<string, map<string, TH1F*> > myMap;
74     myMap MapFit_d0; // Map to store d0 difference [pt, eta]
75     myMap MapFit_dz; // Map to store dz difference [pt, eta]
76     myMap MapFit_pt; // Map to store pt difference [pt, eta]
77     myMap MapFit_d0_pull; // Map to store d0 pulls [pt, eta]
78     myMap MapFit_dz_pull; // Map to store dz pulls [pt, eta]
79     myMap MapFit_pt_pull; // Map to store pt pulls [pt, eta]
80     myMap MapFit_d0err; // Map to store generalTrack d0err [pt, eta]
81     myMap MapFit_dzerr; // Map to store generalTrack dzerr [pt, eta]
82     myMap MapFit_pterr; // Map to store generalTrack pterr [pt, eta]
83    
84     TH2F* H_trked_eta_phi_pur;
85     TH2F* H_trked_eta_phi_eff;
86     TH2F* H_trked_eta_pt_pur;
87     TH2F* H_trked_eta_pt_eff;
88     TH2F* H_trked_normchi2_pt_pur_eta[2];
89     TH2F* H_trked_normchi2_pt_eff_eta[2];
90    
91     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
92     // @ Histograms containing ratios @
93     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
94     TH1F* H_ptres_eta_pt[5];
95     TH1F* H_ptpull_eta_pt[5];
96    
97     TH1F* H_dzres_eta_pt[5];
98     TH1F* H_dzpull_eta_pt[5];
99    
100     TH1F* H_d0res_eta_pt[5];
101     TH1F* H_d0pull_eta_pt[5];
102    
103     TH1F* H_ptres_pt_eta[2];
104     TH1F* H_ptpull_pt_eta[2];
105    
106     TH1F* H_dzres_pt_eta[2];
107     TH1F* H_dzpull_pt_eta[2];
108    
109     TH1F* H_d0res_pt_eta[2];
110     TH1F* H_d0pull_pt_eta[2];
111    
112     TH1F* H_skew_pt_eta[2];
113    
114     TH1F* H_purity_trked_eta;
115     TH1F* H_purity_trked_pt_eta[2];
116     TH1F* H_purity_trked_pt_whole_eta;
117     TH1F* H_purity_trked_phi[3];
118    
119     TH1F* H_efficiency_trked_eta;
120     TH1F* H_efficiency_trked_pt_eta[2];
121     TH1F* H_efficiency_trked_pt_whole_eta;
122     TH1F* H_efficiency_trked_phi[3];
123    
124     TH1F* H_efficiency_trked_eta_algo;
125     TH1F* H_efficiency_trked_pt_eta_algo[2];
126     TH1F* H_efficiency_trked_pt_whole_eta_algo;
127     TH1F* H_efficiency_trked_phi_algo[3];
128    
129     TH1F* H_tk_pur_eff;
130    
131     TH2F* H_purity_trked_normchi2_pt_eta[2];
132     TH2F* H_efficiency_trked_normchi2_pt_eta[2];
133    
134    
135     // ######################
136     // # Histogram VERTICES #
137     // ######################
138    
139     // @@@@@@ Pixel vertices @@@@@@
140    
141     // @@@@@@ Tracker vertices @@@@@@
142    
143     // @@@@@@ Common pixel&tracker vertices @@@@@@
144     TH1F* H_vx_counters;
145    
146     map<string, TH1F*> map_H_vx_x_diff;
147     map<string, TH1F*> map_H_vx_y_diff;
148     map<string, TH1F*> map_H_vx_z_diff;
149     map<string, TH1F*> map_H_vx_x_diff_pull;
150     map<string, TH1F*> map_H_vx_y_diff_pull;
151     map<string, TH1F*> map_H_vx_z_diff_pull;
152     map<string, TH1F*> map_H_vx_x_diff_orig;
153     map<string, TH1F*> map_H_vx_y_diff_orig;
154     map<string, TH1F*> map_H_vx_z_diff_orig;
155     map<string, TH1F*> map_H_vx_x_err2;
156     map<string, TH1F*> map_H_vx_y_err2;
157     map<string, TH1F*> map_H_vx_z_err2;
158    
159     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
160     // @ Histograms containing ratios @
161     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
162     TH1F* H_vx_x_diff_vs_trk;
163     TH1F* H_vx_y_diff_vs_trk;
164     TH1F* H_vx_z_diff_vs_trk;
165     TH1F* H_vx_x_diff_pull_vs_trk;
166     TH1F* H_vx_y_diff_pull_vs_trk;
167     TH1F* H_vx_z_diff_pull_vs_trk;
168     TH1F* H_vx_x_diff_orig_vs_trk;
169     TH1F* H_vx_y_diff_orig_vs_trk;
170     TH1F* H_vx_z_diff_orig_vs_trk;
171     TH1F* H_vx_x_diff_noerr_vs_trk;
172     TH1F* H_vx_y_diff_noerr_vs_trk;
173     TH1F* H_vx_z_diff_noerr_vs_trk;
174     TH1F* H_vx_pur_eff;
175    
176    
177     // @@@@@@ Local variables @@@@@@
178     double EtaThr;
179     double EtaThrBin[3];
180     double PtThr[6];
181     double zRange; // Unit: [um]
182     double d0Range; // Unit: [cm]
183     double dzRange; // Unit: [cm]
184     double ptRange; // Unit: [GeV/c]
185     double d0errRange; // Unit: [um]^2
186     double dzerrRange; // Unit: [um]^2
187     double pterrRange; // Unit: [MeV/c]^2
188     double pullsRange;
189    
190     unsigned int MinEntriesEffPur; // Histogram minimum number of entries for efficiency and purity
191    
192     double MaxEtaTkTrk; // Max eta that a track should have in order to be taken into account
193     double RangePt; // pt bins. Unit: [GeV/c]
194     double PtStep; // pt bins. Unit: [GeV/c]
195     double RangeEta; // eta bins
196     double EtaStep; // eta bins
197     double MinPtTk; // Mini Pt that a track must have in order to be taken into account. Unit: [GeV/c]
198    
199    
200     void InitVariables ()
201     {
202     MinEntriesEffPur = 12;
203    
204     MaxEtaTkTrk = 2.3;
205     RangePt = 8.0;
206     PtStep = 0.25;
207     RangeEta = 2.3;
208     EtaStep = 0.23;
209     MinPtTk = 0.75;
210    
211     cout << "\n@@@@@@ CONFIGURATION PARAMETERS @@@@@@" << endl;
212    
213     cout << "@@@@@@ GENERAL PARAMETERS @@@@@@" << endl;
214     cout << "MinEntriesEffPur --> " << MinEntriesEffPur << endl;
215    
216     cout << "@@@@@@ TRACK PARAMETERS @@@@@@" << endl;
217     cout << "MaxEtaTkTrk --> " << MaxEtaTkTrk << endl;
218     cout << "RangePt --> " << RangePt << endl;
219     cout << "PtStep --> " << PtStep << endl;
220     cout << "RangeEta --> " << RangeEta << endl;
221     cout << "EtaStep --> " << EtaStep << endl;
222     cout << "MinPtTk --> " << MinPtTk << endl;
223    
224     cout << "@@@@@@ VERTEX PARAMETERS @@@@@@" << endl;
225    
226     EtaThr = 1.5;
227     EtaThrBin[0] = 0;
228     EtaThrBin[1] = EtaThr;
229     EtaThrBin[2] = MaxEtaTkTrk;
230    
231     PtThr[0] = MinPtTk;
232     for (int i = 1; i < 4; i++) PtThr[i] = PtThr[i-1] + 0.5;
233     PtThr[4] = MinPtTk;
234     PtThr[5] = RangePt;
235    
236     zRange = 6000.; // Unit: [um]
237    
238     d0Range = 0.6; // Unit: [cm]
239     dzRange = 0.6; // Unit: [cm]
240     ptRange = 2.0; // Unit: [GeV/c]
241     d0errRange = 60000.; // Unit: [um]^2
242     dzerrRange = 60000.; // Unit: [um]^2
243     pterrRange = 4000.; // Unit: [MeV/c]^2
244    
245     pullsRange = 10.;
246     }
247    
248    
249     bool LoadHistos (string datatype, TFile* filename)
250     {
251     if (datatype.compare("Tk") == 0)
252     {
253     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
254     // @ Histograms containing numerators or denominators @
255     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
256     H_pix_trk_pt = (TH1F*)filename->Get("MyProcess/PixelTracks/Pixeltracks pt distribution");
257     H_pix_trk_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/Pixeltracks eta distribution");
258     H_pix_trk_doubleCounter_pt = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter pixeltracks pt distribution");
259     H_pix_trk_doubleCounter_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter pixeltracks eta distribution");
260    
261     H_pix_trk_eta_phi = (TH2F*)filename->Get("MyProcess/PixelTracks/Pixeltracks distribution in eta and phi");
262     H_pix_trk_eta_pt = (TH2F*)filename->Get("MyProcess/PixelTracks/Pixeltracks distribution in eta and pt");
263    
264     H_trk_trk_pt = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks pt distribution");
265     H_trk_trk_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks eta distribution");
266     H_trk_trk_doubleCounter_pt = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter trackertracks pt distribution");
267     H_trk_trk_doubleCounter_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter trackertracks eta distribution");
268    
269     H_trk_trk_eta_phi = (TH2F*)filename->Get("MyProcess/PixelTracks/Trackertracks distribution in eta and phi");
270     H_trk_trk_eta_pt = (TH2F*)filename->Get("MyProcess/PixelTracks/Trackertracks distribution in eta and pt");
271    
272     H_trk_trk_eta_phi_algo = (TH2F*)filename->Get("MyProcess/PixelTracks/Trackertracks distribution in eta and phi algo");
273     H_trk_trk_eta_pt_algo = (TH2F*)filename->Get("MyProcess/PixelTracks/Trackertracks distribution in eta and pt algo");
274    
275     H_trk_eff_pt_num = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs pt numerator");
276     H_trk_eff_pt_den = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs pt denominator");
277     H_trk_eff_eta_num = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs eta numerator");
278     H_trk_eff_eta_den = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs eta denominator");
279    
280     for (int i = 0; i < 2; i++)
281     {
282     stringstream i_str;
283     i_str << i + 1;
284     TString hName;
285    
286     hName = "MyProcess/PixelTracks/Pixeltracks vs chi2DoF and pt cuts" + i_str.str();
287     H_pix_trk_normchi2_pt_eta[i] = (TH2F*)filename->Get(hName);
288    
289     hName = "MyProcess/PixelTracks/Trackertracks vs chi2DoF and pt cuts thr" + i_str.str();
290     H_trk_trk_normchi2_pt_eta[i] = (TH2F*)filename->Get(hName);
291    
292     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks vs chi2DoF and pt cuts pur thr" + i_str.str();
293     H_trked_normchi2_pt_pur_eta[i] = (TH2F*)filename->Get(hName);
294    
295     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks vs chi2DoF and pt cuts eff thr" + i_str.str();
296     H_trked_normchi2_pt_eff_eta[i] = (TH2F*)filename->Get(hName);
297     }
298    
299     H_trked_eta_phi_pur = (TH2F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks distribution in eta and phi pur");
300     H_trked_eta_phi_eff = (TH2F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks distribution in eta and phi eff");
301     H_trked_eta_pt_pur = (TH2F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks distribution in eta and pt pur");
302     H_trked_eta_pt_eff = (TH2F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Matched pixeltracks distribution in eta and pt eff");
303    
304     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
305     // @ Histograms containing ratios @
306     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
307     H_pix_trk_doubleCounter_pt_rel = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter pixeltracks pt relative distribution");
308     H_pix_trk_doubleCounter_eta_rel = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter pixeltracks eta relative distribution");
309     H_pix_trk_doubleCounter_pt_rel->Reset();
310     H_pix_trk_doubleCounter_eta_rel->Reset();
311    
312     H_trk_trk_doubleCounter_pt_rel = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter trackertracks pt relative distribution");
313     H_trk_trk_doubleCounter_eta_rel = (TH1F*)filename->Get("MyProcess/PixelTracks/Doublecounter trackertracks eta relative distribution");
314     H_trk_trk_doubleCounter_pt_rel->Reset();
315     H_trk_trk_doubleCounter_eta_rel->Reset();
316    
317     H_trk_eff_pt = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs pt");
318     H_trk_eff_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/Trackertracks efficiency vs eta");
319     H_trk_eff_pt->Reset();
320     H_trk_eff_eta->Reset();
321    
322     for (int i = 0; i < 2; i++)
323     {
324     stringstream i_str;
325     i_str << i + 1;
326     TString hName;
327    
328     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Purity pixeltracks vs chi2DoF and pt cuts thr" + i_str.str();
329     H_purity_trked_normchi2_pt_eta[i] = (TH2F*)filename->Get(hName);
330     H_purity_trked_normchi2_pt_eta[i]->Reset();
331    
332     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Purity pixeltracks vs pt thr" + i_str.str();
333     H_purity_trked_pt_eta[i] = (TH1F*)filename->Get(hName);
334     H_purity_trked_pt_eta[i]->Reset();
335    
336     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Efficiency pixeltracks vs chi2DoF and pt cuts thr" + i_str.str();
337     H_efficiency_trked_normchi2_pt_eta[i] = (TH2F*)filename->Get(hName);
338     H_efficiency_trked_normchi2_pt_eta[i]->Reset();
339    
340     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Efficiency pixeltracks vs pt thr" + i_str.str();
341     H_efficiency_trked_pt_eta[i] = (TH1F*)filename->Get(hName);
342     H_efficiency_trked_pt_eta[i]->Reset();
343    
344     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Algoefficiency pixeltracks vs pt thr" + i_str.str();
345     H_efficiency_trked_pt_eta_algo[i] = (TH1F*)filename->Get(hName);
346     H_efficiency_trked_pt_eta_algo[i]->Reset();
347     }
348    
349     for (int i = 0; i < 3; i++)
350     {
351     stringstream i_str;
352     i_str << i + 1;
353     TString hName;
354    
355     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Purity pixeltracks vs phi" + i_str.str();
356     H_purity_trked_phi[i] = (TH1F*)filename->Get(hName);
357     H_purity_trked_phi[i]->Reset();
358    
359     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Efficiency pixeltracks vs phi" + i_str.str();
360     H_efficiency_trked_phi[i] = (TH1F*)filename->Get(hName);
361     H_efficiency_trked_phi[i]->Reset();
362    
363     hName = "MyProcess/PixelTracks/CommonPxTkPlots/Algoefficiency pixeltracks vs phi" + i_str.str();
364     H_efficiency_trked_phi_algo[i] = (TH1F*)filename->Get(hName);
365     H_efficiency_trked_phi_algo[i]->Reset();
366     }
367    
368     H_purity_trked_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Purity pixeltracks vs eta");
369     H_purity_trked_pt_whole_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Purity pixeltracks vs pt");
370     H_purity_trked_eta->Reset();
371     H_purity_trked_pt_whole_eta->Reset();
372    
373     H_efficiency_trked_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Efficiency pixeltracks vs eta");
374     H_efficiency_trked_pt_whole_eta = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Efficiency pixeltracks vs pt");
375     H_efficiency_trked_eta->Reset();
376     H_efficiency_trked_pt_whole_eta->Reset();
377    
378     H_efficiency_trked_eta_algo = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Algoefficiency pixeltracks vs eta");
379     H_efficiency_trked_pt_whole_eta_algo = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Algoefficiency pixeltracks vs pt");
380     H_efficiency_trked_eta_algo->Reset();
381     H_efficiency_trked_pt_whole_eta_algo->Reset();
382    
383     H_tk_pur_eff = (TH1F*)filename->Get("MyProcess/PixelTracks/CommonPxTkPlots/Purity Efficiency Tracks");
384     H_tk_pur_eff->Reset();
385    
386     return true;
387     }
388     else if (datatype.compare("Vx") == 0)
389     {
390     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
391     // @ Histograms containing numerators or denominators @
392     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
393     H_vx_counters = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Counters Match Purity Efficiency for Vertices");
394    
395     stringstream String;
396     stringstream HistoName;
397     String.precision(3);
398     HistoName.precision(3);
399    
400     for (int i = 1; i < 100; i++)
401     {
402     String << i;
403    
404     HistoName.str("");
405     HistoName << "MyProcess/VxNtrks/x_diff_Ntrk_" << String.str();
406     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
407     map_H_vx_x_diff[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
408    
409     HistoName.str("");
410     HistoName << "MyProcess/VxNtrks/y_diff_Ntrk_" << String.str();
411     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
412     map_H_vx_y_diff[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
413    
414     HistoName.str("");
415     HistoName << "MyProcess/VxNtrks/z_diff_Ntrk_" << String.str();
416     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
417     map_H_vx_z_diff[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
418    
419     HistoName.str("");
420     HistoName << "MyProcess/VxNtrks/x_diff_pull_Ntrk_" << String.str();
421     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
422     map_H_vx_x_diff_pull[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
423    
424     HistoName.str("");
425     HistoName << "MyProcess/VxNtrks/y_diff_pull_Ntrk_" << String.str();
426     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
427     map_H_vx_y_diff_pull[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
428    
429     HistoName.str("");
430     HistoName << "MyProcess/VxNtrks/z_diff_pull_Ntrk_" << String.str();
431     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
432     map_H_vx_z_diff_pull[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
433    
434     HistoName.str("");
435     HistoName << "MyProcess/VxNtrks/x_diff_orig_Ntrk_" << String.str();
436     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
437     map_H_vx_x_diff_orig[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
438    
439     HistoName.str("");
440     HistoName << "MyProcess/VxNtrks/y_diff_orig_Ntrk_" << String.str();
441     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
442     map_H_vx_y_diff_orig[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
443    
444     HistoName.str("");
445     HistoName << "MyProcess/VxNtrks/z_diff_orig_Ntrk_" << String.str();
446     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
447     map_H_vx_z_diff_orig[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
448    
449     HistoName.str("");
450     HistoName << "MyProcess/VxNtrks/x_err2_Ntrk_" << String.str();
451     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
452     map_H_vx_x_err2[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
453    
454     HistoName.str("");
455     HistoName << "MyProcess/VxNtrks/y_err2_Ntrk_" << String.str();
456     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
457     map_H_vx_y_err2[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
458    
459     HistoName.str("");
460     HistoName << "MyProcess/VxNtrks/z_err2_Ntrk_" << String.str();
461     if ((TH1F*)filename->Get(HistoName.str().c_str()) != NULL)
462     map_H_vx_z_err2[String.str().c_str()] = (TH1F*)filename->Get(HistoName.str().c_str());
463    
464     String.str("");
465     }
466    
467     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
468     // @ Histograms containing ratios @
469     // @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
470     H_vx_pur_eff = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Purity Efficiency Vertices");
471     H_vx_x_diff_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position difference along X vs Ntrks");
472     H_vx_y_diff_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position difference along Y vs Ntrks");
473     H_vx_z_diff_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position difference along Z vs Ntrks");
474     H_vx_x_diff_pull_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position pull along X vs Ntrks");
475     H_vx_y_diff_pull_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position pull along Y vs Ntrks");
476     H_vx_z_diff_pull_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex position pull along Z vs Ntrks");
477     H_vx_x_diff_orig_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Original vertex position difference along X vs Ntrks");
478     H_vx_y_diff_orig_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Original vertex position difference along Y vs Ntrks");
479     H_vx_z_diff_orig_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Original vertex position difference along Z vs Ntrks");
480     H_vx_x_diff_noerr_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex resolution along X vs Ntrks");
481     H_vx_y_diff_noerr_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex resolution along Y vs Ntrks");
482     H_vx_z_diff_noerr_vs_trk = (TH1F*)filename->Get("MyProcess/PixelVertices/CommonPxTkPlots/Vertex resolution along Z vs Ntrks");
483     H_vx_x_diff_vs_trk->Reset();
484     H_vx_y_diff_vs_trk->Reset();
485     H_vx_z_diff_vs_trk->Reset();
486     H_vx_x_diff_pull_vs_trk->Reset();
487     H_vx_y_diff_pull_vs_trk->Reset();
488     H_vx_z_diff_pull_vs_trk->Reset();
489     H_vx_x_diff_orig_vs_trk->Reset();
490     H_vx_y_diff_orig_vs_trk->Reset();
491     H_vx_z_diff_orig_vs_trk->Reset();
492     H_vx_x_diff_noerr_vs_trk->Reset();
493     H_vx_y_diff_noerr_vs_trk->Reset();
494     H_vx_z_diff_noerr_vs_trk->Reset();
495     H_vx_pur_eff->Reset();
496    
497     return true;
498     }
499    
500     return false;
501     }
502    
503    
504     void WriteHistos (string datatype, TFile* filename)
505     {
506     if (datatype.compare("Tk") == 0)
507     {
508     filename->cd("MyProcess/PixelTracks");
509    
510     H_pix_trk_doubleCounter_eta_rel->Write(H_pix_trk_doubleCounter_eta_rel->GetName(), TObject::kWriteDelete);
511     H_pix_trk_doubleCounter_pt_rel->Write(H_pix_trk_doubleCounter_pt_rel->GetName(), TObject::kWriteDelete);
512    
513     H_trk_trk_doubleCounter_eta_rel->Write(H_trk_trk_doubleCounter_eta_rel->GetName(), TObject::kWriteDelete);
514     H_trk_trk_doubleCounter_pt_rel->Write(H_trk_trk_doubleCounter_pt_rel->GetName(), TObject::kWriteDelete);
515    
516     H_trk_eff_pt->Write(H_trk_eff_pt->GetName(), TObject::kWriteDelete);
517     H_trk_eff_eta->Write(H_trk_eff_eta->GetName(), TObject::kWriteDelete);
518    
519     filename->cd("MyProcess/PixelTracks/CommonPxTkPlots");
520    
521     for (int i = 0; i < 2; i++)
522     {
523     H_purity_trked_normchi2_pt_eta[i]->Write(H_purity_trked_normchi2_pt_eta[i]->GetName(), TObject::kWriteDelete);
524     H_purity_trked_pt_eta[i]->Write(H_purity_trked_pt_eta[i]->GetName(), TObject::kWriteDelete);
525    
526     H_efficiency_trked_normchi2_pt_eta[i]->Write(H_efficiency_trked_normchi2_pt_eta[i]->GetName(), TObject::kWriteDelete);
527     H_efficiency_trked_pt_eta[i]->Write(H_efficiency_trked_pt_eta[i]->GetName(), TObject::kWriteDelete);
528     H_efficiency_trked_pt_eta_algo[i]->Write(H_efficiency_trked_pt_eta_algo[i]->GetName(), TObject::kWriteDelete);
529     }
530    
531     for (int i = 0; i < 3; i++)
532     {
533     H_purity_trked_phi[i]->Write(H_purity_trked_phi[i]->GetName(), TObject::kWriteDelete);
534    
535     H_efficiency_trked_phi[i]->Write(H_efficiency_trked_phi[i]->GetName(), TObject::kWriteDelete);
536     H_efficiency_trked_phi_algo[i]->Write(H_efficiency_trked_phi_algo[i]->GetName(), TObject::kWriteDelete);
537     }
538    
539     H_purity_trked_eta->Write(H_purity_trked_eta->GetName(), TObject::kWriteDelete);
540     H_purity_trked_pt_whole_eta->Write(H_purity_trked_pt_whole_eta->GetName(), TObject::kWriteDelete);
541    
542     H_efficiency_trked_eta->Write(H_efficiency_trked_eta->GetName(), TObject::kWriteDelete);
543     H_efficiency_trked_pt_whole_eta->Write(H_efficiency_trked_pt_whole_eta->GetName(), TObject::kWriteDelete);
544    
545     H_efficiency_trked_eta_algo->Write(H_efficiency_trked_eta_algo->GetName(), TObject::kWriteDelete);
546     H_efficiency_trked_pt_whole_eta_algo->Write(H_efficiency_trked_pt_whole_eta_algo->GetName(), TObject::kWriteDelete);
547    
548     H_tk_pur_eff->Write(H_tk_pur_eff->GetName(), TObject::kWriteDelete);
549     }
550     else if (datatype.compare("Vx") == 0)
551     {
552     filename->cd("MyProcess/PixelVertices/CommonPxTkPlots");
553    
554     H_vx_x_diff_vs_trk->Write(H_vx_x_diff_vs_trk->GetName(), TObject::kWriteDelete);
555     H_vx_y_diff_vs_trk->Write(H_vx_y_diff_vs_trk->GetName(), TObject::kWriteDelete);
556     H_vx_z_diff_vs_trk->Write(H_vx_z_diff_vs_trk->GetName(), TObject::kWriteDelete);
557     H_vx_x_diff_pull_vs_trk->Write(H_vx_x_diff_pull_vs_trk->GetName(), TObject::kWriteDelete);
558     H_vx_y_diff_pull_vs_trk->Write(H_vx_y_diff_pull_vs_trk->GetName(), TObject::kWriteDelete);
559     H_vx_z_diff_pull_vs_trk->Write(H_vx_z_diff_pull_vs_trk->GetName(), TObject::kWriteDelete);
560     H_vx_x_diff_orig_vs_trk->Write(H_vx_x_diff_orig_vs_trk->GetName(), TObject::kWriteDelete);
561     H_vx_y_diff_orig_vs_trk->Write(H_vx_y_diff_orig_vs_trk->GetName(), TObject::kWriteDelete);
562     H_vx_z_diff_orig_vs_trk->Write(H_vx_z_diff_orig_vs_trk->GetName(), TObject::kWriteDelete);
563     H_vx_x_diff_noerr_vs_trk->Write(H_vx_x_diff_noerr_vs_trk->GetName(), TObject::kWriteDelete);
564     H_vx_y_diff_noerr_vs_trk->Write(H_vx_y_diff_noerr_vs_trk->GetName(), TObject::kWriteDelete);
565     H_vx_z_diff_noerr_vs_trk->Write(H_vx_z_diff_noerr_vs_trk->GetName(), TObject::kWriteDelete);
566     H_vx_pur_eff->Write(H_vx_pur_eff->GetName(), TObject::kWriteDelete);
567     }
568     }
569    
570    
571     bool FitAlgorithm (TF1* FitFcn, TH1F* Histo)
572     {
573     double constant;
574     double mean;
575     double sigma;
576    
577     constant = Histo->GetMaximum();
578     mean = Histo->GetMean();
579     sigma = Histo->GetRMS();
580    
581     FitFcn->SetParameters(constant, mean, sigma);
582     Histo->Fit("Gauss", "R");
583    
584     if ((isnan(FitFcn->GetParameter(2)) == false) && (isnan(FitFcn->GetParError(2)) == false)) return true;
585     return false;
586     }
587    
588    
589     void FitAndFillMaps (string mapname,
590     double rangemin,
591     double rangemax)
592     {
593     double DeconvolvErrsq;
594     double ErrDeconvolvErrsq;
595     stringstream PtCutStr;
596     stringstream EtaCutStr;
597     PtCutStr.precision(3);
598     EtaCutStr.precision(3);
599     double Val, Err;
600    
601     TF1* Gauss = new TF1("Gauss", "gaus", rangemin, rangemax);
602    
603     if (mapname.compare("d0pulls") == 0) // @@@ d0pulls MAP @@@
604     {
605     TH1F* d0pulls[5];
606     d0pulls[0] = new TH1F("d0pulls_0", "d0pulls_0", 100, -pullsRange/2., pullsRange/2.);
607     d0pulls[1] = new TH1F("d0pulls_1", "d0pulls_1", 100, -pullsRange/2., pullsRange/2.);
608     d0pulls[2] = new TH1F("d0pulls_2", "d0pulls_2", 100, -pullsRange/2., pullsRange/2.);
609     d0pulls[3] = new TH1F("d0pulls_3", "d0pulls_3", 100, -pullsRange/2., pullsRange/2.);
610     d0pulls[4] = new TH1F("d0pulls_4", "d0pulls_4", 100, -pullsRange/2., pullsRange/2.);
611    
612     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
613     {
614     PtCutStr.str("");
615     PtCutStr << PtCut;
616    
617     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
618     {
619     EtaCutStr.str("");
620     EtaCutStr << EtaCut;
621    
622     if ((MapFit_d0_pull.find(PtCutStr.str().c_str()) != MapFit_d0_pull.end()) &&
623     (MapFit_d0_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0_pull[PtCutStr.str().c_str()].end()))
624     {
625     for (int i = 0; i < 2; i++)
626     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
627     d0pulls[i]->Add(MapFit_d0_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
628     }
629     }
630    
631     for (int i = 0; i < 2; i++)
632     {
633     if (FitAlgorithm (Gauss,d0pulls[i]) == true)
634     {
635     H_d0pull_pt_eta[i]->SetBinContent(H_d0pull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParameter(2));
636     H_d0pull_pt_eta[i]->SetBinError(H_d0pull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParError(2));
637     }
638    
639     d0pulls[i]->Reset();
640     }
641     }
642    
643     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
644     {
645     EtaCutStr.str("");
646     EtaCutStr << EtaCut;
647    
648     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
649     {
650     PtCutStr.str("");
651     PtCutStr << PtCut;
652    
653     if ((MapFit_d0_pull.find(PtCutStr.str().c_str()) != MapFit_d0_pull.end()) &&
654     (MapFit_d0_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0_pull[PtCutStr.str().c_str()].end()))
655     {
656     for (int i = 0; i < 5; i++)
657     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
658     d0pulls[i]->Add(MapFit_d0_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
659     }
660     }
661    
662     for (int i = 0; i < 5; i++)
663     {
664     if (FitAlgorithm (Gauss,d0pulls[i]) == true)
665     {
666     H_d0pull_eta_pt[i]->SetBinContent(H_d0pull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParameter(2));
667     H_d0pull_eta_pt[i]->SetBinError(H_d0pull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParError(2));
668     }
669    
670     d0pulls[i]->Reset();
671     }
672     }
673    
674     for (int i = 0; i < 5; i++) delete d0pulls[i];
675     }
676     else if (mapname.compare("d0") == 0) // @@@ d0 MAP @@@
677     {
678     TH1F* d0[5];
679     d0[0] = new TH1F("d0_0", "d0_0", 100, -d0Range/2., d0Range/2.);
680     d0[1] = new TH1F("d0_1", "d0_1", 100, -d0Range/2., d0Range/2.);
681     d0[2] = new TH1F("d0_2", "d0_2", 100, -d0Range/2., d0Range/2.);
682     d0[3] = new TH1F("d0_3", "d0_3", 100, -d0Range/2., d0Range/2.);
683     d0[4] = new TH1F("d0_4", "d0_4", 100, -d0Range/2., d0Range/2.);
684    
685     TH1F* d0err[5];
686     d0err[0] = new TH1F("d0err_0", "d0err_0", 200, 0., d0errRange);
687     d0err[1] = new TH1F("d0err_1", "d0err_1", 200, 0., d0errRange);
688     d0err[2] = new TH1F("d0err_2", "d0err_2", 200, 0., d0errRange);
689     d0err[3] = new TH1F("d0err_3", "d0err_3", 200, 0., d0errRange);
690     d0err[4] = new TH1F("d0err_4", "d0err_4", 200, 0., d0errRange);
691    
692     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
693     {
694     PtCutStr.str("");
695     PtCutStr << PtCut;
696    
697     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
698     {
699     EtaCutStr.str("");
700     EtaCutStr << EtaCut;
701    
702     if ((MapFit_d0.find(PtCutStr.str().c_str()) != MapFit_d0.end()) &&
703     (MapFit_d0[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0[PtCutStr.str().c_str()].end()) &&
704     (MapFit_d0err.find(PtCutStr.str().c_str()) != MapFit_d0err.end()) &&
705     (MapFit_d0err[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0err[PtCutStr.str().c_str()].end()))
706     {
707     for (int i = 0; i < 2; i++)
708     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
709     {
710     d0[i]->Add(MapFit_d0[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
711     d0err[i]->Add(MapFit_d0err[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
712     }
713     }
714     }
715    
716     for (int i = 0; i < 2; i++)
717     {
718     DeconvolvErrsq = d0err[i]->GetMean() / 100000000.;
719     ErrDeconvolvErrsq = d0err[i]->GetRMS() / sqrt(d0err[i]->GetEntries()) / 100000000.;
720    
721     if ((FitAlgorithm (Gauss,d0[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
722     {
723     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
724     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
725    
726     H_d0res_pt_eta[i]->SetBinContent(H_d0res_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Val);
727     H_d0res_pt_eta[i]->SetBinError(H_d0res_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Err);
728     }
729    
730     d0[i]->Reset();
731     d0err[i]->Reset();
732     }
733     }
734    
735     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
736     {
737     EtaCutStr.str("");
738     EtaCutStr << EtaCut;
739    
740     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
741     {
742     PtCutStr.str("");
743     PtCutStr << PtCut;
744    
745     if ((MapFit_d0.find(PtCutStr.str().c_str()) != MapFit_d0.end()) &&
746     (MapFit_d0[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0[PtCutStr.str().c_str()].end()) &&
747     (MapFit_d0err.find(PtCutStr.str().c_str()) != MapFit_d0err.end()) &&
748     (MapFit_d0err[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_d0err[PtCutStr.str().c_str()].end()))
749     {
750     for (int i = 0; i < 5; i++)
751     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
752     {
753     d0[i]->Add(MapFit_d0[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
754     d0err[i]->Add(MapFit_d0err[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
755     }
756     }
757     }
758    
759     for (int i = 0; i < 5; i++)
760     {
761     DeconvolvErrsq = d0err[i]->GetMean() / 100000000.;
762     ErrDeconvolvErrsq = d0err[i]->GetRMS() / sqrt(d0err[i]->GetEntries()) / 100000000.;
763    
764     if ((FitAlgorithm (Gauss,d0[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
765     {
766     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
767     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
768    
769     H_d0res_eta_pt[i]->SetBinContent(H_d0res_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Val);
770     H_d0res_eta_pt[i]->SetBinError(H_d0res_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Err);
771     }
772    
773     d0[i]->Reset();
774     d0err[i]->Reset();
775     }
776     }
777    
778     for (int i = 0; i < 5; i++) { delete d0[i]; delete d0err[i]; }
779     }
780     else if (mapname.compare("dzpulls") == 0) // @@@ dzpulls MAP @@@
781     {
782     TH1F* dzpulls[5];
783     dzpulls[0] = new TH1F("dzpulls_0", "dzpulls_0", 100, -pullsRange/2., pullsRange/2.);
784     dzpulls[1] = new TH1F("dzpulls_1", "dzpulls_1", 100, -pullsRange/2., pullsRange/2.);
785     dzpulls[2] = new TH1F("dzpulls_2", "dzpulls_2", 100, -pullsRange/2., pullsRange/2.);
786     dzpulls[3] = new TH1F("dzpulls_3", "dzpulls_3", 100, -pullsRange/2., pullsRange/2.);
787     dzpulls[4] = new TH1F("dzpulls_4", "dzpulls_4", 100, -pullsRange/2., pullsRange/2.);
788    
789     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
790     {
791     PtCutStr.str("");
792     PtCutStr << PtCut;
793    
794     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
795     {
796     EtaCutStr.str("");
797     EtaCutStr << EtaCut;
798    
799     if ((MapFit_dz_pull.find(PtCutStr.str().c_str()) != MapFit_dz_pull.end()) &&
800     (MapFit_dz_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dz_pull[PtCutStr.str().c_str()].end()))
801     {
802     for (int i = 0; i < 2; i++)
803     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
804     dzpulls[i]->Add(MapFit_dz_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
805     }
806     }
807    
808     for (int i = 0; i < 2; i++)
809     {
810     if (FitAlgorithm (Gauss,dzpulls[i]) == true)
811     {
812     H_dzpull_pt_eta[i]->SetBinContent(H_dzpull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParameter(2));
813     H_dzpull_pt_eta[i]->SetBinError(H_dzpull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParError(2));
814     }
815    
816     dzpulls[i]->Reset();
817     }
818     }
819    
820     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
821     {
822     EtaCutStr.str("");
823     EtaCutStr << EtaCut;
824    
825     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
826     {
827     PtCutStr.str("");
828     PtCutStr << PtCut;
829    
830     if ((MapFit_dz_pull.find(PtCutStr.str().c_str()) != MapFit_dz_pull.end()) &&
831     (MapFit_dz_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dz_pull[PtCutStr.str().c_str()].end()))
832     {
833     for (int i = 0; i < 5; i++)
834     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
835     dzpulls[i]->Add(MapFit_dz_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
836     }
837     }
838    
839     for (int i = 0; i < 5; i++)
840     {
841     if (FitAlgorithm (Gauss,dzpulls[i]) == true)
842     {
843     H_dzpull_eta_pt[i]->SetBinContent(H_dzpull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParameter(2));
844     H_dzpull_eta_pt[i]->SetBinError(H_dzpull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParError(2));
845     }
846    
847     dzpulls[i]->Reset();
848     }
849     }
850    
851     for (int i = 0; i < 5; i++) delete dzpulls[i];
852     }
853     else if (mapname.compare("dz") == 0) // @@@ dz MAP @@@
854     {
855     TH1F* dz[5];
856     dz[0] = new TH1F("dz_0", "dz_0", 100, -dzRange/2., dzRange/2.);
857     dz[1] = new TH1F("dz_1", "dz_1", 100, -dzRange/2., dzRange/2.);
858     dz[2] = new TH1F("dz_2", "dz_2", 100, -dzRange/2., dzRange/2.);
859     dz[3] = new TH1F("dz_3", "dz_3", 100, -dzRange/2., dzRange/2.);
860     dz[4] = new TH1F("dz_4", "dz_4", 100, -dzRange/2., dzRange/2.);
861    
862     TH1F* dzerr[5];
863     dzerr[0] = new TH1F("dzerr_0", "dzerr_0", 200, 0., dzerrRange);
864     dzerr[1] = new TH1F("dzerr_1", "dzerr_1", 200, 0., dzerrRange);
865     dzerr[2] = new TH1F("dzerr_2", "dzerr_2", 200, 0., dzerrRange);
866     dzerr[3] = new TH1F("dzerr_3", "dzerr_3", 200, 0., dzerrRange);
867     dzerr[4] = new TH1F("dzerr_4", "dzerr_4", 200, 0., dzerrRange);
868    
869     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
870     {
871     PtCutStr.str("");
872     PtCutStr << PtCut;
873    
874     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
875     {
876     EtaCutStr.str("");
877     EtaCutStr << EtaCut;
878    
879     if ((MapFit_dz.find(PtCutStr.str().c_str()) != MapFit_dz.end()) &&
880     (MapFit_dz[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dz[PtCutStr.str().c_str()].end()) &&
881     (MapFit_dzerr.find(PtCutStr.str().c_str()) != MapFit_dzerr.end()) &&
882     (MapFit_dzerr[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dzerr[PtCutStr.str().c_str()].end()))
883     {
884     for (int i = 0; i < 2; i++)
885     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
886     {
887     dz[i]->Add(MapFit_dz[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
888     dzerr[i]->Add(MapFit_dzerr[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
889     }
890     }
891     }
892    
893     for (int i = 0; i < 2; i++)
894     {
895     DeconvolvErrsq = dzerr[i]->GetMean() / 100000000.;
896     ErrDeconvolvErrsq = dzerr[i]->GetRMS() / sqrt(dzerr[i]->GetEntries()) / 100000000.;
897    
898     if ((FitAlgorithm (Gauss,dz[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
899     {
900     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
901     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
902    
903     H_dzres_pt_eta[i]->SetBinContent(H_dzres_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Val);
904     H_dzres_pt_eta[i]->SetBinError(H_dzres_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Err);
905     }
906    
907     dz[i]->Reset();
908     dzerr[i]->Reset();
909     }
910     }
911    
912     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
913     {
914     EtaCutStr.str("");
915     EtaCutStr << EtaCut;
916    
917     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
918     {
919     PtCutStr.str("");
920     PtCutStr << PtCut;
921    
922     if ((MapFit_dz.find(PtCutStr.str().c_str()) != MapFit_dz.end()) &&
923     (MapFit_dz[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dz[PtCutStr.str().c_str()].end()) &&
924     (MapFit_dzerr.find(PtCutStr.str().c_str()) != MapFit_dzerr.end()) &&
925     (MapFit_dzerr[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_dzerr[PtCutStr.str().c_str()].end()))
926     {
927     for (int i = 0; i < 5; i++)
928     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
929     {
930     dz[i]->Add(MapFit_dz[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
931     dzerr[i]->Add(MapFit_dzerr[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
932     }
933     }
934     }
935    
936     for (int i = 0; i < 5; i++)
937     {
938     DeconvolvErrsq = dzerr[i]->GetMean() / 100000000.;
939     ErrDeconvolvErrsq = dzerr[i]->GetRMS() / sqrt(dzerr[i]->GetEntries()) / 100000000.;
940    
941     if ((FitAlgorithm (Gauss,dz[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
942     {
943     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
944     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
945    
946     H_dzres_eta_pt[i]->SetBinContent(H_dzres_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Val);
947     H_dzres_eta_pt[i]->SetBinError(H_dzres_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Err);
948     }
949    
950     dz[i]->Reset();
951     dzerr[i]->Reset();
952     }
953     }
954    
955     for (int i = 0; i < 5; i++) { delete dz[i]; delete dzerr[i]; }
956     }
957     else if (mapname.compare("ptpulls") == 0) // @@@ ptpulls MAP @@@
958     {
959     TH1F* ptpulls[5];
960     ptpulls[0] = new TH1F("ptpulls_0", "ptpulls_0", 100, -pullsRange/2., pullsRange/2.);
961     ptpulls[1] = new TH1F("ptpulls_1", "ptpulls_1", 100, -pullsRange/2., pullsRange/2.);
962     ptpulls[2] = new TH1F("ptpulls_2", "ptpulls_2", 100, -pullsRange/2., pullsRange/2.);
963     ptpulls[3] = new TH1F("ptpulls_3", "ptpulls_3", 100, -pullsRange/2., pullsRange/2.);
964     ptpulls[4] = new TH1F("ptpulls_4", "ptpulls_4", 100, -pullsRange/2., pullsRange/2.);
965    
966     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
967     {
968     PtCutStr.str("");
969     PtCutStr << PtCut;
970    
971     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
972     {
973     EtaCutStr.str("");
974     EtaCutStr << EtaCut;
975    
976     if ((MapFit_pt_pull.find(PtCutStr.str().c_str()) != MapFit_pt_pull.end()) &&
977     (MapFit_pt_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pt_pull[PtCutStr.str().c_str()].end()))
978     {
979     for (int i = 0; i < 2; i++)
980     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
981     ptpulls[i]->Add(MapFit_pt_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
982     }
983     }
984    
985     if ((PtCut == RangePt) || (H_ptpull_pt_eta[0]->GetXaxis()->FindBin(PtCut-PtStep) != H_ptpull_pt_eta[0]->GetXaxis()->FindBin(PtCut)))
986     {
987     for (int i = 0; i < 2; i++)
988     {
989     if (FitAlgorithm (Gauss,ptpulls[i]) == true)
990     {
991     H_ptpull_pt_eta[i]->SetBinContent(H_ptpull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParameter(2));
992     H_ptpull_pt_eta[i]->SetBinError(H_ptpull_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Gauss->GetParError(2));
993     }
994    
995     ptpulls[i]->Reset();
996     }
997     }
998     }
999    
1000     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
1001     {
1002     EtaCutStr.str("");
1003     EtaCutStr << EtaCut;
1004    
1005     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
1006     {
1007     PtCutStr.str("");
1008     PtCutStr << PtCut;
1009    
1010     if ((MapFit_pt_pull.find(PtCutStr.str().c_str()) != MapFit_pt_pull.end()) &&
1011     (MapFit_pt_pull[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pt_pull[PtCutStr.str().c_str()].end()))
1012     {
1013     for (int i = 0; i < 5; i++)
1014     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
1015     ptpulls[i]->Add(MapFit_pt_pull[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
1016     }
1017     }
1018    
1019     for (int i = 0; i < 5; i++)
1020     {
1021     if (FitAlgorithm (Gauss,ptpulls[i]) == true)
1022     {
1023     H_ptpull_eta_pt[i]->SetBinContent(H_ptpull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParameter(2));
1024     H_ptpull_eta_pt[i]->SetBinError(H_ptpull_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Gauss->GetParError(2));
1025     }
1026    
1027     ptpulls[i]->Reset();
1028     }
1029     }
1030    
1031     for (int i = 0; i < 5; i++) delete ptpulls[i];
1032     }
1033     else if (mapname.compare("pt") == 0) // @@@ pt MAP @@@
1034     {
1035     TH1F* pt[5];
1036     pt[0] = new TH1F("pt_0", "pt_0", 100, -ptRange/2., ptRange/2.);
1037     pt[1] = new TH1F("pt_1", "pt_1", 100, -ptRange/2., ptRange/2.);
1038     pt[2] = new TH1F("pt_2", "pt_2", 100, -ptRange/2., ptRange/2.);
1039     pt[3] = new TH1F("pt_3", "pt_3", 100, -ptRange/2., ptRange/2.);
1040     pt[4] = new TH1F("pt_4", "pt_4", 100, -ptRange/2., ptRange/2.);
1041    
1042     TH1F* pterr[5];
1043     pterr[0] = new TH1F("pterr_0", "pterr_0", 800, 0., pterrRange);
1044     pterr[1] = new TH1F("pterr_1", "pterr_1", 800, 0., pterrRange);
1045     pterr[2] = new TH1F("pterr_2", "pterr_2", 800, 0., pterrRange);
1046     pterr[3] = new TH1F("pterr_3", "pterr_3", 800, 0., pterrRange);
1047     pterr[4] = new TH1F("pterr_4", "pterr_4", 800, 0., pterrRange);
1048    
1049     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
1050     {
1051     PtCutStr.str("");
1052     PtCutStr << PtCut;
1053    
1054     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
1055     {
1056     EtaCutStr.str("");
1057     EtaCutStr << EtaCut;
1058    
1059     if ((MapFit_pt.find(PtCutStr.str().c_str()) != MapFit_pt.end()) &&
1060     (MapFit_pt[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pt[PtCutStr.str().c_str()].end()) &&
1061     (MapFit_pterr.find(PtCutStr.str().c_str()) != MapFit_pterr.end()) &&
1062     (MapFit_pterr[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pterr[PtCutStr.str().c_str()].end()))
1063     {
1064     for (int i = 0; i < 2; i++)
1065     if ((EtaCut > EtaThrBin[i]) && (EtaCut <= EtaThrBin[i+1]))
1066     {
1067     pt[i]->Add(MapFit_pt[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
1068     pterr[i]->Add(MapFit_pterr[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
1069     }
1070     }
1071     }
1072    
1073     if ((PtCut == RangePt) || (H_ptres_pt_eta[0]->GetXaxis()->FindBin(PtCut-PtStep) != H_ptres_pt_eta[0]->GetXaxis()->FindBin(PtCut)))
1074     {
1075     for (int i = 0; i < 2; i++)
1076     {
1077     DeconvolvErrsq = pterr[i]->GetMean() / 1000000.;
1078     ErrDeconvolvErrsq = pterr[i]->GetRMS() / sqrt(pterr[i]->GetEntries()) / 1000000.;
1079    
1080     if ((FitAlgorithm (Gauss,pt[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
1081     {
1082     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
1083     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
1084    
1085     H_ptres_pt_eta[i]->SetBinContent(H_ptres_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Val / (PtCut-PtStep/2.));
1086     H_ptres_pt_eta[i]->SetBinError(H_ptres_pt_eta[i]->GetXaxis()->FindBin(PtCut-PtStep), Err / (PtCut-PtStep/2.));
1087    
1088     H_skew_pt_eta[i]->Fill(PtCut-PtStep, pt[i]->GetSkewness());
1089     }
1090    
1091     pt[i]->Reset();
1092     pterr[i]->Reset();
1093     }
1094     }
1095     }
1096    
1097     for (double EtaCut = EtaStep; EtaCut <= RangeEta; EtaCut += EtaStep)
1098     {
1099     EtaCutStr.str("");
1100     EtaCutStr << EtaCut;
1101    
1102     for (double PtCut = MinPtTk+PtStep; PtCut <= RangePt; PtCut += PtStep)
1103     {
1104     PtCutStr.str("");
1105     PtCutStr << PtCut;
1106    
1107     if ((MapFit_pt.find(PtCutStr.str().c_str()) != MapFit_pt.end()) &&
1108     (MapFit_pt[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pt[PtCutStr.str().c_str()].end()) &&
1109     (MapFit_pterr.find(PtCutStr.str().c_str()) != MapFit_pterr.end()) &&
1110     (MapFit_pterr[PtCutStr.str().c_str()].find(EtaCutStr.str().c_str()) != MapFit_pterr[PtCutStr.str().c_str()].end()))
1111     {
1112     for (int i = 0; i < 5; i++)
1113     if ((PtCut > PtThr[i]) && (PtCut <= PtThr[i+1]))
1114     {
1115     pt[i]->Add(MapFit_pt[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
1116     pterr[i]->Add(MapFit_pterr[PtCutStr.str().c_str()][EtaCutStr.str().c_str()]);
1117     }
1118     }
1119     }
1120    
1121     for (int i = 0; i < 5; i++)
1122     {
1123     DeconvolvErrsq = pterr[i]->GetMean() / 1000000.;
1124     ErrDeconvolvErrsq = pterr[i]->GetRMS() / sqrt(pterr[i]->GetEntries()) / 1000000.;
1125    
1126     if ((FitAlgorithm (Gauss,pt[i]) == true) && (Gauss->GetParameter(2)*Gauss->GetParameter(2) > DeconvolvErrsq))
1127     {
1128     Val = sqrt(Gauss->GetParameter(2)*Gauss->GetParameter(2) - DeconvolvErrsq);
1129     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) + powf(ErrDeconvolvErrsq,2.));
1130    
1131     H_ptres_eta_pt[i]->SetBinContent(H_ptres_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Val / (PtThr[i+1] - (PtThr[i+1]-PtThr[i])/2.));
1132     H_ptres_eta_pt[i]->SetBinError(H_ptres_eta_pt[i]->GetXaxis()->FindBin(EtaCut-EtaStep), Err / (PtThr[i+1] - (PtThr[i+1]-PtThr[i])/2.));
1133     }
1134    
1135     pt[i]->Reset();
1136     pterr[i]->Reset();
1137     }
1138     }
1139    
1140     for (int i = 0; i < 5; i++) { delete pt[i]; delete pterr[i]; }
1141     }
1142    
1143     delete Gauss;
1144     }
1145    
1146    
1147     void MakePlots(string datatype)
1148     {
1149     double SumNum;
1150     double SumDen;
1151    
1152     if (datatype.compare("Vx") == 0)
1153     {
1154     // #################
1155     // ### VERTEXING ###
1156     // #################
1157    
1158     // ########################
1159     // # PURITY & EFFICIENCTY #
1160     // ########################
1161     SumNum = H_vx_counters->GetBinContent(H_vx_counters->GetXaxis()->FindBin("MatchEff"));
1162     SumDen = H_vx_counters->GetBinContent(H_vx_counters->GetXaxis()->FindBin("Strip"));
1163     H_vx_pur_eff->Fill("Efficiency", (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1164     H_vx_pur_eff->SetBinError(H_vx_pur_eff->GetXaxis()->FindBin("Efficiency"), (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1165    
1166     SumNum = H_vx_counters->GetBinContent(H_vx_counters->GetXaxis()->FindBin("MatchPur"));
1167     SumDen = H_vx_counters->GetBinContent(H_vx_counters->GetXaxis()->FindBin("Pixel"));
1168     H_vx_pur_eff->Fill("Purity", (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1169     H_vx_pur_eff->SetBinError(H_vx_pur_eff->GetXaxis()->FindBin("Purity"), (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1170     // }
1171    
1172     // #######################################
1173     // # Vertex Resolutions and Pulls vs Trk #
1174     // #######################################
1175     // if ((SubDirPixVxCom.get() != NULL) && (SubDirPixVxNtrks.get() != NULL))
1176     // {
1177     stringstream String;
1178     double Val,Err;
1179    
1180     TF1* Gauss = new TF1("Gauss", "gaus", -zRange/2., zRange/2.);
1181    
1182     for (map<string,TH1F*>::iterator it = map_H_vx_x_diff.begin(); it != map_H_vx_x_diff.end(); it++)
1183     {
1184     String.str("");
1185     String << atoi(it->first.c_str());
1186    
1187     Gauss->SetRange(-zRange/2., zRange/2.);
1188    
1189     // #############################
1190     // # Compute simple difference #
1191     // #############################
1192     if (map_H_vx_x_diff.find(String.str().c_str()) != map_H_vx_x_diff.end() && (FitAlgorithm (Gauss, map_H_vx_x_diff[String.str().c_str()]) == true))
1193     {
1194     H_vx_x_diff_vs_trk->SetBinContent(H_vx_x_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1195     H_vx_x_diff_vs_trk->SetBinError(H_vx_x_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1196     }
1197     if (map_H_vx_y_diff.find(String.str().c_str()) != map_H_vx_y_diff.end() && (FitAlgorithm (Gauss, map_H_vx_y_diff[String.str().c_str()]) == true))
1198     {
1199     H_vx_y_diff_vs_trk->SetBinContent(H_vx_y_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1200     H_vx_y_diff_vs_trk->SetBinError(H_vx_y_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1201     }
1202     if (map_H_vx_z_diff.find(String.str().c_str()) != map_H_vx_z_diff.end() && (FitAlgorithm (Gauss, map_H_vx_z_diff[String.str().c_str()]) == true))
1203     {
1204     H_vx_z_diff_vs_trk->SetBinContent(H_vx_z_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1205     H_vx_z_diff_vs_trk->SetBinError(H_vx_z_diff_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1206     }
1207    
1208     // #################################################
1209     // # Compute difference deconvolving vx resolution #
1210     // #################################################
1211     if ((map_H_vx_x_err2.find(String.str().c_str()) != map_H_vx_x_err2.end()) &&
1212     (map_H_vx_x_diff.find(String.str().c_str()) != map_H_vx_x_diff.end()) &&
1213     (FitAlgorithm (Gauss, map_H_vx_x_diff[String.str().c_str()]) == true) &&
1214     (powf(Gauss->GetParameter(2),2.) > map_H_vx_x_err2[String.str().c_str()]->GetMean()))
1215     {
1216     Val = sqrt(powf(Gauss->GetParameter(2),2.) - map_H_vx_x_err2[String.str().c_str()]->GetMean());
1217    
1218     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) +
1219     powf(map_H_vx_x_err2[String.str().c_str()]->GetRMS()/sqrt(map_H_vx_x_err2[String.str().c_str()]->GetEntries()),2.));
1220    
1221     H_vx_x_diff_noerr_vs_trk->SetBinContent(H_vx_x_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Val);
1222     H_vx_x_diff_noerr_vs_trk->SetBinError(H_vx_x_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Err);
1223     }
1224     if ((map_H_vx_y_err2.find(String.str().c_str()) != map_H_vx_y_err2.end()) &&
1225     (map_H_vx_y_diff.find(String.str().c_str()) != map_H_vx_y_diff.end()) &&
1226     (FitAlgorithm (Gauss, map_H_vx_y_diff[String.str().c_str()]) == true) &&
1227     (powf(Gauss->GetParameter(2),2.) > map_H_vx_y_err2[String.str().c_str()]->GetMean()))
1228     {
1229     Val = sqrt(powf(Gauss->GetParameter(2),2.) - map_H_vx_y_err2[String.str().c_str()]->GetMean());
1230    
1231     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) +
1232     powf(map_H_vx_y_err2[String.str().c_str()]->GetRMS()/sqrt(map_H_vx_y_err2[String.str().c_str()]->GetEntries()),2.));
1233    
1234     H_vx_y_diff_noerr_vs_trk->SetBinContent(H_vx_y_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Val);
1235     H_vx_y_diff_noerr_vs_trk->SetBinError(H_vx_y_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Err);
1236     }
1237     if ((map_H_vx_z_err2.find(String.str().c_str()) != map_H_vx_z_err2.end()) &&
1238     (map_H_vx_z_diff.find(String.str().c_str()) != map_H_vx_z_diff.end()) &&
1239     (FitAlgorithm (Gauss, map_H_vx_z_diff[String.str().c_str()]) == true) &&
1240     (powf(Gauss->GetParameter(2),2.) > map_H_vx_z_err2[String.str().c_str()]->GetMean()))
1241     {
1242     Val = sqrt(powf(Gauss->GetParameter(2),2.) - map_H_vx_z_err2[String.str().c_str()]->GetMean());
1243    
1244     Err = 1./(2.*Val) * sqrt(powf(2.*Gauss->GetParameter(2)*Gauss->GetParError(2),2.) +
1245     powf(map_H_vx_z_err2[String.str().c_str()]->GetRMS()/sqrt(map_H_vx_z_err2[String.str().c_str()]->GetEntries()),2.));
1246    
1247     H_vx_z_diff_noerr_vs_trk->SetBinContent(H_vx_z_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Val);
1248     H_vx_z_diff_noerr_vs_trk->SetBinError(H_vx_z_diff_noerr_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Err);
1249     }
1250    
1251     Gauss->SetRange(-pullsRange/2., pullsRange/2.);
1252    
1253     // #################
1254     // # Compute pulls #
1255     // #################
1256     if (map_H_vx_x_diff_pull.find(String.str().c_str()) != map_H_vx_x_diff_pull.end() && (FitAlgorithm (Gauss, map_H_vx_x_diff_pull[String.str().c_str()]) == true))
1257     {
1258     H_vx_x_diff_pull_vs_trk->SetBinContent(H_vx_x_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1259     H_vx_x_diff_pull_vs_trk->SetBinError(H_vx_x_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1260     }
1261     if (map_H_vx_y_diff_pull.find(String.str().c_str()) != map_H_vx_y_diff_pull.end() && (FitAlgorithm (Gauss, map_H_vx_y_diff_pull[String.str().c_str()]) == true))
1262     {
1263     H_vx_y_diff_pull_vs_trk->SetBinContent(H_vx_y_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1264     H_vx_y_diff_pull_vs_trk->SetBinError(H_vx_y_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1265     }
1266     if (map_H_vx_z_diff_pull.find(String.str().c_str()) != map_H_vx_z_diff_pull.end() && (FitAlgorithm (Gauss, map_H_vx_z_diff_pull[String.str().c_str()]) == true))
1267     {
1268     H_vx_z_diff_pull_vs_trk->SetBinContent(H_vx_z_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1269     H_vx_z_diff_pull_vs_trk->SetBinError(H_vx_z_diff_pull_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1270     }
1271     }
1272    
1273     Gauss->SetRange(-zRange/2., zRange/2.);
1274    
1275     for (map<string,TH1F*>::iterator it = map_H_vx_x_diff_orig.begin(); it != map_H_vx_x_diff_orig.end(); it++)
1276     {
1277     String.str("");
1278     String << atoi(it->first.c_str());
1279    
1280     // ####################################################
1281     // # Compute simple difference with original vertices #
1282     // ####################################################
1283     if (map_H_vx_x_diff_orig.find(String.str().c_str()) != map_H_vx_x_diff_orig.end() && (FitAlgorithm (Gauss, map_H_vx_x_diff_orig[String.str().c_str()]) == true))
1284     {
1285     H_vx_x_diff_orig_vs_trk->SetBinContent(H_vx_x_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1286     H_vx_x_diff_orig_vs_trk->SetBinError(H_vx_x_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1287     }
1288     if (map_H_vx_y_diff_orig.find(String.str().c_str()) != map_H_vx_y_diff_orig.end() && (FitAlgorithm (Gauss, map_H_vx_y_diff_orig[String.str().c_str()]) == true))
1289     {
1290     H_vx_y_diff_orig_vs_trk->SetBinContent(H_vx_y_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1291     H_vx_y_diff_orig_vs_trk->SetBinError(H_vx_y_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1292     }
1293     if (map_H_vx_z_diff_orig.find(String.str().c_str()) != map_H_vx_z_diff_orig.end() && (FitAlgorithm (Gauss, map_H_vx_z_diff_orig[String.str().c_str()]) == true))
1294     {
1295     H_vx_z_diff_orig_vs_trk->SetBinContent(H_vx_z_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParameter(2));
1296     H_vx_z_diff_orig_vs_trk->SetBinError(H_vx_z_diff_orig_vs_trk->GetXaxis()->FindBin(atoi(it->first.c_str())), Gauss->GetParError(2));
1297     }
1298     }
1299    
1300     delete Gauss;
1301     }
1302     else if (datatype.compare("Tk") == 0)
1303     {
1304     // ################
1305     // ### TRACKING ###
1306     // ################
1307    
1308     double TotalSumNum;
1309     double TotalSumDen;
1310    
1311     // ################################
1312     // # 3-PIXEL-HIT EFFICIENCY VS PT #
1313     // ################################
1314     for (int i = 1; i <= H_trk_eff_pt->GetNbinsX(); i++)
1315     {
1316     H_trk_eff_pt->SetBinContent(i, (H_trk_eff_pt_den->GetBinContent(i) > H_trk_eff_pt_num->GetBinContent(i)) ?
1317     H_trk_eff_pt_num->GetBinContent(i) / H_trk_eff_pt_den->GetBinContent(i) : 0.);
1318     H_trk_eff_pt->SetBinError(i, (H_trk_eff_pt_den->GetBinContent(i) > H_trk_eff_pt_num->GetBinContent(i)) ?
1319     sqrt((H_trk_eff_pt_num->GetBinContent(i)*(H_trk_eff_pt_den->GetBinContent(i) - H_trk_eff_pt_num->GetBinContent(i))) /
1320     powf(H_trk_eff_pt_den->GetBinContent(i),3.)) : 0.);
1321     }
1322    
1323     // #################################
1324     // # 3-PIXEL-HIT EFFICIENCY VS ETA #
1325     // #################################
1326     for (int i = 1; i <= H_trk_eff_eta->GetNbinsX(); i++)
1327     {
1328     H_trk_eff_eta->SetBinContent(i, (H_trk_eff_eta_den->GetBinContent(i) > H_trk_eff_eta_num->GetBinContent(i)) ?
1329     H_trk_eff_eta_num->GetBinContent(i) / H_trk_eff_eta_den->GetBinContent(i) : 0.);
1330     H_trk_eff_eta->SetBinError(i, (H_trk_eff_eta_den->GetBinContent(i) > H_trk_eff_eta_num->GetBinContent(i)) ?
1331     sqrt((H_trk_eff_eta_num->GetBinContent(i)*(H_trk_eff_eta_den->GetBinContent(i) - H_trk_eff_eta_num->GetBinContent(i))) /
1332     powf(H_trk_eff_eta_den->GetBinContent(i),3.)) : 0.);
1333     }
1334    
1335     // ################################
1336     // # PIXEL-TRACKS DOUBLE COUNTING #
1337     // ################################
1338     for (int i = 1; i <= H_pix_trk_doubleCounter_eta->GetNbinsX(); i++)
1339     {
1340     H_pix_trk_doubleCounter_eta_rel->SetBinContent(i, (H_pix_trk_eta->GetBinContent(i) > H_pix_trk_doubleCounter_eta->GetBinContent(i)) ?
1341     H_pix_trk_doubleCounter_eta->GetBinContent(i) / H_pix_trk_eta->GetBinContent(i) : 0.);
1342     H_pix_trk_doubleCounter_eta_rel->SetBinError(i, (H_pix_trk_eta->GetBinContent(i) > H_pix_trk_doubleCounter_eta->GetBinContent(i)) ?
1343     sqrt((H_pix_trk_doubleCounter_eta->GetBinContent(i)*(H_pix_trk_eta->GetBinContent(i) - H_pix_trk_doubleCounter_eta->GetBinContent(i))) /
1344     powf(H_pix_trk_eta->GetBinContent(i),3.)) : 0.);
1345     }
1346     for (int i = 1; i <= H_pix_trk_doubleCounter_pt->GetNbinsX(); i++)
1347     {
1348     H_pix_trk_doubleCounter_pt_rel->SetBinContent(i, (H_pix_trk_pt->GetBinContent(i) > H_pix_trk_doubleCounter_pt->GetBinContent(i)) ?
1349     H_pix_trk_doubleCounter_pt->GetBinContent(i) / H_pix_trk_pt->GetBinContent(i) : 0.);
1350     H_pix_trk_doubleCounter_pt_rel->SetBinError(i, (H_pix_trk_pt->GetBinContent(i) > H_pix_trk_doubleCounter_pt->GetBinContent(i)) ?
1351     sqrt((H_pix_trk_doubleCounter_pt->GetBinContent(i)*(H_pix_trk_pt->GetBinContent(i) - H_pix_trk_doubleCounter_pt->GetBinContent(i))) /
1352     powf(H_pix_trk_pt->GetBinContent(i),3.)) : 0.);
1353     }
1354    
1355     // ##################################
1356     // # TRACKER-TRACKS DOUBLE COUNTING #
1357     // ##################################
1358     for (int i = 1; i <= H_trk_trk_doubleCounter_eta->GetNbinsX(); i++)
1359     {
1360     H_trk_trk_doubleCounter_eta_rel->SetBinContent(i, (H_trk_trk_eta->GetBinContent(i) > H_trk_trk_doubleCounter_eta->GetBinContent(i)) ?
1361     H_trk_trk_doubleCounter_eta->GetBinContent(i) / H_trk_trk_eta->GetBinContent(i) : 0.);
1362     H_trk_trk_doubleCounter_eta_rel->SetBinError(i, (H_trk_trk_eta->GetBinContent(i) > H_trk_trk_doubleCounter_eta->GetBinContent(i)) ?
1363     sqrt((H_trk_trk_doubleCounter_eta->GetBinContent(i)*(H_trk_trk_eta->GetBinContent(i) - H_trk_trk_doubleCounter_eta->GetBinContent(i))) /
1364     powf(H_trk_trk_eta->GetBinContent(i),3.)) : 0.);
1365     }
1366     for (int i = 1; i <= H_trk_trk_doubleCounter_pt->GetNbinsX(); i++)
1367     {
1368     H_trk_trk_doubleCounter_pt_rel->SetBinContent(i, (H_trk_trk_pt->GetBinContent(i) > H_trk_trk_doubleCounter_pt->GetBinContent(i)) ?
1369     H_trk_trk_doubleCounter_pt->GetBinContent(i) / H_trk_trk_pt->GetBinContent(i) : 0.);
1370     H_trk_trk_doubleCounter_pt_rel->SetBinError(i, (H_trk_trk_pt->GetBinContent(i) > H_trk_trk_doubleCounter_pt->GetBinContent(i)) ?
1371     sqrt((H_trk_trk_doubleCounter_pt->GetBinContent(i)*(H_trk_trk_pt->GetBinContent(i) - H_trk_trk_doubleCounter_pt->GetBinContent(i))) /
1372     powf(H_trk_trk_pt->GetBinContent(i),3.)) : 0.);
1373     }
1374    
1375     // #####################
1376     // # PURITY vs CHI2-Pt #
1377     // #####################
1378     for (int k = 0; k < 2; k++)
1379     for (int i = 1; i <= H_purity_trked_normchi2_pt_eta[k]->GetNbinsX(); i++)
1380     for (int j = 1; j <= H_purity_trked_normchi2_pt_eta[k]->GetNbinsY(); j++)
1381     if (H_pix_trk_normchi2_pt_eta[k]->GetBinContent(i,j) >= MinEntriesEffPur)
1382     H_purity_trked_normchi2_pt_eta[k]->SetBinContent(i,j, H_trked_normchi2_pt_pur_eta[k]->GetBinContent(i,j) / H_pix_trk_normchi2_pt_eta[k]->GetBinContent(i,j));
1383    
1384     // #################
1385     // # PURITY vs ETA #
1386     // #################
1387     TotalSumNum = 0.;
1388     TotalSumDen = 0.;
1389     for (int i = 1; i <= H_pix_trk_eta_phi->GetNbinsX(); i++)
1390     {
1391     SumNum = 0.;
1392     SumDen = 0.;
1393     for (int j = 1; j <= H_pix_trk_eta_phi->GetNbinsY(); j++)
1394     {
1395     SumNum = SumNum + H_trked_eta_phi_pur->GetBinContent(i,j);
1396     SumDen = SumDen + H_pix_trk_eta_phi->GetBinContent(i,j);
1397     }
1398    
1399     TotalSumNum = TotalSumNum + SumNum;
1400     TotalSumDen = TotalSumDen + SumDen;
1401    
1402     H_purity_trked_eta->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1403     H_purity_trked_eta->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1404     }
1405     H_tk_pur_eff->Fill("Purity", (TotalSumDen >= MinEntriesEffPur) ? TotalSumNum / TotalSumDen : 0.);
1406     H_tk_pur_eff->SetBinError(H_tk_pur_eff->GetXaxis()->FindBin("Purity"), (TotalSumDen >= MinEntriesEffPur) ? sqrt((TotalSumNum*(TotalSumDen-TotalSumNum))/powf(TotalSumDen,3.)) : 0.);
1407    
1408     // ################
1409     // # PURITY vs Pt #
1410     // ################
1411     for (int j = H_pix_trk_eta_pt->GetYaxis()->FindBin(MinPtTk); j <= H_pix_trk_eta_pt->GetNbinsY(); j++)
1412     {
1413     SumNum = 0.;
1414     SumDen = 0.;
1415     for (int i = H_pix_trk_eta_pt->GetXaxis()->FindBin(-EtaThr); i <= H_pix_trk_eta_pt->GetXaxis()->FindBin(EtaThr); i++)
1416     {
1417     SumNum = SumNum + H_trked_eta_pt_pur->GetBinContent(i,j);
1418     SumDen = SumDen + H_pix_trk_eta_pt->GetBinContent(i,j);
1419     }
1420     H_purity_trked_pt_eta[0]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1421     H_purity_trked_pt_eta[0]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1422    
1423     SumNum = 0.;
1424     SumDen = 0.;
1425     for (int i = 1; i < H_pix_trk_eta_pt->GetXaxis()->FindBin(-EtaThr); i++)
1426     {
1427     SumNum = SumNum + H_trked_eta_pt_pur->GetBinContent(i,j);
1428     SumDen = SumDen + H_pix_trk_eta_pt->GetBinContent(i,j);
1429     }
1430     for (int i = H_pix_trk_eta_pt->GetXaxis()->FindBin(EtaThr) + 1; i <= H_pix_trk_eta_pt->GetNbinsX(); i++)
1431     {
1432     SumNum = SumNum + H_trked_eta_pt_pur->GetBinContent(i,j);
1433     SumDen = SumDen + H_pix_trk_eta_pt->GetBinContent(i,j);
1434     }
1435     H_purity_trked_pt_eta[1]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1436     H_purity_trked_pt_eta[1]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1437     }
1438    
1439     for (int j = H_pix_trk_eta_pt->GetYaxis()->FindBin(MinPtTk); j <= H_pix_trk_eta_pt->GetNbinsY(); j++)
1440     {
1441     SumNum = 0.;
1442     SumDen = 0.;
1443     for (int i = H_pix_trk_eta_pt->GetXaxis()->FindBin(-MaxEtaTkTrk); i <= H_pix_trk_eta_pt->GetXaxis()->FindBin(MaxEtaTkTrk); i++)
1444     {
1445     SumNum = SumNum + H_trked_eta_pt_pur->GetBinContent(i,j);
1446     SumDen = SumDen + H_pix_trk_eta_pt->GetBinContent(i,j);
1447     }
1448     H_purity_trked_pt_whole_eta->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1449     H_purity_trked_pt_whole_eta->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1450     }
1451    
1452     // #################
1453     // # PURITY vs Phi #
1454     // #################
1455     for (int i = 1; i <= H_pix_trk_eta_phi->GetNbinsY(); i++)
1456     {
1457     SumNum = 0.;
1458     SumDen = 0.;
1459     for (int j = H_pix_trk_eta_phi->GetXaxis()->FindBin(-MaxEtaTkTrk); j <= H_pix_trk_eta_phi->GetXaxis()->FindBin(-EtaThr); j++)
1460     {
1461     SumNum = SumNum + H_trked_eta_phi_pur->GetBinContent(j,i);
1462     SumDen = SumDen + H_pix_trk_eta_phi->GetBinContent(j,i);
1463     }
1464     H_purity_trked_phi[0]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1465     H_purity_trked_phi[0]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1466    
1467     SumNum = 0.;
1468     SumDen = 0.;
1469     for (int j = H_pix_trk_eta_phi->GetXaxis()->FindBin(-EtaThr); j <= H_pix_trk_eta_phi->GetXaxis()->FindBin(EtaThr); j++)
1470     {
1471     SumNum = SumNum + H_trked_eta_phi_pur->GetBinContent(j,i);
1472     SumDen = SumDen + H_pix_trk_eta_phi->GetBinContent(j,i);
1473     }
1474     H_purity_trked_phi[1]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1475     H_purity_trked_phi[1]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1476    
1477     SumNum = 0.;
1478     SumDen = 0.;
1479     for (int j = H_pix_trk_eta_phi->GetXaxis()->FindBin(EtaThr); j <= H_pix_trk_eta_phi->GetXaxis()->FindBin(MaxEtaTkTrk); j++)
1480     {
1481     SumNum = SumNum + H_trked_eta_phi_pur->GetBinContent(j,i);
1482     SumDen = SumDen + H_pix_trk_eta_phi->GetBinContent(j,i);
1483     }
1484     H_purity_trked_phi[2]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1485     H_purity_trked_phi[2]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1486     }
1487    
1488     // #########################
1489     // # EFFICIENCY vs CHI2-Pt #
1490     // #########################
1491     for (int k = 0; k < 2; k++)
1492     for (int i = 1; i <= H_efficiency_trked_normchi2_pt_eta[k]->GetNbinsX(); i++)
1493     for (int j = 1; j <= H_efficiency_trked_normchi2_pt_eta[k]->GetNbinsY(); j++)
1494     if (H_trk_trk_normchi2_pt_eta[k]->GetBinContent(i,j) >= MinEntriesEffPur)
1495     H_efficiency_trked_normchi2_pt_eta[k]->SetBinContent(i,j, H_trked_normchi2_pt_eff_eta[k]->GetBinContent(i,j) / H_trk_trk_normchi2_pt_eta[k]->GetBinContent(i,j));
1496    
1497     // #####################
1498     // # EFFICIENCY vs ETA #
1499     // #####################
1500     TotalSumNum = 0.;
1501     TotalSumDen = 0.;
1502     for (int i = 1; i <= H_trk_trk_eta_phi->GetNbinsX(); i++)
1503     {
1504     SumNum = 0.;
1505     SumDen = 0.;
1506     for (int j = 1; j <= H_trk_trk_eta_phi->GetNbinsY(); j++)
1507     {
1508     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(i,j);
1509     SumDen = SumDen + H_trk_trk_eta_phi->GetBinContent(i,j);
1510     }
1511    
1512     TotalSumNum = TotalSumNum + SumNum;
1513     TotalSumDen = TotalSumDen + SumDen;
1514    
1515     H_efficiency_trked_eta->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1516     H_efficiency_trked_eta->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1517     }
1518     H_tk_pur_eff->Fill("Efficiency", (TotalSumDen >= MinEntriesEffPur) ? TotalSumNum / TotalSumDen : 0.);
1519     H_tk_pur_eff->SetBinError(H_tk_pur_eff->GetXaxis()->FindBin("Efficiency"), (TotalSumDen >= MinEntriesEffPur) ? sqrt((TotalSumNum*(TotalSumDen-TotalSumNum))/powf(TotalSumDen,3.)) : 0.);
1520    
1521     // ##########################
1522     // # ALGO-EFFICIENCY vs ETA #
1523     // ##########################
1524     TotalSumNum = 0.;
1525     TotalSumDen = 0.;
1526     for (int i = 1; i <= H_trk_trk_eta_phi_algo->GetNbinsX(); i++)
1527     {
1528     SumNum = 0.;
1529     SumDen = 0.;
1530     for (int j = 1; j <= H_trk_trk_eta_phi_algo->GetNbinsY(); j++)
1531     {
1532     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(i,j);
1533     SumDen = SumDen + H_trk_trk_eta_phi_algo->GetBinContent(i,j);
1534     }
1535    
1536     TotalSumNum = TotalSumNum + SumNum;
1537     TotalSumDen = TotalSumDen + SumDen;
1538    
1539     H_efficiency_trked_eta_algo->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1540     H_efficiency_trked_eta_algo->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1541     }
1542     H_tk_pur_eff->Fill("Algo-Efficiency", (TotalSumDen >= MinEntriesEffPur) ? TotalSumNum / TotalSumDen : 0.);
1543     H_tk_pur_eff->SetBinError(H_tk_pur_eff->GetXaxis()->FindBin("Algo-Efficiency"), (TotalSumDen >= MinEntriesEffPur) ? sqrt((TotalSumNum*(TotalSumDen-TotalSumNum))/powf(TotalSumDen,3.)) : 0.);
1544    
1545     // ####################
1546     // # EFFICIENCY vs Pt #
1547     // ####################
1548     for (int j = H_trk_trk_eta_pt->GetYaxis()->FindBin(MinPtTk); j <= H_trk_trk_eta_pt->GetNbinsY(); j++)
1549     {
1550     SumNum = 0.;
1551     SumDen = 0.;
1552     for (int i = H_trk_trk_eta_pt->GetXaxis()->FindBin(-EtaThr); i <= H_trk_trk_eta_pt->GetXaxis()->FindBin(EtaThr); i++)
1553     {
1554     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1555     SumDen = SumDen + H_trk_trk_eta_pt->GetBinContent(i,j);
1556     }
1557     H_efficiency_trked_pt_eta[0]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1558     H_efficiency_trked_pt_eta[0]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1559    
1560     SumNum = 0.;
1561     SumDen = 0.;
1562     for (int i = 1; i < H_trk_trk_eta_pt->GetXaxis()->FindBin(-EtaThr); i++)
1563     {
1564     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1565     SumDen = SumDen + H_trk_trk_eta_pt->GetBinContent(i,j);
1566     }
1567     for (int i = H_trk_trk_eta_pt->GetXaxis()->FindBin(EtaThr) + 1; i <= H_trk_trk_eta_pt->GetNbinsX(); i++)
1568     {
1569     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1570     SumDen = SumDen + H_trk_trk_eta_pt->GetBinContent(i,j);
1571     }
1572     H_efficiency_trked_pt_eta[1]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1573     H_efficiency_trked_pt_eta[1]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1574     }
1575    
1576     for (int j = H_trk_trk_eta_pt->GetYaxis()->FindBin(MinPtTk); j <= H_trk_trk_eta_pt->GetNbinsY(); j++)
1577     {
1578     SumNum = 0.;
1579     SumDen = 0.;
1580     for (int i = H_trk_trk_eta_pt->GetXaxis()->FindBin(-MaxEtaTkTrk); i <= H_trk_trk_eta_pt->GetXaxis()->FindBin(MaxEtaTkTrk); i++)
1581     {
1582     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1583     SumDen = SumDen + H_trk_trk_eta_pt->GetBinContent(i,j);
1584     }
1585     H_efficiency_trked_pt_whole_eta->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1586     H_efficiency_trked_pt_whole_eta->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1587     }
1588    
1589     // #########################
1590     // # ALGO-EFFICIENCY vs Pt #
1591     // #########################
1592     for (int j = H_trk_trk_eta_pt_algo->GetYaxis()->FindBin(MinPtTk); j <= H_trk_trk_eta_pt_algo->GetNbinsY(); j++)
1593     {
1594     SumNum = 0.;
1595     SumDen = 0.;
1596     for (int i = H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(-EtaThr); i <= H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(EtaThr); i++)
1597     {
1598     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1599     SumDen = SumDen + H_trk_trk_eta_pt_algo->GetBinContent(i,j);
1600     }
1601     H_efficiency_trked_pt_eta_algo[0]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1602     H_efficiency_trked_pt_eta_algo[0]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1603    
1604     SumNum = 0.;
1605     SumDen = 0.;
1606     for (int i = 1; i < H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(-EtaThr); i++)
1607     {
1608     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1609     SumDen = SumDen + H_trk_trk_eta_pt_algo->GetBinContent(i,j);
1610     }
1611     for (int i = H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(EtaThr) + 1; i <= H_trk_trk_eta_pt_algo->GetNbinsX(); i++)
1612     {
1613     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1614     SumDen = SumDen + H_trk_trk_eta_pt_algo->GetBinContent(i,j);
1615     }
1616     H_efficiency_trked_pt_eta_algo[1]->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1617     H_efficiency_trked_pt_eta_algo[1]->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1618     }
1619    
1620     for (int j = H_trk_trk_eta_pt_algo->GetYaxis()->FindBin(MinPtTk); j <= H_trk_trk_eta_pt_algo->GetNbinsY(); j++)
1621     {
1622     SumNum = 0.;
1623     SumDen = 0.;
1624     for (int i = H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(-MaxEtaTkTrk); i <= H_trk_trk_eta_pt_algo->GetXaxis()->FindBin(MaxEtaTkTrk); i++)
1625     {
1626     SumNum = SumNum + H_trked_eta_pt_eff->GetBinContent(i,j);
1627     SumDen = SumDen + H_trk_trk_eta_pt_algo->GetBinContent(i,j);
1628     }
1629     H_efficiency_trked_pt_whole_eta_algo->SetBinContent(j, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1630     H_efficiency_trked_pt_whole_eta_algo->SetBinError(j, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1631     }
1632    
1633     // #####################
1634     // # EFFICIENCY vs Phi #
1635     // #####################
1636     for (int i = 1; i <= H_trk_trk_eta_phi->GetNbinsY(); i++)
1637     {
1638     SumNum = 0.;
1639     SumDen = 0.;
1640     for (int j = H_trk_trk_eta_phi->GetXaxis()->FindBin(-MaxEtaTkTrk); j <= H_trk_trk_eta_phi->GetXaxis()->FindBin(-EtaThr); j++)
1641     {
1642     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1643     SumDen = SumDen + H_trk_trk_eta_phi->GetBinContent(j,i);
1644     }
1645     H_efficiency_trked_phi[0]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1646     H_efficiency_trked_phi[0]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1647    
1648     SumNum = 0.;
1649     SumDen = 0.;
1650     for (int j = H_trk_trk_eta_phi->GetXaxis()->FindBin(-EtaThr); j <= H_trk_trk_eta_phi->GetXaxis()->FindBin(EtaThr); j++)
1651     {
1652     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1653     SumDen = SumDen + H_trk_trk_eta_phi->GetBinContent(j,i);
1654     }
1655     H_efficiency_trked_phi[1]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1656     H_efficiency_trked_phi[1]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1657    
1658     SumNum = 0.;
1659     SumDen = 0.;
1660     for (int j = H_trk_trk_eta_phi->GetXaxis()->FindBin(EtaThr); j <= H_trk_trk_eta_phi->GetXaxis()->FindBin(MaxEtaTkTrk); j++)
1661     {
1662     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1663     SumDen = SumDen + H_trk_trk_eta_phi->GetBinContent(j,i);
1664     }
1665     H_efficiency_trked_phi[2]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1666     H_efficiency_trked_phi[2]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1667     }
1668    
1669     // ##########################
1670     // # ALGO-EFFICIENCY vs Phi #
1671     // ##########################
1672     for (int i = 1; i <= H_trk_trk_eta_phi_algo->GetNbinsY(); i++)
1673     {
1674     SumNum = 0.;
1675     SumDen = 0.;
1676     for (int j = H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(-MaxEtaTkTrk); j <= H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(-EtaThr); j++)
1677     {
1678     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1679     SumDen = SumDen + H_trk_trk_eta_phi_algo->GetBinContent(j,i);
1680     }
1681     H_efficiency_trked_phi_algo[0]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1682     H_efficiency_trked_phi_algo[0]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1683    
1684     SumNum = 0.;
1685     SumDen = 0.;
1686     for (int j = H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(-EtaThr); j <= H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(EtaThr); j++)
1687     {
1688     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1689     SumDen = SumDen + H_trk_trk_eta_phi_algo->GetBinContent(j,i);
1690     }
1691     H_efficiency_trked_phi_algo[1]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1692     H_efficiency_trked_phi_algo[1]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1693    
1694     SumNum = 0.;
1695     SumDen = 0.;
1696     for (int j = H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(EtaThr); j <= H_trk_trk_eta_phi_algo->GetXaxis()->FindBin(MaxEtaTkTrk); j++)
1697     {
1698     SumNum = SumNum + H_trked_eta_phi_eff->GetBinContent(j,i);
1699     SumDen = SumDen + H_trk_trk_eta_phi_algo->GetBinContent(j,i);
1700     }
1701     H_efficiency_trked_phi_algo[2]->SetBinContent(i, (SumDen >= MinEntriesEffPur) ? SumNum / SumDen : 0.);
1702     H_efficiency_trked_phi_algo[2]->SetBinError(i, (SumDen >= MinEntriesEffPur) ? sqrt((SumNum*(SumDen-SumNum))/powf(SumDen,3.)) : 0.);
1703     }
1704    
1705     FitAndFillMaps ("d0pulls", -pullsRange/2., pullsRange/2.);
1706     FitAndFillMaps ("d0", -d0Range/2., d0Range/2.);
1707     FitAndFillMaps ("dzpulls", -pullsRange/2., pullsRange/2.);
1708     FitAndFillMaps ("dz", -dzRange/2., dzRange/2.);
1709     FitAndFillMaps ("ptpulls", -pullsRange/2., pullsRange/2.);
1710     FitAndFillMaps ("pt", -ptRange/2., ptRange/2.);
1711     }
1712     }
1713    
1714    
1715     int main (int argc, char** argv)
1716     {
1717     stringstream DataType;
1718     stringstream FileName;
1719    
1720     if (argc == 3)
1721     {
1722     FileName << argv[1];
1723     TFile* File = TFile::Open(FileName.str().c_str(),"UPDATE");
1724    
1725     DataType << argv[2];
1726    
1727     if (LoadHistos(DataType.str(), File) == true)
1728     {
1729     InitVariables();
1730     MakePlots(DataType.str());
1731     WriteHistos(DataType.str(), File);
1732     File->Close();
1733    
1734     return 1;
1735     }
1736     else
1737     {
1738     cout << "Error while loading histograms from: " << FileName.str().c_str() << endl;
1739     return 0;
1740     }
1741     }
1742     else
1743     {
1744     cout << "@@@ Synopsis @@@" << endl;
1745     cout << "./MyPixelAnalysisRatios filename.root datatype[Tk or Vx]" << endl;
1746    
1747     return 0;
1748     }
1749     }