attenuation_calibration_comparisons_withfb.C
Go to the documentation of this file.
1 //Author: Prabhjot Singh (prabhjot@fnal.gov).
2 //Dated April 19 2016 03:00 am
3 //This program is a template macro to plot and compare
4 //a. uncorrected PE/cm profiles,
5 //b. threshold and shielding corrections,
6 //c. attenution fitting plots for any randomly chosen plane and cell,
7 //d. calibration constants, coeff_exp, attenlength, background and chisq,
8 //e. list and number of calibrated channels for the following cases
9 //1. ND MC,
10 //2. ND Data,
11 //3. FD MC,
12 //4. FD Data,
13 //5. Tricell vs. Trajectory (two cases in one time)
14 //6. Any two cases in one time e.g.
15 // a. different periods in same analysis,
16 // b. new calibration vs. old calibration,
17 // c. one analyis vs. another
18 //7. Many vs. one case, e.g. all 2nd analysis cases vs. 1st analysis e.g. period1, period2, epoch3b, epoch3c vs. 1st analysis.
19 //
20 //8. You have to set the booleans given at the start of the program to select one or more cases.
21 //
22 //Medbh Campbell (medbh.campbell.15@ucl.ac.uk), 25/10/16: Making additions/changes in order to plot the new FD fibre brightness MC results against SA data
23 
24 #include "HistogramAttr.h"
25 #include "Header.h"
26 
27 bool uncorrectedPEplots = false; // set true if you want plots of uncorrected PE/cm vs. W
28 bool thresholdshadowingplots = true; // set true if you want plots of threshold and shadowing correction factors.
29 bool correctedPEplots = true; // set true if you want plots of corrected PE/cm vs. W. That is, attenuation fitting plots and thier ratios.
30 bool calibrationconstantsplots = false; // set it to true if you wans to make plots of calibration constants
31 bool zerosplots = false; // set true if you want plots comparing below threshold hits for each FB
32  // NB only recent threshold files can do this (i.e. 3rd analysis onwards)
33  // Also, it takes 9 thresh files from BEFORE the fit
34 
35 bool nd = false; // set true for the ND comparisons
36 bool fd = true; // set true for the FD comparisons
37 bool mc = true; // set true for the Monte Carlo based comparisons
38 bool data = false; // set true for the Data based comparisons
39 bool tricellvstrajectory = false; // set true for the tricell vs. trajectory comparisons
40 bool onetoone = false; // set true for one vs. one case e.g. point 6 above
41 bool manytoone = false; // set true for many vs. one case e..g. point 7 above
42 bool manytomany = true; // set true when comparing a fibre brightness calibration to a different fibre calibration
43 bool nddatamuoncatcher = false; // set true for the ND data muon catcher plots
44 
45 bool fb = true; // set true to compare FD fibre brightness MC results against SA data
46 
47 int fbcolors [9] = {100, 94, 90, 80, 418, 66, 60, 221, 6}; // if plotting fibre brightness, we need 9 easily distinguished colours
48 //and ideally we don't want to redefine them each time. I've tried to avoid anything violently fluorescent.
49 
50 attenuation_calibration_comparisons_withfb(){
51  gStyle->SetOptStat(" ");
52  //-----------------------------------------------------------------------------------------------------
53  //set of all .root files for uncorrected PE/cm vs. w plots
54  //char* c1 = "/nova/ana/users/prabhjot/gridjob/Work_with_Adam_and_Brian/nd/NDPCHitData_S150224/mainbody_and_muoncatcher/mainboby/ConsolidateViewsMode_false/Fromgrid_aftercosmiccalib/histofiles/calibNDPCHitData_hist_mainbody_ConsolidateViewsMode_false_S150224.root"; // base file
55  //char* c2 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_period2ddactivity.root"; // 2nd file
56 
57  /*TFile * zerosfile0 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb0/tmp_data/old_hadd_running_total.root");
58  TFile * zerosfile1 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb1/tmp_data/old_hadd_running_total.root");
59  TFile * zerosfile2 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb2/tmp_data/old_hadd_running_total.root");
60  TFile * zerosfile3 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb3/tmp_data/old_hadd_running_total.root");
61  TFile * zerosfile4 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb4/tmp_data/old_hadd_running_total.root");
62  TFile * zerosfile5 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb5/tmp_data/old_hadd_running_total.root");
63  TFile * zerosfile6 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb6/tmp_data/old_hadd_running_total.root");
64  TFile * zerosfile7 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb7/tmp_data/old_hadd_running_total.root");
65  TFile * zerosfile8 = new TFile ("/nova/ana/users/mcampbel/calibration_thirdanalysis/fdmc/gain100/thresh_for_zeros_study/fb8/tmp_data/old_hadd_running_total.root");*/
66 
67  //ND COMPARISONS
68 
69  //List of SA Prabhjot's file locations
70 
71  //ndmc
72  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/calib_hist_2ndAna_fulldetector_ndmc_idealcondtions.root";
73  //char* cthreshold1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/threshold_correction.root";
74  //char* cthreshold1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/comparison_macros_template/old/nd/MC/correction_2ndanalysis.root";
75  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/FitAttenuation/Atten_Fit_Results.root";
76  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/ndmc/Ideal_Conditions/FitAttenuation/final/calib_atten_consts_mc_uniq_sort.txt";
77 
78  //nddata period2ddactivity
79  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_period2ddactivity.root";
80  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody/FitAttenuation/Atten_Fit_Results.root";
81  //char* cattenuation_muoncatcher1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root";
82  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Period2_DDActivity/mainbody_plus_muoncatcher_final_results/calib_atten_consts_uniq_sort.txt";
83 
84  //nddata epoch3b ddactivity
85  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_epoch3b.root";
86  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/mainbody/FitAttenuation/Atten_Fit_Results.root";
87  //char* cattenuation_muoncatcher1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root";
88  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3b_DDActivity/mainbody_plus_muoncatcher_final_results/calib_atten_consts_uniq_sort.txt";
89 
90  //nddata epoch3c ddactivity
91  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3c_DDActivity/mainbody/calib_hist_2ndAna_fulldetector_nddata_mainbody_epoch3c.root";
92  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3c_DDActivity/mainbody/FitAttenuation/Atten_Fit_Results.root";
93  //char* cattenuation_muoncatcher1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3c_DDActivity/muoncatcher/FitAttenuation/Atten_Fit_Results.root";
94  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/nddata/Epoch3c_DDActivity/mainbody_plus_muoncatcher_final_results/calib_atten_consts_uniq_sort.txt";
95 
96  //List of Medbh's 2016 miniproduction file locations
97 
98  //ndmc
99  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/calib_hist_hadded.root";
100  //char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/consolidated/threshold_comparison.root";
101  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/attenfit/Atten_Fit_Results.root";
102  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/ndmc/attenfit/calib_atten_consts_uniq_sort.txt";
103 
104  //nddata period1
105  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period1/calib_hist_hadded.root";
106  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period1/Atten_Fit_Results.root";
107  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata_muoncatcher/period1/Atten_Fit_Results.root";
108  //char* cconstants1 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period1/attenfit/final/calib_atten_consts_uniq_sort.txt";
109 
110  //nddata period2
111  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period2/calib_hist_hadded.root";
112  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period2/Atten_Fit_Results.root";
113  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata_muoncatcher/period2/Atten_Fit_Results.root";
114  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period2/attenfit/final/calib_atten_consts_uniq_sort.txt";
115 
116  //nddata period3
117  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period3/calib_hist_hadded.root";
118  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period3/Atten_Fit_Results.root";
119  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata_muoncatcher/period3/Atten_Fit_Results.root";
120  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/period3/attenfit/final/calib_atten_consts_uniq_sort.txt";
121 
122  //nddata period2ddactivity
123  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period2ddactivity/calib_hists_nddata_miniprod_period2dd.root";
124  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period2ddactivity/Atten_Fit_Results.root";
125  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata_muoncatcher/removebeamspill/period2ddactivity/Atten_Fit_Results.root";
126  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period2ddactivity/main_and_muoncatcher/calib_atten_consts_uniq_sort.txt";
127 
128  //nddata period3ddactivity
129  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period3ddactivity/calib_hist_nddata_miniprod_period3dd.root";
130  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period3ddactivity/Atten_Fit_Results.root";
131  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata_muoncatcher/removebeamspill/period3ddactivity/Atten_Fit_Results.root";
132  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/nddata/removebeamspills/period3ddactivity/attenfit/final/calib_atten_consts_uniq_sort.txt";
133 
134  //List of Medbh's third analysis files (Jan/Feb 2017)
135  //ndmc
136  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/ndmc/uncorrected_attenprofs_fbandconsolidatedbyview/calib_hists_thirdanalysis_ndmc_fb.root";
137  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/ndmc/thresh_fb/thresh_fits_thirdanalysis_nd_fb0-8.root";
138  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/ndmc/thresh_nofb/threshana_ta_ndmc_nofb.root";
139  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/ndmc/corrected_attenprofs_mainbody_consolidated/Atten_Fit_Results.root";
140  //char* cattenuation_muoncatcher2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/ndmc/corrected_attenprofs_muoncatcher/Atten_Fit_Results.root";
141 
142  //-----------------------------------------------------------------------------------------------------
143  //FD COMPARISONS
144 
145  //List of Prabhjot's SA file locations
146 
147  //fdmc low gain
148  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain100/calib_hist_2ndAna_fulldetector_fdmc_gain100.root";
149  //char* cthreshold1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/comparison/correction_2ndanalysislowgain.root";
150  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain100/FitAttenuation/Atten_Fit_Results.root";
151  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain100/FitAttenuation/final/calib_atten_consts_mc_uniq_sort.txt";
152 
153  //fdmc high gain
154  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain140/calib_hist_2ndAna_fulldetector_fdmc_gain140.root";
155  //char* cthreshold1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/comparison/correction_2ndanalysishighgain.root";
156  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain140/FitAttenuation/Atten_Fit_Results.root";
157  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fdmc/Gain140/FitAttenuation/final/calib_atten_consts_mc_uniq_sort.txt";
158 
159  //fddata epoch 3b
160  //char* c1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fddata/Epoch3b/calib_hist_2ndAna_fulldetector_fddata_Epoch3b.root";
161  //char* cattenuation1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fddata/Epoch3b/FitAttenuation/Atten_Fit_Final_Results.root";
162  //char* cconstants1 = "/nova/ana/users/prabhjot/AttenuationCalibration_2ndAnalysis/fddata/Epoch3b/FitAttenuation/final/calib_atten_consts_uniq_sort.txt";
163 
164  //List of Medbh's 2016 miniproduction file locations
165 
166  //fdmc low gain
167  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain100/calib_hist_hadded.root";
168  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain100/nofb/threshold_comparison.root";
169  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain100/attenfit/Atten_Fit_Results.root";
170  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain100/attenfit/calib_atten_consts_uniq_sort.txt";
171 
172  //fdmc high gain
173  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain140/calib_hist_hadded.root";
174  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain140/nofb/threshold_comparison.root";
175  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain140/attenfit/Atten_Fit_Results.root";
176  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fdmc_gain140/attenfit/calib_atten_consts_uniq_sort.txt";
177 
178  //fd data epoch 2b
179  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch2b/calib_hist_hadded.root";
180  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch2b/attenfit_dec22/final/Atten_Fit_Final_Results.root";
181  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch2b/attenfit_dec22/final/calib_atten_consts_uniq_sort.txt";
182 
183  //fd data epoch 3b
184  //char* c2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch3b/calib_hists_fddata_miniprod_epoch3b.root";
185  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch3b/attenfit_dec22/final/Atten_Fit_Final_Results.root";
186  //char* cconstants2 = "/nova/ana/users/mcampbel/calibration_miniproduction/fddata/epoch3b/attenfit_dec22/final/calib_atten_consts_uniq_sort.txt";
187 
188  //List of Medbh's third analysis files (Jan/Feb 2017)
189  //fdmc low gain
190  char* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain100/uncorrected_attenprofs_fbandconsolidatedbyview/calib_hists_thirdanalysis_fdmc_gain100.root";
191  char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain100/thresh/thresh_fits_thirdanalysis_fd_gain100_fb0-8.root";
192  char* cattenuation1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain100/corrected_attenprofs/Atten_Fit_Results.root";
193 
194  //fdmc high gain
195  //char* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain140/uncorrected_attenprofs_fbandconsolidatedbyview/calib_hists_thirdanalysis_fdmc_gain140.root";
196  //char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain140/thresh/thresh_fits_thirdanalysis_fd_gain140_fb0-8.root";
197  //char* cattenuation1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/original_mc/fdmc/gain140/corrected_attenprofs/Atten_Fit_Results.root";
198 
199  //updated light levels fdmc low gain
200  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain100/uncorrected_attenprofs/calibhists_updatedlightlevels_fdmc_gain100.root";
201  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain100/thresh/thresh_fits_thirdanalysis_updatedlightlevels_fd_gain100_fb0-8.root";
202  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain100/corrected_attenprofs/Atten_Fit_Results_0.root";
203 
204  //updated light levels fdmc high gain
205  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/uncorrected_attenprofs/calibhists_updatedlightlevels_fdmc_gain140.root";
206  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/thresh/thresh_fits_thirdanalysis_updatedlightlevels_fd_gain140_fb0-8.root";
207  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/corrected_attenprofs/Atten_Fit_Results_0.root";
208 
209  //updated light levels, real conditions fdmc low gain
210  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain100/uncorrected_attenprofs/calibhists_realandupdatedlightlevels_fdmc_gain100.root";
211  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain100/thresh/thresh_fits_thirdanalysis_realandupdatedlightlevels_fd_gain100_fb0-8.root";
212  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain100/corrected_attenprofs/Atten_Fit_Results.root";
213 
214  //updated light levels, real conditions fdmc high gain
215  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain140/uncorrected_attenprofs/calibhists_realandupdatedlightlevels_fdmc_gain140.root";
216  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain140/thresh/thresh_fits_thirdanalysis_realandupdatedlightlevels_fd_gain140_fb0-8.root";
217  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain140/corrected_attenprofs/Atten_Fit_Results.root";
218 
219  //updated light levels, real conditions, TRAJECTORY HITS (not xy) fdmc high gain
220  char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/trajectoryhits/fdmc_gain140/uncorrected_attenprofs/calibhists_trajhits_fdmc_gain140.root";
221  char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/trajectoryhits/fdmc_gain140/thresh/thresh_fits_thirdanalysis_trajhits_fd_gain140_fb0-8.root";
222  char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/trajectoryhits/fdmc_gain140/corrected_attenprofs/Atten_Fit_Results.root";
223 
224  ///nova/ana/users/mcampbel/calibration_thirdanalysis/trajectoryhits/fdmc_gain140/
225  //----------------------------------------------------------------------------------------------------
226 
227  //.root file for uncorrected PE/cm for tricell vs trajectory study
228  char* c3 = "";
229 
230  //---------------------------------------------------------------------------------------------------
231 
232  TString hist1; // base histogram from the base file. You want to compare other histograms wrt base histogram.
233  TString hist2; // histogram from the second file that you want to compare with the base histogram.
234  TString det; // name of the detector.
235  TString mcdata; // labeling for the Monte Carlo or the data .
236  TString analysis1; // labeling for the base file.
237  TString analysis2; // labelling for the second file.
238  TString var; // variable labelling on the x axis. It is W (cm).
239  TString plane; // plane number.
240  TString cell; // cell number.
241 
242  //analysis1 = "2nd_analysis_period2";
243  //analysis1 = "2nd_analysis";
244  analysis1 = "3rd_analysis";
245  analysis2 = "Updated light levels";
246 
247  hist1 = "h_WPE_corr_xy";
248  hist2 = "h_WPE_corr_xy"; //_corr_xy_ for the tricell, _corr_traj_ for the trajectory
249 
250  var = "W";//W or Cell
251  plane = "000";
252  cell = "007";
253 
254  //ND MC two cases comparison
255  if(nd && onetoone && !fb){
256  det = "nd";
257  if(mc) mcdata = "mc";
258  if(data) mcdata = "data";
259 
260  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
261  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
262 
263  if(mc && thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
264 
265  if(data && nddatamuoncatcher && correctedPEplots)corrected_PE(cattenuation_muoncatcher1, cattenuation_muoncatcher2, det, mcdata, analysis1, analysis2, plane, cell);
266  else if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
267 
268  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
269  }//end of if(nd && mc && onetoone)
270 
271  //FD two case comparison
272  if(fd && onetoone && !fb){
273  det = "fd";
274  if(mc) mcdata = "mc";
275  if(data) mcdata = "data";
276 
277  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
278  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
279 
280  if(mc && thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
281 
282  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
283 
284  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
285  }//end of if(fd && mc && onetoone)
286 
287  if(fb && fd){
288  det = "fd";
289  mcdata = "mc";
290 
291  if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
292  if(thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
293  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
294 
295  }
296 
297  if(fb && nd){
298  det = "nd";
299  mcdata = "mc";
300 
301  if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
302  if(thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
303  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
304  }
305 
306  //if(zerosplots) zerosstudy(zerosfile0, zerosfile1, zerosfile2, zerosfile3, zerosfile4, zerosfile5, zerosfile6, zerosfile7, zerosfile8);
307 
308 
309 }//attenuation_calibration_comparisons()
310 /////////////////////////////////////////////////////////////////////////////////
311 
312 //----------------------------------------------------------------------------------------------------------------------
313 uncorrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
314  TFile* f1 = new TFile(c1, "read");
315  TFile* f2 = new TFile(c2, "read");
316  TString viewstr;
317  TString viewStr;
318  for(int view = 0; view < 2; ++view){
319  if(nddatamuoncatcher) viewstr = view ? "muy" : "mux";
320  else if(!nddatamuoncatcher) viewstr = view ? "y" : "x";
321  if(nddatamuoncatcher) viewStr = view ? "muY" : "muX";
322  else if(!nddatamuoncatcher) viewStr = view ? "Y" : "X";
323 
324  TH2F* h2D_1[9];
325  TH2F* h2D_2[9];
326  TProfile* ProfileX_1 [9];
327  TProfile* ProfileX_2 [9];
328 
329  f1->cd();
330  if(!manytomany){
331  if(f1->GetListOfKeys()->Contains("calib"))
332  h2D_1[0] = (TH2F*)f1->Get("calib/"+hist1+"_in"+viewStr);
333  if(f1->GetListOfKeys()->Contains("make")){
334  h2D_1[0] = (TH2F*)f1->Get("make/"+hist1+"_in"+viewStr);
335  std::cout << "h1 = make/" << hist1 << "_in" << viewStr << std::endl;
336  }
337  ProfileX_1[0] = h2D_1[0]->ProfileX();
338  if(det=="nd" && mcdata=="data") ProfileX_1[0]->Rebin();
339  } else {
340  TString histdir;
341  if(f1->GetListOfKeys()->Contains("calib")) histdir = "calib/";
342  else if(f1->GetListOfKeys()->Contains("make")) histdir = "make/";
343  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
344  TString fbStr = Form("_fb%i", fiberbrightness);
345  h2D_1[fiberbrightness] = (TH2F*)f1->Get(histdir+hist1+"_in"+viewStr+fbStr);
346  ProfileX_1[fiberbrightness] = h2D_1[fiberbrightness]->ProfileX();
347  }
348  }
349 
350  f2->cd();
351  if(fb){
352  TString histdir;
353  if(f2->GetListOfKeys()->Contains("calib")) histdir = "calib/";
354  else if(f2->GetListOfKeys()->Contains("make")) histdir = "make/";
355  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
356  TString fbStr = Form("_fb%i", fiberbrightness);
357  h2D_2[fiberbrightness] = (TH2F*)f2->Get(histdir+hist2+"_in"+viewStr+fbStr);
358  ProfileX_2[fiberbrightness] = h2D_2[fiberbrightness]->ProfileX();
359  }
360  } else {
361  if(f2->GetListOfKeys()->Contains("calib"))
362  h2D_2[0] = (TH2F*)f2->Get("calib/"+hist2+"_in"+viewStr);
363  else if(f2->GetListOfKeys()->Contains("make"))
364  h2D_2[0] = (TH2F*)f2->Get("make/"+hist2+"_in"+viewStr);
365  ProfileX_2[0] = h2D_2[0]->ProfileX();
366  if(det=="nd" && mcdata=="data") ProfileX_2[0]->Rebin();
367  }
368 
369  new TCanvas;
370  gPad->SetGrid();
371 
372  HistogramAttr1D(ProfileX_1[0], "W (cm)", "Uncorrected PE/cm", ProfileX_1[0]->GetXaxis()->GetBinCenter(1), ProfileX_1[0]->GetXaxis()->GetBinCenter(ProfileX_1[0]->GetXaxis()->GetNbins()), 0, 60, kBlack);//old limits: 15, 80, kBlue
373  ProfileX_1[0]->Draw();
374  ProfileX_1[0]->SetLineColor(kBlack);
375 
376  if(fb){
377  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
378  ProfileX_2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
379  ProfileX_2[fiberbrightness]->Draw("same");
380  if(manytomany){
381  ProfileX_1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
382  ProfileX_1[fiberbrightness]->SetLineStyle(kDotted);
383  if(fiberbrightness > 0) ProfileX_1[fiberbrightness]->Draw("same");
384  }
385  }
386  } else {
387  ProfileX_2[0]->SetLineColor(kRed);
388  ProfileX_2[0]->Draw("same");
389  }
390 
391  TLegend *legend = new TLegend(.12,.5,.88,.90);//was .12, .62, .64, .90
392  legend->SetBorderSize(0);
393  legend->SetFillStyle(0);
394  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
395  if(!manytomany) legend->AddEntry(ProfileX_1[0], analysis1+" "+viewStr+" view");
396  else {
397  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
398  TString fbStr = Form("FB %i", fiberbrightness);
399  legend->AddEntry(ProfileX_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr);
400  }
401  }
402  if(fb) {
403  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
404  TString fbStr = Form("FB %i", fiberbrightness);
405  legend->AddEntry(ProfileX_2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr);
406  }
407  } else {
408  legend->AddEntry(ProfileX_2[0], analysis2+" "+viewStr+" view");
409  }
410 
411  legend->Draw();
412 
413  gPad->Print("Images/uncorrected_PEcm_vs_w_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
414 
415  if(manytomany){
416  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
417  new TCanvas;
418  ProfileX_1[fiberbrightness]->SetMinimum(0);
419  ProfileX_1[fiberbrightness]->SetMaximum(60);
420  ProfileX_1[fiberbrightness]->SetLineColor(kBlack);
421  ProfileX_1[fiberbrightness]->SetLineStyle(1);
422  ProfileX_1[fiberbrightness]->SetMarkerColor(kBlack);
423  ProfileX_1[fiberbrightness]->Draw();
424  ProfileX_2[fiberbrightness]->SetLineColor(kRed);
425  ProfileX_2[fiberbrightness]->SetMarkerColor(kRed);
426  ProfileX_2[fiberbrightness]->Draw("same");
427  TString fbstr;
428  fbstr.Form("fb%i", fiberbrightness);
429  TLegend * oneFBlegend = new TLegend(0.12, 0.7, 0.5, 0.88);
430  oneFBlegend->SetBorderSize(0);
431  oneFBlegend->AddEntry(ProfileX_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
432  oneFBlegend->AddEntry(ProfileX_2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
433  oneFBlegend->Draw();
434  gPad->SaveAs("Images/uncorrected_PEcm_vs_w_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view_"+fbstr+".pdf");
435  }
436  }
437 
438  //ratios
439  new TCanvas;
440  gPad->SetGrid();
441  TString hrationame[9];
442  TH1D* hratio[9];
443  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
444  TString fbStr;
445  if(fb) fbStr.Form("_fb%i", fiberbrightness);
446  hrationame[fiberbrightness].Form("hratio%s_%s_%s%s", det.Data(), mcdata.Data(), viewstr.Data(), fbStr.Data());
447  hratio[fiberbrightness] = new TH1D(hrationame[fiberbrightness], "",ProfileX_1[0]->GetXaxis()->GetNbins(), ProfileX_1[0]->GetBinCenter(0), ProfileX_1[0]->GetBinCenter(ProfileX_1[0]->GetXaxis()->GetNbins()));
448  }
449 
450  if(fb) {
451  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
452  if(manytomany) Ratiohist(ProfileX_1[fiberbrightness], ProfileX_2[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, "W (cm)", "Ratio", 0.8, 1.2);//old limits: 0.85, 1.1
453  else Ratiohist(ProfileX_1[0], ProfileX_2[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, "W (cm)", "Ratio", 0.7, 1.6);//old limits: 0.85, 1.1
454  }
455  } else {
456  Ratiohist(ProfileX_1[0], ProfileX_2[0], hratio[0], kRed, 24, "W (cm)", "Ratio", 0.85, 1.10);
457  }
458 
459  TLegend *ratiolegend = new TLegend(.12,.62,.64,.90);
460  ratiolegend->SetBorderSize(0);
461  ratiolegend->SetFillStyle(0);
462  //ratiolegend->SetTextSize(0.04);
463  ratiolegend->AddEntry((TObject*)0, det+" "+mcdata," ");
464  if(fb) {
465  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
466  TString fbStr = Form(" FB %i", fiberbrightness);
467  ratiolegend->AddEntry(hratio[fiberbrightness],analysis2+" to "+analysis1+" "+viewStr+" view"+fbStr);
468  }
469  } else {
470  ratiolegend->AddEntry(hratio[0],analysis2+" to "+analysis1+" "+viewStr+" view");
471  }
472  ratiolegend->Draw();
473  gPad->Print("Images/uncorrected_PEcm_vs_w_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
474  }//end of loop over views
475 } // end of uncorrected_PE()
476 //-------------------------------------------------------------------------------------------------------------------------------
477 void uncorrected_PE_tricellvstrajectory(char* c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
478  uncorrected_PE(c1, c1, det, mcdata, analysis1, analysis2, hist1, hist2);
479 }//end of uncorrected_PE_tricellvstrajectory()
480 //-------------------------------------------------------------------------------------------------------------------------------
481 void thresholdshadowing(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var){
482  std::cout << "Begin threshold/shadowing correction files" << std::endl;
483  TFile* f1 = new TFile(c1, "read");
484  TFile* f2 = new TFile(c2, "read");
485  std::cout << "Found both files" << std::endl;
486  int nbrightnesses = 9;
487  for(int view = 0; view < 2; ++view){
488  const TString viewstr = view ? "y" : "x";
489  const TString viewStr = view ? "Y" : "X";
490 
491  TH1D* correction_1[9];
492  TH1D* corrections_to_compare[9];
493  TH1D* correction_2;
494 
495  f1->cd();
496  if(!manytomany){
497  if(var=="W")correction_1[0] = (TH1D*)f1->Get("a"+viewstr+"view");
498  else if(var=="Cell") correction_1[0] = (TH1D*)f1->Get("c"+viewstr+"view");
499  correction_1[0]->SetLineColor(kBlack);
500  std::cout << "Found non-fb histogram" << std::endl;
501  }
502 
503  if(fb){
504  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++) {
505  const TString fbStr;
506  fbStr.Form("_fb%i", fiberbrightness);
507 
508  if(manytomany) {
509  if(var=="W")correction_1[fiberbrightness] = (TH1D*)f1->Get("a"+viewstr+"view"+fbStr);
510  else if(var=="Cell")correction_1[fiberbrightness] = (TH1D*)f1->Get("c"+viewstr+"view"+fbStr);
511  correction_1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
512  correction_1[fiberbrightness]->SetLineStyle(kDotted);
513  }
514 
515  if(var=="W")corrections_to_compare[fiberbrightness] = (TH1D*)f2->Get("a"+viewstr+"view"+fbStr);
516  else if(var=="Cell")corrections_to_compare[fiberbrightness] = (TH1D*)f2->Get("c"+viewstr+"view"+fbStr);
517  corrections_to_compare[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
518  }
519  }
520  else {
521  f2->cd();
522  correction_2 = (TH1D*)f2->Get("a"+viewstr+"view");
523  correction_2->SetLineColor(kRed);
524  }
525  std::cout << "Found fb histograms" << std::endl;
526 
527  TLegend *legend = new TLegend(.5,.62,.88,.90);
528  legend->SetBorderSize(0);
529  legend->SetFillStyle(0);
530  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
531  if(!manytomany) legend->AddEntry(correction_1[0], analysis1+" "+viewStr+" view");
532  else {
533  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
534  const TString fbStr;
535  fbStr.Form("FB %i", fiberbrightness);
536  legend->AddEntry(correction_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr);
537  }
538  }
539  if(fb){
540  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
541  const TString fbStr;
542  fbStr.Form("FB %i", fiberbrightness);
543  legend->AddEntry(corrections_to_compare[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr);
544  }
545  }
546  else legend->AddEntry(correction_2, analysis2+" "+viewStr+" view");
547 
548  new TCanvas;
549  gPad->SetGrid();
550  if(correction_1[0]->FindObject("stats")){
551  TPaveStats* pave1 = (TPaveStats*)correction_1[0]->FindObject("stats");
552  pave1->Delete();
553  }
554  if(!fb && correction_2->FindObject("stats")){
555  TPaveStats* pave2 = (TPaveStats*)correction_2->FindObject("stats");
556  pave2->Delete();
557  }
558  correction_1[0]->SetMinimum(1.1);
559  correction_1[0]->SetMaximum(1.5);
560  correction_1[0]->Draw();
561  if(fb){
562  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
563  corrections_to_compare[fiberbrightness]->Draw("same");
564  if(manytomany && fiberbrightness > 0) correction_1[fiberbrightness]->Draw("same");
565  }
566  }
567  else correction_2->Draw("same");
568  TString labelx;
569  if(var=="W") labelx = "W (cm)";
570  if(var=="Cell") labelx = "Cell";
571 
572  float ylow = 1.0;
573  float yhigh = 1.8;
574  /*if(view == 0) {
575  ylow = 1.2;
576  yhigh = 1.5;
577  }*/
578 
579  //std::cout << "Drawing hist... " << std::endl;
580  //std::cout << "title: " << labelx.Data();
581  //std::cout << " bins:" << correction_1->GetXaxis()->GetBinCenter(1);
582  //std::cout << "-" << correction_1->GetXaxis()->GetBinCenter(correction_1->GetXaxis()->GetNbins());
583  HistogramAttr1D(correction_1[0], labelx.Data(), "correction factor", correction_1[0]->GetXaxis()->GetBinCenter(1), correction_1[0]->GetXaxis()->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()), ylow, yhigh, kBlack);//Old limits: 0.5, 2.0
584  //std::cout << "Done!" << std::endl;
585 
586  legend->Draw();
587  gPad->Print("Images/thresholdshadowing_vs_"+var+"_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
588 
589  if(manytomany){
590  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; ++fiberbrightness){
591  new TCanvas;
592  correction_1[fiberbrightness]->SetMinimum(1.1);
593  correction_1[fiberbrightness]->SetMaximum(1.6);
594  correction_1[fiberbrightness]->Draw();
595  corrections_to_compare[fiberbrightness]->Draw("same");
596  TString fbstr;
597  fbstr.Form("fb%i", fiberbrightness);
598  TLegend * oneFBlegend = new TLegend(0.5, 0.7, 0.9, 0.9);
599  oneFBlegend->AddEntry(correction_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
600  oneFBlegend->AddEntry(corrections_to_compare[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
601  oneFBlegend->Draw();
602  gPad->SaveAs("Images/thresholdshadowing_vs_"+var+"_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view_"+fbstr+".pdf");
603  }
604  }
605 
606  new TCanvas;
607  gPad->SetGrid();
608 
609  //ratios
610  TH1D* hratio[9];
611  if(fb){
612  TLegend *legend = new TLegend(.5,.62,.88,.90);
613  legend->SetBorderSize(0);
614  legend->SetFillStyle(0);
615  //legend->SetTextSize(0.04);
616  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
617 
618  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
619  const TString fbStr = Form("FB %i", fiberbrightness);
620  const TString hrationame = Form("hratio%s_%s_%s_FB%i_%s", det.Data(), mcdata.Data(), viewstr.Data(), fiberbrightness, var.Data());
621  std::cout << hrationame << std::endl;
622  std::cout << correction_1[0]->GetXaxis()->GetNbins() << " bins, ";
623  std::cout << "first bin center " << correction_1[0]->GetBinCenter(0);
624  std::cout << ", last bin center " << correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()) << std::endl;
625  hratio[fiberbrightness] = new TH1D(hrationame, "",correction_1[0]->GetXaxis()->GetNbins(), correction_1[0]->GetBinCenter(0), correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()));
626  if(fiberbrightness==0){
627  std::cout << "FB=0" << std::endl;
628  hratio[0]->SetMinimum(0.95);
629  hratio[0]->SetMaximum(1.1);
630  }
631  std::cout << "made hratio[" << fiberbrightness << "], ";
632  if(manytomany) Ratiohist(correction_1[fiberbrightness], corrections_to_compare[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, labelx.Data(), "Correction factor ratio", 0.97, 1.03);//Old limits:1.8, 1.2, then 0.9, 1.2
633  else Ratiohist(correction_1[0], corrections_to_compare[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, labelx.Data(), "Correction factor ratio", 0.9, 1.2);
634  std::cout << "and filled it, ";
635  legend->AddEntry(hratio[fiberbrightness],analysis2+" "+fbStr+" to "+analysis1+" "+viewStr+" view");
636  std::cout << "and added it to legend." <<std::endl;
637  }
638  legend->Draw();
639  } else {
640  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data(), var.Data());
641  hratio[0] = new TH1D(hrationame, "",correction_1[0]->GetXaxis()->GetNbins(), correction_1[0]->GetBinCenter(0), correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()));
642  Ratiohist(correction_1, correction_2, hratio[0], kRed, 24, labelx.Data(), "Correction factor ratio", 1.8, 1.2);
643 
644  TLegend *legend = new TLegend(.5,.62,.88,.90);
645  legend->SetBorderSize(0);
646  legend->SetFillStyle(0);
647  legend->SetTextSize(0.04);
648  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
649  legend->AddEntry(hratio[0],analysis2+" to "+analysis1+" "+viewStr+" view");
650  legend->Draw();
651  }//end if(!fb)
652 
653  gPad->Print("Images/thresholdshadowinge_vs_"+var+"_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
654 
655  }//end of loop over views
656 }//end of thresholdshadowing()
657 //-------------------------------------------------------------------------------------------------------------------------------
658 void corrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell){
659  int nbrightnesses = 9;
660 
661  if(nddatamuoncatcher){
662  if(plane.Atoi() > 3 || plane.Atoi() < 2){
663  std::cout << "It is a ND muon catcher data file and has only 002, 003 planes. Plane " << plane << " is out of range. " << std::endl;
664  abort();
665  }
666  if(cell.Atoi() > 95 || cell.Atoi() < 0){
667  std::cout << "It is a Data file. Cell " << cell << " is out of range. " << std::endl;
668  abort();
669  }
670  }
671  if(det.Data()=="nd"){
672  if(mcdata.Data()=="mc"){
673  if(plane.Atoi() > 3 || plane.Atoi() < 0){
674  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
675  abort();
676  }
677  if(cell.Atoi() > 95 || cell.Atoi() < 0){
678  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
679  abort();
680  }
681  }//end of MC condition
682  if(mcdata.Data()=="data"){
683  if(plane.Atoi() > 191 || plane.Atoi() < 0){
684  std::cout << "It is a ND data file and has only 191 main body planes. Plane " << plane << " out of range. " << std::endl;
685  abort();
686  }
687  if(cell.Atoi() > 95 || cell.Atoi() < 0){
688  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
689  abort();
690  }
691  }//end of data condition
692  }//end of the ND condition
693 
694  if(det.Data()=="fd"){
695  if(mcdata.Data()=="mc"){
696  if(plane.Atoi() > 1 || plane.Atoi() < 0){
697  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
698  abort();
699  }
700  if(cell.Atoi() > 383 || cell.Atoi() < 0){
701  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
702  abort();
703  }
704  }//end of MC condition
705  else if(mcdata.Data()=="data"){
706  if(plane.Atoi() > 895 || plane.Atoi() < 0){
707  std::cout << "It is a Data file. Plane " << plane << " out of range. " << std::endl;
708  abort();
709  }
710  if(cell.Atoi() > 383 || cell.Atoi() < 0){
711  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
712  abort();
713  }
714  }//end of data condition
715  }//end of FD condition
716 
717  TFile* f1 = new TFile(c1, "read");
718  TFile* f2 = new TFile(c2, "read");
719 
720  TF1* fit1[9];
721  TF1* fit2[9];
722 
723  f1->cd();
724  std::cout << "getting " << c1 << ": plane_" << plane << "/Atten_Fit_" << plane << "_" << cell << std::endl;
725  TH1F* h1[9];
726  if(!manytomany){
727  h1[0] = (TH1F*)f1->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
728  h1[0]->ls();
729  Color_t color1 = kBlack;
730  if(mcdata=="mc"){
731  h1[0]->GetListOfFunctions()->Print();
732  if(plane.Atoi()==0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
733  else if(plane.Atoi()==1) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
734  else if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
735  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
736  if(fit1[0]) fit1[0]->SetLineColor(color1);
737  }//end of MC condition
738  else if(mcdata=="data"){
739  if(nddatamuoncatcher){
740  if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
741  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
742  if(fit1[0]) fit1[0]->SetLineColor(color1);
743  }//end of muon catcher condition
744  else if(plane.Atoi()%2==0) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
745  else if(plane.Atoi()%2!=0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
746  if(fit1[0]) fit1[0]->SetLineColor(color1);
747  }//end of data condition
748  if(!fit1[0]) std::cout << "Can't find fit1[0]" << std::endl;
749  } else {
750  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
751  //std::cout << "making h1[" << fiberbrightness << "]" << std::endl;
752  std::cout << "Filling h1[" << fiberbrightness << "]" << std::endl;
753 
754  TString fbStr = Form("fb%i_", fiberbrightness);
755  h1[fiberbrightness] = (TH1F*)f1->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell);
756 
757  std::cout << "Getting " << (fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell).Data() << std::endl;
758  h1[fiberbrightness]->GetListOfFunctions()->Print();
759  if(!h1[fiberbrightness]) std::cout << "Can't find h1 " << fiberbrightness << std::endl;
760  //h1[fiberbrightness]->SetMarkerStyle[24];
761 
762  //h2[fiberbrightness]->ls();
763  if(plane.Atoi()%2==0){
764  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitY");//X
765  std::cout << "using fitY" <<std::endl;
766  } else {
767  std::cout << "looking for fitX" << std::endl;
768  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");
769  if(fit1[fiberbrightness]) std::cout << "Got fitX" << std::endl;
770  }
771  fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
772  fit1[fiberbrightness]->SetLineStyle(kDotted);
773 
774  if(!fit1[fiberbrightness]) fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");//Y
775  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitMuX");
776 
777  if(fit1[fiberbrightness]){
778  fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
779  fit1[fiberbrightness]->SetLineStyle(kDotted);
780  } else {
781  std::cout << "can't find fit1[" << fiberbrightness << "]"<< std::endl;
782  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("rollFitY");
783  if(fit1[fiberbrightness]) fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
784  }
785  }
786  }
787 
788  f2->cd();
789  TH1F* h2[9];
790  //f2->Get("fb0_plane_000")->ls();
791  if(fb){
792  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
793  std::cout << "making h2[" << fiberbrightness << "]" << std::endl;
794  TString fbStr = Form("fb%i_", fiberbrightness);
795  h2[fiberbrightness] = (TH1F*)f2->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell);
796  //if(fiberbrightness==0)
797  //h2[fiberbrightness]->GetListOfFunctions()->Print();
798  if(!h2[fiberbrightness]) std::cout << "Can't find h2 " << fiberbrightness << std::endl;
799  //h2[fiberbrightness]->ls();
800  if(plane.Atoi()%2==0){
801  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitY");//X
802  std::cout << "using fitY" <<std::endl;
803  }
804  if(!fit2[fiberbrightness]) fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");//Y
805  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitMuX");
806  if(fit2[fiberbrightness]) fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
807  else {//std::cout << "can't find fit2[" << fiberbrightness << "]"<< std::endl;
808  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("rollFitY");
809  if(fit2[fiberbrightness]) fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
810  }
811  }
812  } else {
813  h2[0] = (TH1F*)f2->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
814  }
815  Color_t color1 = kBlack;
816  Color_t color2 = kRed;
817 
818  if(mcdata=="mc" && !fb){
819  if(plane.Atoi()==0) fit2[0] = (TF1*)h2[0]->FindObject("fitX");
820  else if(plane.Atoi()==1) fit2[0] = (TF1*)h2[0]->FindObject("fitY");
821  else if(plane.Atoi()==2) fit2[0] = (TF1*)h2[0]->FindObject("fitMuX");
822  else if(plane.Atoi()==3) fit2[0] = (TF1*)h2[0]->FindObject("fitMuY");
823  if(fit2[0]) fit2[0]->SetLineColor(color2);
824  }//end of MC condition
825  else if(mcdata=="data" && !fb){
826  if(nddatamuoncatcher){
827  if(plane.Atoi()==2) fit2[0] = (TF1*)h2[0]->FindObject("fitMuX");
828  else if(plane.Atoi()==3) fit2[0] = (TF1*)h2[0]->FindObject("fitMuY");
829  if(fit2[0]) fit2[0]->SetLineColor(color2);
830  }//end of muon catcher condition
831  else if(plane.Atoi()%2==0) fit2[0] = (TF1*)h2[0]->FindObject("fitY");
832  else if(plane.Atoi()%2!=0) fit2[0] = (TF1*)h2[0]->FindObject("fitX");
833  if(fit2[0]) fit2[0]->SetLineColor(color2);
834  }//end of data condition
835 
836  new TCanvas;
837  gPad->SetGrid();
838  h1[0]->Draw();
839  HistogramAttr1D(h1[0], "W (cm)", "Mean PE/cm", h1[0]->GetXaxis()->GetBinCenter(1), h1[0]->GetXaxis()->GetBinCenter(h1[0]->GetXaxis()->GetNbins()), 0., 80., color1);//old limits: 0, 150
840  if(!manytomany) h1[0]->SetLineColor(color1);
841  h1[0]->SetMarkerStyle(24);
842  h1[0]->SetMaximum(50);//old max: 120, then 60
843  if(fb) {
844  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
845  std::cout << "Drawing fiberbrightness = " << fiberbrightness << std::endl;
846  if(manytomany) {
847  h1[fiberbrightness]->Draw("same");
848  h1[fiberbrightness]->SetLineStyle(kDotted);
849  h1[fiberbrightness]->SetMarkerColor(fbcolors[fiberbrightness]);
850  h1[fiberbrightness]->SetMarkerStyle(25);
851  h1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
852  }
853  h2[fiberbrightness]->Draw("same");
854  h2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
855  h2[fiberbrightness]->SetMarkerColor(fbcolors[fiberbrightness]);
856  h2[fiberbrightness]->SetMarkerStyle(24);
857  }
858  } else {
859  h2[0]->Draw("same");
860  h2[0]->SetLineColor(color2);
861  h2[0]->SetMarkerColor(color2);
862  h2[0]->SetMarkerStyle(24);
863  }
864 
865  TLegend *legend = new TLegend(.12,.4,.64,.90);
866  legend->SetBorderSize(0);
867  legend->SetFillStyle(0);
868  //legend->SetTextSize(0.04);
869  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
870  if(!manytomany)legend->AddEntry(h1[0], analysis1);
871  else {
872  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
873  TString fbStr = Form(", FB %i", fiberbrightness);
874  legend->AddEntry(h1[fiberbrightness], analysis1+fbStr);
875  }
876  }
877  if(fb) {
878  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
879  TString fbStr = Form(", FB %i", fiberbrightness);
880  legend->AddEntry(h2[fiberbrightness], analysis2+fbStr);
881  }
882  } else {
883  legend->AddEntry(h2[0], analysis2);
884  }
885  legend->Draw();
886  gPad->Print("Images/Atten_Fit_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+plane+"_"+cell+".pdf");
887 
888  if(manytomany) {
889  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
890  new TCanvas();
891  h1[fiberbrightness]->Draw();
892  h2[fiberbrightness]->Draw("same");
893  TLegend * oneFBlegend = new TLegend(0.5, 0.7, 0.9, 0.9);
894  const TString viewStr = (plane.Atoi()%2==0) ? "Y" : "X";
895  oneFBlegend->AddEntry(h1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
896  oneFBlegend->AddEntry(h2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
897  oneFBlegend->Draw();
898  TString fbstr;
899  fbstr.Form("fb%i", fiberbrightness);
900  gPad->Print("Images/Atten_Fit_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+plane+"_"+cell+"_"+fbstr+".pdf");
901  }
902  }
903 
904  //ratio of data points
905  new TCanvas;
906  gPad->SetGrid();
907  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
908  TH1D* hratio[9];
909  TString hrationamearr[9];
910  if(fb) {
911  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
912  hrationamearr[fiberbrightness] = Form("%s_fb%i", hrationame.Data(), fiberbrightness);
913  //hratio[fiberbrightness] = new TH1D(hrationamearr[fiberbrightness], "",h1->GetXaxis()->GetNbins(), h1->GetBinCenter(0), h1->GetBinCenter(h1->GetXaxis()->GetNbins()));
914  hratio[fiberbrightness] = new TH1D(hrationamearr[fiberbrightness], "",h1[0]->GetXaxis()->GetNbins(), -800, 800);
915  if(manytomany) Ratiohist(h1[fiberbrightness], h2[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 2.0);//old limits: 0.3, 1
916  else Ratiohist(h1[0], h2[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 2.0);//old limits: 0.3, 1
917  }
918  } else {
919  hratio[0] = = new TH1D(hrationame, "",h1[0]->GetXaxis()->GetNbins(), h1[0]->GetBinCenter(0), h1[0]->GetBinCenter(h1[0]->GetXaxis()->GetNbins()));
920  Ratiohist(h1[0], h2[0], hratio[0], kRed, 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 1.0);
921  }
922 
923  TLegend *legend = new TLegend(.12,.62,.64,.90);
924  legend->SetBorderSize(0);
925  legend->SetFillStyle(0);
926  legend->SetTextSize(0.04);
927  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
928  if(fb) {
929  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
930  TString fbStr = Form(" FB %i", fiberbrightness);
931  legend->AddEntry(hratio[fiberbrightness],analysis2+fbStr+" to "+analysis1);
932  }
933  } else {
934  legend->AddEntry(hratio[0],analysis2+" to "+analysis1);
935  }
936  legend->Draw();
937  gPad->Print("Images/Atten_Fit_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");
938 
939  //ratio of fitting curve
940  new TCanvas;
941  gPad->SetGrid();
942  const TString hratiofitname = Form("hratiofit%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
943  TString hratiofitnamearr[9];
944  float x_limit;
945  int x_bins;
946  if(det=="nd"){
947  x_limit = 250.;
948  x_bins = 2*x_limit;
949  }
950  else if(det=="fd"){
951  x_limit = 800.;
952  x_bins = 2*x_limit;
953  }
954 
955  TH1D* hfit = new TH1D("hfit", "", x_bins, -x_limit, +x_limit);
956  std::cout << "l592 about to use HistogramAttr1D" << std::endl;
957  HistogramAttr1D(hfit, "W (cm)", "Ratio of Mean PE/cm fitting curves", -x_limit, +x_limit, 0.5, 1.5, kWhite);//old limits: 0.55, 0.68
958  hfit->Draw();
959  TH1D* hratiofit[9];
960  if(fb){
961  std::cout << "running over " << nbrightnesses << "brightnesses" << std::endl;
962  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
963  hratiofitnamearr[fiberbrightness] = Form("%s_fb%i", hratiofitname.Data(), fiberbrightness);
964  hratiofit[fiberbrightness] = new TH1D(hratiofitnamearr[fiberbrightness], "", x_bins, -x_limit, +x_limit);
965  Ratiofit(fit1[0], fit2[fiberbrightness], hratiofit[fiberbrightness], plane.Data(), cell.Data(), fbcolors[fiberbrightness], 24, "W (cm)", "Ratio of Mean PE/cm fitting curves", 0.7, 2.0);//old limits: 0.55, 0.68
966  }
967  } else {
968  hratiofit[0] = new TH1D(hratiofitname, "", x_bins, -x_limit, +x_limit);
969  Ratiofit(fit1[0], fit2[0], hratiofit[0], plane.Data(), cell.Data(), color2, 24, "W (cm)", "Ratio of Mean PE/cm fitting curves", 0.55, 0.68);
970  }
971 
972  TLegend *legend = new TLegend(.12,.62,.64,.90);
973  legend->SetBorderSize(0);
974  legend->SetFillStyle(0);
975  legend->SetTextSize(0.04);
976  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
977 
978  if(fb) {
979  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
980  TString fbStr = Form(" FB %i", fiberbrightness);
981  legend->AddEntry(hratiofit[fiberbrightness],analysis2+fbStr+" to "+analysis1);
982  }
983  } else {
984  legend->AddEntry(hratiofit,analysis2+" to "+analysis1);
985  }
986  legend->Draw();
987  gPad->Print("Images/Atten_Fit_fitratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");
988 }//end of corrected_PE()
989 
990 //---------------------------------------------------------------------------------------------------
991 void calibrationconstants(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2){
992  ifstream in_1;
993  ifstream in_2;
994 
995  ofstream out;
996 
997  FILE* fout = fopen("uncalibrated_channels_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+".csv","w");
998  fprintf(fout, "%s%50s%s\n", analysis1.Data(), " ", analysis2.Data());
999  fprintf(fout, "channel,coeff_exp,atten_length,background,chisq, %9s, channel,coeff_exp,atten_length,background,chisq\n", " ");
1000 
1001  in_1.open(c1);
1002  in_2.open(c2);
1003 
1004 
1005  //x=coeff_exp, y=atten_length, z=background, a=edgx_low, b=edgx_high, c=coeff_low
1006  Float_t x_1, y_1, z_1, a_1, b_1, c_1, d_1, e_1, g_1;
1007  Float_t x_2, y_2, z_2, a_2, b_2, c_2, d_2, e_2, g_2;
1008  Float_t chisq = 0.2;
1009 
1010  Int_t channels;
1011  Int_t totalchannels1 = 0;
1012  Int_t totalchannels2 = 0;
1013  Int_t calibratedchannels1 = 0;
1014  Int_t calibratedchannels2 = 0;
1015  Int_t f_1;
1016  Int_t f_2;
1017  Int_t bins1;
1018 
1019  Double_t limit1D = 3.;
1020  Double_t limit2D = 1000.;
1021 
1022  if(det=="nd" && mcdata=="mc" ) { channels = 3063; }
1023  if(det=="nd" && mcdata=="data" ) { channels = 1003063; }
1024  if(det=="fd" && mcdata=="mc" ) { channels = 1383; }
1025  if(det=="fd" && mcdata=="data" ) { channels = 895383; }
1026 
1027 
1028  TH1F* hcoeffexp = new TH1F("hcoeffexp", "coeffexp difference distribution", 120, -limit1D, +limit1D );
1029  TH1F* hattenlength = new TH1F("hattenlength", "attenlength difference distribution", 600, -limit1D, +limit1D );
1030  TH1F* hbackground = new TH1F("hbackground", "background difference distribution", 120, -limit1D, +limit1D );
1031  TH1F* hchisq = new TH1F("hchisq", "chisq difference distribution", 600, -limit1D, +limit1D );
1032 
1033  TH2F* h2coeffexp = new TH2F("h2coeffexp", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1034  TH2F* h2attenlength = new TH2F("h2attenlength"," " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1035  TH2F* h2background = new TH2F("h2background", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1036  TH2F* h2chisq = new TH2F("h2chisq"," ", 200,0.,2., 200,0.,2. );
1037 
1038  while(1){
1039  in_1 >> x_1 >> y_1 >> z_1 >> a_1 >> b_1 >> c_1 >> d_1 >> e_1 >> f_1 >> g_1;
1040  in_2 >> x_2 >> y_2 >> z_2 >> a_2 >> b_2 >> c_2 >> d_2 >> e_2 >> f_2 >> g_2;
1041  if(!in_1.good()) break;
1042  if(!in_2.good()) break;
1043 
1044  if(f_1<=channels && e_1 < chisq && e_2 < chisq){
1045  hcoeffexp->Fill((x_1 - x_2)/x_1);
1046  hattenlength->Fill((y_1 - y_2)/y_1);
1047  hbackground->Fill((z_1 - z_2)/z_1);
1048 
1049  h2coeffexp->Fill(x_2, x_1);
1050  h2attenlength->Fill(y_2, y_1);
1051  h2background->Fill(z_2, z_1);
1052  }//end of loop over channels and chisq conditions.
1053  if(f_1<=channels){
1054  totalchannels1++;
1055  totalchannels2++;
1056  if(e_1<chisq) calibratedchannels1++;
1057  if(e_2<chisq) calibratedchannels2++;
1058 
1059  hchisq->Fill((e_1 - e_2)/e_1);
1060  h2chisq->Fill(e_2, e_1);
1061 
1062  if(e_1>=chisq || e_2>=chisq){
1063  fprintf(fout, "%i,%g,%g,%g,%g,%20s,%i,%g,%g,%g,%g\n",f_1,x_1,y_1,z_1,e_1, " ", f_2,x_2,y_2,z_2,e_2);
1064  out << f_1 << setw(22) << x_1 << setw(22) << y_1 << setw(22) << z_1 << setw(22) << e_1 << f_2 << setw(22) << x_2 << setw(22) << y_2 << setw(22) << z_2 << setw(22) << e_2 << endl;
1065  }
1066  }//end of only channels condition for chisq plots. For chisq plots we have to consider all channels.
1067  }//end of while condition.
1068 
1069  cout << " " << endl;
1070  cout << "Total channels in "+analysis1+" = " << totalchannels1 << " and calibrated channels = " << calibratedchannels1 << " = " << ((float)calibratedchannels1/(float)totalchannels1)*100. << " %." << endl;
1071  cout << "Total channels in "+analysis2+" = " << totalchannels2 << " and calibrated channels = " << calibratedchannels2 << " = " << ((float)calibratedchannels2/(float)totalchannels2)*100. << " %." << endl;
1072  cout << " " << endl;
1073 
1074  //coeffexp plots
1075  new TCanvas;
1076  gPad->SetGrid();
1077  std::cout << "l712 about to use HistogramAttr1D" << std::endl;
1078  HistogramAttr1D(hcoeffexp, "#Deltacoeffexp/coeffexp", "# of channels", 0.5*hcoeffexp->GetXaxis()->GetBinCenter(1), 0.5*hcoeffexp->GetXaxis()->GetBinCenter(hcoeffexp->GetXaxis()->GetNbins()), 0., 7e5, kBlack);
1079  hcoeffexp->Draw();
1080  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1081  legend->SetBorderSize(0);
1082  legend->SetFillStyle(0);
1083  legend->SetTextSize(0.05);
1084  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1085  legend->AddEntry(hcoeffexp,"#splitline{"+analysis1.Data()+"}{vs "+analysis2.Data()+"}","P");
1086  legend->Draw();
1087  gPad->Print("Images/Delta_coeffexp_"+analysis1.Data()+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1088 
1089  new TCanvas;
1090  gPad->SetGrid();
1091  gPad->SetLogz();
1092 
1093  HistogramAttr2D(h2coeffexp, "coeffexp_{"+analysis2+"}", "coeffexp_{"+analysis1+"}", " ", 0, 100, 0, 100);
1094  h2coeffexp->Draw("colz");
1095  gPad->Print("Images/coeffexp2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1096 
1097  //attenlength plots
1098  new TCanvas;
1099  gPad->SetGrid();
1100  HistogramAttr1D(hattenlength, "#Deltaattenlength/attenlength", "# of channels", 0.5*hattenlength->GetXaxis()->GetBinCenter(1), 0.5*hattenlength->GetXaxis()->GetBinCenter(hattenlength->GetXaxis()->GetNbins()), 0., 166500, kBlack);
1101  hattenlength->Draw();
1102  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1103  legend->SetBorderSize(0);
1104  legend->SetFillStyle(0);
1105  legend->SetTextSize(0.05);
1106  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1107  legend->AddEntry(hattenlength,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1108  legend->Draw();
1109  gPad->Print("Images/Delta_attenlength_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1110 
1111  new TCanvas;
1112  gPad->SetGrid();
1113  gPad->SetLogz();
1114  HistogramAttr2D(h2attenlength, "attenlength_{"+analysis2+"} (cm)", "attenlength_{"+analysis1+"} (cm)", " ", 0, 1000, 0, 1000);
1115  h2attenlength->Draw("colz");
1116  gPad->Print("Images/attenlength2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1117 
1118  //backgrounds plots
1119  new TCanvas;
1120  gPad->SetGrid();
1121  HistogramAttr1D(hbackground, "#Deltabackground/background", "# of channels", 0.5*hbackground->GetXaxis()->GetBinCenter(1), 0.5*hbackground->GetXaxis()->GetBinCenter(hbackground->GetXaxis()->GetNbins()), 0., 3e4, kBlack);
1122  hbackground->Draw();
1123  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1124  legend->SetBorderSize(0);
1125  legend->SetFillStyle(0);
1126  legend->SetTextSize(0.05);
1127  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1128  legend->AddEntry(hbackground,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1129  legend->Draw();
1130  gPad->Print("Images/Delta_background_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1131 
1132  new TCanvas;
1133  gPad->SetGrid();
1134  gPad->SetLogz();
1135  HistogramAttr2D(h2background, "background_{"+analysis2+"}", "background_{"+analysis1+"}", " ", -50, 50, -50, 50);
1136  h2background->Draw("colz");
1137  gPad->Print("Images/background2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1138 
1139  //chisq plots
1140  new TCanvas;
1141  gPad->SetGrid();
1142  HistogramAttr1D(hchisq, "#Deltachisq/chisq", "# of channels", 0.5*hchisq->GetXaxis()->GetBinCenter(1), 0.5*hchisq->GetXaxis()->GetBinCenter(hchisq->GetXaxis()->GetNbins()), 0., 4200., kBlack);
1143  hchisq->Draw();
1144  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1145  legend->SetBorderSize(0);
1146  legend->SetFillStyle(0);
1147  legend->SetTextSize(0.05);
1148  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1149  legend->AddEntry(hchisq,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1150  legend->Draw();
1151  gPad->Print("Images/Delta_chisq_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1152 
1153  new TCanvas;
1154  gPad->SetGrid();
1155  gPad->SetLogz();
1156  TString chisq_tstring = "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1.Data+"}";
1157  char * chisq_char = chisq_tstring.Data();
1158  HistogramAttr2D(h2chisq, "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1+"}", " ", 0, 0.5, 0, 0.5);
1159  h2chisq->Draw("colz");
1160  gPad->Print("Images/chisq2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1161 
1162  in_1.close();
1163  in_2.close();
1164  out.close();
1165  fclose(fout);
1166 }//end of calibrationconstants()
1167 //---------------------------------------------------------------------------------------------------
1168 
1169 void zerosstudy(TFile * f0, TFile * f1, TFile * f2, TFile * f3, TFile * f4, TFile * f5, TFile * f6, TFile * f7, TFile * f8)
1170 {
1171  const int nbrightnesses = 9;
1172  //int nviews = 2;
1173  TFile * file[nbrightnesses] = {f0, f1, f2, f3, f4, f5, f6, f7, f8};
1174  TLegend * lzeros = new TLegend(0.5, 0.1, 0.88, 0.5);
1175  lzeros->SetBorderSize(0);
1176  lzeros->SetFillStyle(0);
1177  lzeros->SetTextSize(0.05);
1178 
1179  TH1F * hZerosVsWX[nbrightnesses];
1180  TH1F * hZerosVsWY[nbrightnesses];
1181  TH1F * hNCalibVsWX[nbrightnesses];
1182  TH1F * hNCalibVsWY[nbrightnesses];
1183  TH1F * hZerosVsCellX[nbrightnesses];
1184  TH1F * hZerosVsCellY[nbrightnesses];
1185  TH1F * hNCalibVsCellX[nbrightnesses];
1186  TH1F * hNCalibVsCellY[nbrightnesses];
1187 
1188  TCanvas * cZerosVsWX = new TCanvas;
1189  TCanvas * cZerosVsWY = new TCanvas;
1190  TCanvas * cCalibVsWX = new TCanvas;
1191  TCanvas * cCalibVsWY = new TCanvas;
1192  TCanvas * cZerosVsCellX = new TCanvas;
1193  TCanvas * cZerosVsCellY = new TCanvas;
1194  TCanvas * cCalibVsCellX = new TCanvas;
1195  TCanvas * cCalibVsCellY = new TCanvas;
1196 
1197  for(int fb = 0; fb < nbrightnesses; ++fb){
1198  TString fbstr;
1199  fbstr.Form("fb%i", fb);
1200  TString FBStr;
1201  FBStr.Form("FB %i", fb);
1202 
1203  std::cout << "Looking for nZerosVsW_xview_" << fbstr << std::endl;
1204  hZerosVsWX[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsW_xview_"+fbstr);
1205  hZerosVsWY[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsW_yview_"+fbstr);
1206  //hCalibVsWX[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsW_xview_"+fbstr);
1207  //hCalibVsWY[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsW_yview_"+fbstr);
1208  hZerosVsCellX[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsCell_xview_"+fbstr);
1209  hZerosVsCellY[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsCell_yview_"+fbstr);
1210  //hCalibVsCellX[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsCell_xview_"+fbstr);
1211  //hCalibVsCellY[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsCell_yview_"+fbstr);
1212 
1213  hZerosVsWX[fb]->SetTitle("XView");
1214  hZerosVsWY[fb]->SetTitle("YView");
1215 
1216  hZerosVsWX[fb]->SetLineColor(fbcolors[fb]);
1217  hZerosVsWY[fb]->SetLineColor(fbcolors[fb]);
1218  //hCalibVsWX[fb]->SetLineColor(fbcolors[fb]);
1219  //hCalibVsWY[fb]->SetLineColor(fbcolors[fb]);
1220  hZerosVsCellX[fb]->SetLineColor(fbcolors[fb]);
1221  hZerosVsCellY[fb]->SetLineColor(fbcolors[fb]);
1222  //hCalibVsCellX[fb]->SetLineColor(fbcolors[fb]);
1223  //hCalibVsCellY[fb]->SetLineColor(fbcolors[fb]);
1224 
1225  lzeros->AddEntry(hZerosVsWX[fb], FBStr.Data(), "l");
1226 
1227  cZerosVsWX->cd();
1228  hZerosVsWX[fb]->Draw("same");
1229  cZerosVsWY->cd();
1230  hZerosVsWY[fb]->Draw("same");
1231  cCalibVsWX->cd();
1232  //hCalibVsWX[fb]->Draw("same");
1233  cCalibVsWY->cd();
1234  //hCalibVsWY[fb]->Draw("same");
1235  cZerosVsCellX->cd();
1236  hZerosVsCellX[fb]->Draw("same");
1237  cZerosVsCellY->cd();
1238  hZerosVsCellY[fb]->Draw("same");
1239  cCalibVsCellX->cd();
1240  //hCalibVsCellX[fb]->Draw("same");
1241  cCalibVsCellY->cd();
1242  //hCalibVsCellY[fb]->Draw("same");
1243 
1244  }
1245 
1246  cZerosVsWX->cd();
1247  lzeros->Draw();
1248  cZerosVsWX->SaveAs("zeros_W_X.pdf");
1249 
1250  cZerosVsWY->cd();
1251  lzeros->Draw();
1252  cZerosVsWY->SaveAs("zeros_W_Y.pdf");
1253 
1254  cCalibVsWX->cd();
1255  lzeros->Draw();
1256  cCalibVsWX->SaveAs("calib_W_X.pdf");
1257 
1258  cCalibVsWY->cd();
1259  lzeros->Draw();
1260  cCalibVsWY->SaveAs("calib_W_Y.pdf");
1261 
1262  cZerosVsWX->cd();
1263  lzeros->Draw();
1264  cZerosVsWX->SaveAs("zeros_W_X.pdf");
1265 
1266  cZerosVsWY->cd();
1267  lzeros->Draw();
1268  cZerosVsWY->SaveAs("zeros_W_Y.pdf");
1269 
1270  cCalibVsWX->cd();
1271  lzeros->Draw();
1272  cCalibVsWX->SaveAs("calib_W_X.pdf");
1273 
1274  cCalibVsWY->cd();
1275  lzeros->Draw();
1276  cCalibVsWY->SaveAs("calib_W_Y.pdf");
1277 
1278 }
1279 
1280 // LocalWords: ProfileX
f4
Float_t f4
Definition: comparison_ascii.C:56
kRed
enum BeamMode kRed
PandAna.Demos.tute_pid_validation.var
list var
Definition: tute_pid_validation.py:17
Ratiohist
void Ratiohist(TH1 *h1, TH1 *h2, TH1 *hratio, Color_t color, Style_t mstyle, TString labelx, TString labely, Double_t lowy, Double_t highy)
Definition: HistogramAttr.cxx:45
hist2
TH1D * hist2
Definition: plotHisto.C:12
zerosplots
bool zerosplots
Definition: attenuation_calibration_comparisons_withfb.C:31
mc
bool mc
Definition: attenuation_calibration_comparisons_withfb.C:37
calibrationconstants
void calibrationconstants(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2)
Definition: attenuation_calibration_comparisons_withfb.C:991
submit_syst.fout
fout
Definition: submit_syst.py:244
onetoone
bool onetoone
Definition: attenuation_calibration_comparisons_withfb.C:40
manytomany
bool manytomany
Definition: attenuation_calibration_comparisons_withfb.C:42
Get
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
f2
Float_t f2
Definition: comparison_ascii.C:56
Munits::c_2
static constexpr Double_t c_2
Definition: Munits.h:245
make_mec_shifts_plots.legend
legend
Definition: make_mec_shifts_plots.py:148
thresholdshadowingplots
bool thresholdshadowingplots
Definition: attenuation_calibration_comparisons_withfb.C:28
data
const XML_Char const XML_Char * data
Definition: expat.h:268
demo5.c2
c2
Definition: demo5.py:33
manytoone
bool manytoone
Definition: attenuation_calibration_comparisons_withfb.C:41
fclose
fclose(fg1)
tricellvstrajectory
bool tricellvstrajectory
Definition: attenuation_calibration_comparisons_withfb.C:39
f3
Float_t f3
Definition: comparison_ascii.C:56
calibrationconstantsplots
bool calibrationconstantsplots
Definition: attenuation_calibration_comparisons_withfb.C:30
channels
std::map< ToFCounter, std::vector< unsigned int > > channels
Definition: ToFPositionRecoAnalysis_module.cc:81
HistogramAttr2D
void HistogramAttr2D(TH2 *h, char *titlx_x, char *titlx_y, char *titlx_z, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy)
Definition: HistogramAttr.cxx:24
check_time_usage.float
float
Definition: check_time_usage.py:7
uncorrectedPEplots
bool uncorrectedPEplots
Definition: attenuation_calibration_comparisons_withfb.C:27
GetXaxis
correl_xv GetXaxis() -> SetDecimals()
Definition: plot_lowE_sig_spread_genie.C:84
hist1
TH1D * hist1
Definition: plotHisto.C:9
fb
bool fb
Definition: attenuation_calibration_comparisons_withfb.C:45
f1
Float_t f1
Definition: comparison_ascii.C:56
attenuation_calibration_comparisons_withfb
attenuation_calibration_comparisons_withfb()
Definition: attenuation_calibration_comparisons_withfb.C:50
fbcolors
int fbcolors[9]
Definition: attenuation_calibration_comparisons_withfb.C:47
thresholdshadowing
void thresholdshadowing(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var)
Definition: attenuation_calibration_comparisons_withfb.C:481
om::cout
OStream cout
Definition: OStream.cxx:6
h2
TH1F * h2
Definition: plot.C:45
Ratiofit
void Ratiofit(TF1 *f1, TF1 *f2, TH1 *hratio, TString plane, TString cell, Color_t color, Style_t mstyle, TString labelx, TString labely, Double_t lowy, Double_t highy)
Definition: HistogramAttr.cxx:65
correctedPEplots
bool correctedPEplots
Definition: attenuation_calibration_comparisons_withfb.C:29
h1
TH1F * h1
Definition: nslices_potspill.C:4
fd
bool fd
Definition: attenuation_calibration_comparisons_withfb.C:36
nddatamuoncatcher
bool nddatamuoncatcher
Definition: attenuation_calibration_comparisons_withfb.C:43
uncorrected_PE
uncorrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
Definition: attenuation_calibration_comparisons_withfb.C:313
file
TFile * file
Definition: cellShifts.C:17
demo5.c1
c1
Definition: demo5.py:24
HistogramAttr1D
HistogramAttr1D(TH1 *h, char *title_x, char *title_y, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy, Color_t lcolor)
Definition: HistogramAttr.cxx:5
uncorrected_PE_tricellvstrajectory
void uncorrected_PE_tricellvstrajectory(char *c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
Definition: attenuation_calibration_comparisons_withfb.C:477
corrected_PE
void corrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell)
Definition: attenuation_calibration_comparisons_withfb.C:658
zerosstudy
void zerosstudy(TFile *f0, TFile *f1, TFile *f2, TFile *f3, TFile *f4, TFile *f5, TFile *f6, TFile *f7, TFile *f8)
Definition: attenuation_calibration_comparisons_withfb.C:1169
nd
bool nd
Definition: attenuation_calibration_comparisons_withfb.C:35
Generated on Fri Feb 26 2021 14:13:25 for NOvA by   doxygen 1.8.11