attenuation_calibration_comparisons_withfb_temp.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 = false; // 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 
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* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain100/uncorrected_attenprofs/calibhists_updatedlightlevels_fdmc_gain100.root";
201  //char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain100/thresh/thresh_fits_thirdanalysis_updatedlightlevels_fd_gain100_fb0-8.root";
202  //char* cattenuation1 = "/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* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/uncorrected_attenprofs/calibhists_updatedlightlevels_fdmc_gain140.root";
206  char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/thresh/thresh_fits_thirdanalysis_updatedlightlevels_fd_gain140_fb0-8.root";
207  //char* cattenuation1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/updatedlightlevels/fdmc_gain140/corrected_attenprofs/Atten_Fit_Results.root";
208 
209  //updated light levels, real conditions fdmc low gain LOW STATS
210  //char* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain100/uncorrected_attenprofs/calibhists_realandupdatedlightlevels_fdmc_gain100.root";
211  //char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain100/thresh/thresh_fits_thirdanalysis_realandupdatedlightlevels_fd_gain100_fb0-8.root";
212  //char* cattenuation1 = "/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 LOW STATS
215  //char* c1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain140/uncorrected_attenprofs/calibhists_realandupdatedlightlevels_fdmc_gain140.root";
216  //char* cthreshold1 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/real_and_updatedlightlevels/fdmc_gain140/thresh/thresh_fits_thirdanalysis_realandupdatedlightlevels_fd_gain140_fb0-8.root";
217  //char* cattenuation1 = "/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  //updated light levels, real conditions fdmc low gain FULL STATS
225  //char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/period2ONLY__USE_THIS/uncorrected_attenprofs/calib_hist_fdmc_lowgain_highstatsrealconditions.root";
226  //char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/period2ONLY__USE_THIS/threshold/thresh_fits_thirdanalysis_fd_gain100_histats_fb0-8.root";
227  //char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/period2ONLY__USE_THIS/corrected_attenprofs/Atten_Fit_Results.root";
228 
229  //updated light levels, real conditions fdmc high gain FULL STATS
230  char* c2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/highgain_USE_THIS/uncorrected_attenprofs/calib_hists_fdmc_highgain_realconditions_fullstats.root";
231  char* cthreshold2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/highgain_USE_THIS/threshold/thresh_fits_thirdanalysis_fd_gain140_histats_fb0-8.root";
232  char* cattenuation2 = "/nova/ana/users/mcampbel/calibration_thirdanalysis/july6_2017_fdmc/highgain_USE_THIS/corrected_attenprofs/Atten_Fit_Results.root";
233 
234  ///nova/ana/users/mcampbel/calibration_thirdanalysis/trajectoryhits/fdmc_gain140/
235  //----------------------------------------------------------------------------------------------------
236 
237  //.root file for uncorrected PE/cm for tricell vs trajectory study
238  char* c3 = "";
239 
240  //---------------------------------------------------------------------------------------------------
241 
242  TString hist1; // base histogram from the base file. You want to compare other histograms wrt base histogram.
243  TString hist2; // histogram from the second file that you want to compare with the base histogram.
244  TString det; // name of the detector.
245  TString mcdata; // labeling for the Monte Carlo or the data .
246  TString analysis1; // labeling for the base file.
247  TString analysis2; // labelling for the second file.
248  TString var; // variable labelling on the x axis. It is W (cm).
249  TString plane; // plane number.
250  TString cell; // cell number.
251 
252  //analysis1 = "2nd_analysis_period2";
253  //analysis1 = "2nd_analysis";
254  analysis1 = "IdealConditions";
255  analysis2 = "FullStatsRealConditions";
256  //analysis1 = "LowStats";
257  //analysis2 = "FullStats";
258 
259  hist1 = "h_WPE_corr_xy";
260  hist2 = "h_WPE_corr_xy"; //_corr_xy_ for the tricell, _corr_traj_ for the trajectory
261 
262  var = "W";//W or Cell
263  plane = "000";
264  cell = "142";
265 
266  //ND MC two cases comparison
267  if(nd && onetoone && !fb){
268  det = "nd";
269  if(mc) mcdata = "mc";
270  if(data) mcdata = "data";
271 
272  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
273  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
274 
275  if(mc && thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
276 
277  if(data && nddatamuoncatcher && correctedPEplots)corrected_PE(cattenuation_muoncatcher1, cattenuation_muoncatcher2, det, mcdata, analysis1, analysis2, plane, cell);
278  else if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
279 
280  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
281  }//end of if(nd && mc && onetoone)
282 
283  //FD two case comparison
284  if(fd && onetoone && !fb){
285  det = "fd";
286  if(mc) mcdata = "mc";
287  if(data) mcdata = "data";
288 
289  if(tricellvstrajectory && uncorrectedPEplots) uncorrected_PE_tricellvstrajectory(c3, det, mcdata, analysis1, analysis2, hist1, hist2);
290  else if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
291 
292  if(mc && thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
293 
294  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
295 
296  if(calibrationconstantsplots) calibrationconstants(cconstants1, cconstants2, det, mcdata, analysis1, analysis2);
297  }//end of if(fd && mc && onetoone)
298 
299  if(fb && fd){
300  det = "fd";
301  mcdata = "mc";
302 
303  if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
304  if(thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
305  if(correctedPEplots){
306  corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
307  // corrected_PE_mean(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
308  }
309  }
310 
311  if(fb && nd){
312  det = "nd";
313  mcdata = "mc";
314 
315  if(uncorrectedPEplots) uncorrected_PE(c1, c2, det, mcdata, analysis1, analysis2, hist1, hist2);
316  if(thresholdshadowingplots) thresholdshadowing(cthreshold1, cthreshold2, det, mcdata, analysis1, analysis2, var);
317  if(correctedPEplots) corrected_PE(cattenuation1, cattenuation2, det, mcdata, analysis1, analysis2, plane, cell);
318  }
319 
320  //if(zerosplots) zerosstudy(zerosfile0, zerosfile1, zerosfile2, zerosfile3, zerosfile4, zerosfile5, zerosfile6, zerosfile7, zerosfile8);
321 
322 
323 }//attenuation_calibration_comparisons()
324 /////////////////////////////////////////////////////////////////////////////////
325 
326 //----------------------------------------------------------------------------------------------------------------------
327 uncorrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
328  TFile* f1 = new TFile(c1, "read");
329  TFile* f2 = new TFile(c2, "read");
330  TString viewstr;
331  TString viewStr;
332  for(int view = 0; view < 2; ++view){
333  if(nddatamuoncatcher) viewstr = view ? "muy" : "mux";
334  else if(!nddatamuoncatcher) viewstr = view ? "y" : "x";
335  if(nddatamuoncatcher) viewStr = view ? "muY" : "muX";
336  else if(!nddatamuoncatcher) viewStr = view ? "Y" : "X";
337 
338  TH2F* h2D_1[9];
339  TH2F* h2D_2[9];
340  TProfile* ProfileX_1 [9];
341  TProfile* ProfileX_2 [9];
342 
343  f1->cd();
344  if(!manytomany){
345  if(f1->GetListOfKeys()->Contains("calib"))
346  h2D_1[0] = (TH2F*)f1->Get("calib/"+hist1+"_in"+viewStr);
347  if(f1->GetListOfKeys()->Contains("make")){
348  h2D_1[0] = (TH2F*)f1->Get("make/"+hist1+"_in"+viewStr);
349  std::cout << "h1 = make/" << hist1 << "_in" << viewStr << std::endl;
350  }
351  ProfileX_1[0] = h2D_1[0]->ProfileX();
352  if(det=="nd" && mcdata=="data") ProfileX_1[0]->Rebin();
353  } else {
354  TString histdir;
355  if(f1->GetListOfKeys()->Contains("calib")) histdir = "calib/";
356  else if(f1->GetListOfKeys()->Contains("make")) histdir = "make/";
357  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
358  TString fbStr = Form("_fb%i", fiberbrightness);
359  h2D_1[fiberbrightness] = (TH2F*)f1->Get(histdir+hist1+"_in"+viewStr+fbStr);
360  ProfileX_1[fiberbrightness] = h2D_1[fiberbrightness]->ProfileX();
361  }
362  }
363 
364  f2->cd();
365  if(fb){
366  TString histdir;
367  if(f2->GetListOfKeys()->Contains("calib")) histdir = "calib/";
368  else if(f2->GetListOfKeys()->Contains("make")) histdir = "make/";
369  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
370  TString fbStr = Form("_fb%i", fiberbrightness);
371  h2D_2[fiberbrightness] = (TH2F*)f2->Get(histdir+hist2+"_in"+viewStr+fbStr);
372  ProfileX_2[fiberbrightness] = h2D_2[fiberbrightness]->ProfileX();
373  }
374  } else {
375  if(f2->GetListOfKeys()->Contains("calib"))
376  h2D_2[0] = (TH2F*)f2->Get("calib/"+hist2+"_in"+viewStr);
377  else if(f2->GetListOfKeys()->Contains("make"))
378  h2D_2[0] = (TH2F*)f2->Get("make/"+hist2+"_in"+viewStr);
379  ProfileX_2[0] = h2D_2[0]->ProfileX();
380  if(det=="nd" && mcdata=="data") ProfileX_2[0]->Rebin();
381  }
382 
383  new TCanvas;
384  gPad->SetGrid();
385 
386  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
387  ProfileX_1[0]->Draw();
388  ProfileX_1[0]->SetLineColor(kBlack);
389 
390  if(fb){
391  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
392  ProfileX_2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
393  ProfileX_2[fiberbrightness]->Draw("same");
394  if(manytomany){
395  ProfileX_1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
396  ProfileX_1[fiberbrightness]->SetLineStyle(kDotted);
397  if(fiberbrightness > 0) ProfileX_1[fiberbrightness]->Draw("same");
398  }
399  }
400  } else {
401  ProfileX_2[0]->SetLineColor(kRed);
402  ProfileX_2[0]->Draw("same");
403  }
404 
405  TLegend *legend = new TLegend(.12,.5,.88,.90);//was .12, .62, .64, .90
406  legend->SetBorderSize(0);
407  legend->SetFillStyle(0);
408  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
409  if(!manytomany) legend->AddEntry(ProfileX_1[0], analysis1+" "+viewStr+" view");
410  else {
411  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
412  TString fbStr = Form("FB %i", fiberbrightness);
413  legend->AddEntry(ProfileX_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr);
414  }
415  }
416  if(fb) {
417  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
418  TString fbStr = Form("FB %i", fiberbrightness);
419  legend->AddEntry(ProfileX_2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr);
420  }
421  } else {
422  legend->AddEntry(ProfileX_2[0], analysis2+" "+viewStr+" view");
423  }
424 
425  legend->Draw();
426 
427  gPad->Print("Images/uncorrected_PEcm_vs_w_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
428 
429  if(manytomany){
430  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
431  new TCanvas;
432  ProfileX_1[fiberbrightness]->SetMinimum(0);
433  ProfileX_1[fiberbrightness]->SetMaximum(60);
434  ProfileX_1[fiberbrightness]->SetLineColor(kBlack);
435  ProfileX_1[fiberbrightness]->SetLineStyle(1);
436  ProfileX_1[fiberbrightness]->SetMarkerColor(kBlack);
437  ProfileX_1[fiberbrightness]->Draw();
438  ProfileX_2[fiberbrightness]->SetLineColor(kRed);
439  ProfileX_2[fiberbrightness]->SetMarkerColor(kRed);
440  ProfileX_2[fiberbrightness]->Draw("same");
441  TString fbstr;
442  fbstr.Form("fb%i", fiberbrightness);
443  TLegend * oneFBlegend = new TLegend(0.12, 0.7, 0.5, 0.88);
444  oneFBlegend->SetBorderSize(0);
445  oneFBlegend->AddEntry(ProfileX_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
446  oneFBlegend->AddEntry(ProfileX_2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
447  oneFBlegend->Draw();
448  gPad->SaveAs("Images/uncorrected_PEcm_vs_w_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view_"+fbstr+".pdf");
449  }
450  }
451 
452  //ratios
453  new TCanvas;
454  gPad->SetGrid();
455  TString hrationame[9];
456  TH1D* hratio[9];
457  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
458  TString fbStr;
459  if(fb) fbStr.Form("_fb%i", fiberbrightness);
460  hrationame[fiberbrightness].Form("hratio%s_%s_%s%s", det.Data(), mcdata.Data(), viewstr.Data(), fbStr.Data());
461  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()));
462  }
463 
464  if(fb) {
465  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
466  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
467  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
468  }
469  } else {
470  Ratiohist(ProfileX_1[0], ProfileX_2[0], hratio[0], kRed, 24, "W (cm)", "Ratio", 0.85, 1.10);
471  }
472 
473  TLegend *ratiolegend = new TLegend(.12,.62,.64,.90);
474  ratiolegend->SetBorderSize(0);
475  ratiolegend->SetFillStyle(0);
476  //ratiolegend->SetTextSize(0.04);
477  ratiolegend->AddEntry((TObject*)0, det+" "+mcdata," ");
478  if(fb) {
479  for(int fiberbrightness=0;fiberbrightness < 9;fiberbrightness++){
480  TString fbStr = Form(" FB %i", fiberbrightness);
481  ratiolegend->AddEntry(hratio[fiberbrightness],analysis2+" to "+analysis1+" "+viewStr+" view"+fbStr);
482  }
483  } else {
484  ratiolegend->AddEntry(hratio[0],analysis2+" to "+analysis1+" "+viewStr+" view");
485  }
486  ratiolegend->Draw();
487  gPad->Print("Images/uncorrected_PEcm_vs_w_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
488  }//end of loop over views
489 } // end of uncorrected_PE()
490 //-------------------------------------------------------------------------------------------------------------------------------
491 void uncorrected_PE_tricellvstrajectory(char* c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2){
492  uncorrected_PE(c1, c1, det, mcdata, analysis1, analysis2, hist1, hist2);
493 }//end of uncorrected_PE_tricellvstrajectory()
494 //-------------------------------------------------------------------------------------------------------------------------------
495 void thresholdshadowing(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var){
496  std::cout << "Begin threshold/shadowing correction files" << std::endl;
497  TFile* f1 = new TFile(c1, "read");
498  TFile* f2 = new TFile(c2, "read");
499  std::cout << "Found both files" << std::endl;
500  int nbrightnesses = 9;
501  for(int view = 0; view < 2; ++view){
502  const TString viewstr = view ? "y" : "x";
503  const TString viewStr = view ? "Y" : "X";
504 
505  TH1D* correction_1[9];
506  TH1D* corrections_to_compare[9];
507  TH1D* correction_2;
508 
509  f1->cd();
510  if(!manytomany){
511  if(var=="W")correction_1[0] = (TH1D*)f1->Get("a"+viewstr+"view");
512  else if(var=="Cell") correction_1[0] = (TH1D*)f1->Get("c"+viewstr+"view");
513  correction_1[0]->SetLineColor(kBlack);
514  std::cout << "Found non-fb histogram" << std::endl;
515  }
516 
517  if(fb){
518  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++) {
519  const TString fbStr;
520  fbStr.Form("_fb%i", fiberbrightness);
521 
522  if(manytomany) {
523  if(var=="W")correction_1[fiberbrightness] = (TH1D*)f1->Get("a"+viewstr+"view"+fbStr);
524  else if(var=="Cell")correction_1[fiberbrightness] = (TH1D*)f1->Get("c"+viewstr+"view"+fbStr);
525  correction_1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
526  correction_1[fiberbrightness]->SetLineStyle(kDotted);
527  }
528 
529  if(var=="W")corrections_to_compare[fiberbrightness] = (TH1D*)f2->Get("a"+viewstr+"view"+fbStr);
530  else if(var=="Cell")corrections_to_compare[fiberbrightness] = (TH1D*)f2->Get("c"+viewstr+"view"+fbStr);
531  corrections_to_compare[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
532  }
533  }
534  else {
535  f2->cd();
536  correction_2 = (TH1D*)f2->Get("a"+viewstr+"view");
537  correction_2->SetLineColor(kRed);
538  }
539  std::cout << "Found fb histograms" << std::endl;
540 
541  TLegend *legend = new TLegend(.5,.62,.88,.90);
542  legend->SetBorderSize(0);
543  legend->SetFillStyle(0);
544  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
545  if(!manytomany) legend->AddEntry(correction_1[0], analysis1+" "+viewStr+" view");
546  else {
547  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
548  const TString fbStr;
549  fbStr.Form("FB %i", fiberbrightness);
550  legend->AddEntry(correction_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr);
551  }
552  }
553  if(fb){
554  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
555  const TString fbStr;
556  fbStr.Form("FB %i", fiberbrightness);
557  legend->AddEntry(corrections_to_compare[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr);
558  }
559  }
560  else legend->AddEntry(correction_2, analysis2+" "+viewStr+" view");
561 
562  new TCanvas;
563  gPad->SetGrid();
564  if(correction_1[0]->FindObject("stats")){
565  TPaveStats* pave1 = (TPaveStats*)correction_1[0]->FindObject("stats");
566  pave1->Delete();
567  }
568  if(!fb && correction_2->FindObject("stats")){
569  TPaveStats* pave2 = (TPaveStats*)correction_2->FindObject("stats");
570  pave2->Delete();
571  }
572  correction_1[0]->SetMinimum(1.1);
573  correction_1[0]->SetMaximum(1.5);
574  correction_1[0]->Draw();
575  if(fb){
576  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; fiberbrightness++){
577  corrections_to_compare[fiberbrightness]->Draw("same");
578  if(manytomany && fiberbrightness > 0) correction_1[fiberbrightness]->Draw("same");
579  }
580  }
581  else correction_2->Draw("same");
582  TString labelx;
583  if(var=="W") labelx = "W (cm)";
584  if(var=="Cell") labelx = "Cell";
585 
586  float ylow = 1.0;
587  float yhigh = 1.8;
588  /*if(view == 0) {
589  ylow = 1.2;
590  yhigh = 1.5;
591  }*/
592 
593  //std::cout << "Drawing hist... " << std::endl;
594  //std::cout << "title: " << labelx.Data();
595  //std::cout << " bins:" << correction_1->GetXaxis()->GetBinCenter(1);
596  //std::cout << "-" << correction_1->GetXaxis()->GetBinCenter(correction_1->GetXaxis()->GetNbins());
597  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
598  //std::cout << "Done!" << std::endl;
599 
600  legend->Draw();
601  gPad->Print("Images/thresholdshadowing_vs_"+var+"_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
602 
603  if(manytomany){
604  for(int fiberbrightness = 0; fiberbrightness < nbrightnesses; ++fiberbrightness){
605  new TCanvas;
606  correction_1[fiberbrightness]->SetMinimum(1.1);
607  correction_1[fiberbrightness]->SetMaximum(1.6);
608  correction_1[fiberbrightness]->Draw();
609  corrections_to_compare[fiberbrightness]->Draw("same");
610  TString fbstr;
611  fbstr.Form("fb%i", fiberbrightness);
612  TLegend * oneFBlegend = new TLegend(0.5, 0.7, 0.9, 0.9);
613  oneFBlegend->AddEntry(correction_1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
614  oneFBlegend->AddEntry(corrections_to_compare[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
615  oneFBlegend->Draw();
616  gPad->SaveAs("Images/thresholdshadowing_vs_"+var+"_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view_"+fbstr+".pdf");
617  }
618  }
619 
620  new TCanvas;
621  gPad->SetGrid();
622 
623  //ratios
624  TH1D* hratio[9];
625  if(fb){
626  TLegend *legend = new TLegend(.5,.62,.88,.90);
627  legend->SetBorderSize(0);
628  legend->SetFillStyle(0);
629  //legend->SetTextSize(0.04);
630  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
631 
632  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
633  const TString fbStr = Form("FB %i", fiberbrightness);
634  const TString hrationame = Form("hratio%s_%s_%s_FB%i_%s", det.Data(), mcdata.Data(), viewstr.Data(), fiberbrightness, var.Data());
635  std::cout << hrationame << std::endl;
636  std::cout << correction_1[0]->GetXaxis()->GetNbins() << " bins, ";
637  std::cout << "first bin center " << correction_1[0]->GetBinCenter(0);
638  std::cout << ", last bin center " << correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()) << std::endl;
639  hratio[fiberbrightness] = new TH1D(hrationame, "",correction_1[0]->GetXaxis()->GetNbins(), correction_1[0]->GetBinCenter(0), correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()));
640  if(fiberbrightness==0){
641  std::cout << "FB=0" << std::endl;
642  hratio[0]->SetMinimum(0.95);
643  hratio[0]->SetMaximum(1.1);
644  }
645  std::cout << "made hratio[" << fiberbrightness << "], ";
646  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
647  else Ratiohist(correction_1[0], corrections_to_compare[fiberbrightness], hratio[fiberbrightness], fbcolors[fiberbrightness], 24, labelx.Data(), "Correction factor ratio", 0.9, 1.2);
648  std::cout << "and filled it, ";
649  legend->AddEntry(hratio[fiberbrightness],analysis2+" "+fbStr+" to "+analysis1+" "+viewStr+" view");
650  std::cout << "and added it to legend." <<std::endl;
651  }
652  legend->Draw();
653  } else {
654  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), viewstr.Data(), var.Data());
655  hratio[0] = new TH1D(hrationame, "",correction_1[0]->GetXaxis()->GetNbins(), correction_1[0]->GetBinCenter(0), correction_1[0]->GetBinCenter(correction_1[0]->GetXaxis()->GetNbins()));
656  Ratiohist(correction_1, correction_2, hratio[0], kRed, 24, labelx.Data(), "Correction factor ratio", 1.8, 1.2);
657 
658  TLegend *legend = new TLegend(.5,.62,.88,.90);
659  legend->SetBorderSize(0);
660  legend->SetFillStyle(0);
661  legend->SetTextSize(0.04);
662  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
663  legend->AddEntry(hratio[0],analysis2+" to "+analysis1+" "+viewStr+" view");
664  legend->Draw();
665  }//end if(!fb)
666 
667  gPad->Print("Images/thresholdshadowinge_vs_"+var+"_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+viewstr+"view.pdf");
668 
669  }//end of loop over views
670 }//end of thresholdshadowing()
671 //-------------------------------------------------------------------------------------------------------------------------------
672 void corrected_PE(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell){
673  int nbrightnesses = 9;
674 
675  if(nddatamuoncatcher){
676  if(plane.Atoi() > 3 || plane.Atoi() < 2){
677  std::cout << "It is a ND muon catcher data file and has only 002, 003 planes. Plane " << plane << " is out of range. " << std::endl;
678  abort();
679  }
680  if(cell.Atoi() > 95 || cell.Atoi() < 0){
681  std::cout << "It is a Data file. Cell " << cell << " is out of range. " << std::endl;
682  abort();
683  }
684  }
685  if(det.Data()=="nd"){
686  if(mcdata.Data()=="mc"){
687  if(plane.Atoi() > 3 || plane.Atoi() < 0){
688  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
689  abort();
690  }
691  if(cell.Atoi() > 95 || cell.Atoi() < 0){
692  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
693  abort();
694  }
695  }//end of MC condition
696  if(mcdata.Data()=="data"){
697  if(plane.Atoi() > 191 || plane.Atoi() < 0){
698  std::cout << "It is a ND data file and has only 191 main body planes. Plane " << plane << " out of range. " << std::endl;
699  abort();
700  }
701  if(cell.Atoi() > 95 || cell.Atoi() < 0){
702  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
703  abort();
704  }
705  }//end of data condition
706  }//end of the ND condition
707 
708  if(det.Data()=="fd"){
709  if(mcdata.Data()=="mc"){
710  if(plane.Atoi() > 1 || plane.Atoi() < 0){
711  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
712  abort();
713  }
714  if(cell.Atoi() > 383 || cell.Atoi() < 0){
715  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
716  abort();
717  }
718  }//end of MC condition
719  else if(mcdata.Data()=="data"){
720  if(plane.Atoi() > 895 || plane.Atoi() < 0){
721  std::cout << "It is a Data file. Plane " << plane << " out of range. " << std::endl;
722  abort();
723  }
724  if(cell.Atoi() > 383 || cell.Atoi() < 0){
725  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
726  abort();
727  }
728  }//end of data condition
729  }//end of FD condition
730 
731  TFile* f1 = new TFile(c1, "read");
732  TFile* f2 = new TFile(c2, "read");
733 
734  TF1* fit1[9];
735  TF1* fit2[9];
736 
737  f1->cd();
738  std::cout << "getting " << c1 << ": plane_" << plane << "/Atten_Fit_" << plane << "_" << cell << std::endl;
739  TH1F* h1[9];
740  if(!manytomany){
741  h1[0] = (TH1F*)f1->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
742  h1[0]->ls();
743  Color_t color1 = kBlack;
744  if(mcdata=="mc"){
745  h1[0]->GetListOfFunctions()->Print();
746  if(plane.Atoi()==0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
747  else if(plane.Atoi()==1) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
748  else if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
749  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
750  if(fit1[0]) fit1[0]->SetLineColor(color1);
751  }//end of MC condition
752  else if(mcdata=="data"){
753  if(nddatamuoncatcher){
754  if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
755  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
756  if(fit1[0]) fit1[0]->SetLineColor(color1);
757  }//end of muon catcher condition
758  else if(plane.Atoi()%2==0) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
759  else if(plane.Atoi()%2!=0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
760  if(fit1[0]) fit1[0]->SetLineColor(color1);
761  }//end of data condition
762  if(!fit1[0]) std::cout << "Can't find fit1[0]" << std::endl;
763  } else {
764  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
765  //std::cout << "making h1[" << fiberbrightness << "]" << std::endl;
766  std::cout << "Filling h1[" << fiberbrightness << "]" << std::endl;
767 
768  TString fbStr = Form("fb%i_", fiberbrightness);
769  h1[fiberbrightness] = (TH1F*)f1->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell);
770 
771  std::cout << "Getting " << (fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell).Data() << std::endl;
772  if(!h1[fiberbrightness]) std::cout << "Can't find h1 " << fiberbrightness << std::endl;
773  h1[fiberbrightness]->GetListOfFunctions()->Print();
774  //h1[fiberbrightness]->SetMarkerStyle[24];
775 
776  //h2[fiberbrightness]->ls();
777  if(plane.Atoi()%2==0){
778  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitY");//X
779  if(fit1[fiberbrightness]) std::cout << "using fitY" <<std::endl;
780  } else {
781  std::cout << "looking for fitX" << std::endl;
782  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");
783  if(fit1[fiberbrightness]) std::cout << "using fitX" << std::endl;
784  }
785  // fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
786  // fit1[fiberbrightness]->SetLineStyle(kDotted);
787 
788  if(!fit1[fiberbrightness])
789  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");
790  if(!fit1[fiberbrightness])
791  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitY");
792  if(!fit1[fiberbrightness])
793  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitMuX");
794 
795  if(fit1[fiberbrightness]){
796  fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
797  fit1[fiberbrightness]->SetLineStyle(kDotted);
798  } else {
799  std::cout << "can't find fit1[" << fiberbrightness << "]"<< std::endl;
800  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("rollFitY");
801  if(fit1[fiberbrightness]) fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
802  }
803  }
804  }
805 
806  f2->cd();
807  TH1F* h2[9];
808  //f2->Get("fb0_plane_000")->ls();
809  if(fb){
810  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
811  std::cout << "making h2[" << fiberbrightness << "]" << std::endl;
812  TString fbStr = Form("fb%i_", fiberbrightness);
813  h2[fiberbrightness] = (TH1F*)f2->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cell);
814  //if(fiberbrightness==0)
815  //h2[fiberbrightness]->GetListOfFunctions()->Print();
816  if(!h2[fiberbrightness]) std::cout << "Can't find h2 " << fiberbrightness << std::endl;
817  //h2[fiberbrightness]->ls();
818  if(plane.Atoi()%2==0){
819  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitY");//X
820  std::cout << "using fitY" <<std::endl;
821  } else {
822  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");
823  }
824 
825  if(!fit2[fiberbrightness])
826  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");
827  if(!fit2[fiberbrightness])
828  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitY");
829  if(!fit2[fiberbrightness])
830  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitMuX");
831 
832  if(fit2[fiberbrightness]) fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
833  else {//std::cout << "can't find fit2[" << fiberbrightness << "]"<< std::endl;
834  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("rollFitY");
835  if(fit2[fiberbrightness]) fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
836  }
837  }
838  } else {
839  h2[0] = (TH1F*)f2->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
840  }
841  Color_t color1 = kBlack;
842  Color_t color2 = kRed;
843 
844  if(mcdata=="mc" && !fb){
845  if(plane.Atoi()==0) fit2[0] = (TF1*)h2[0]->FindObject("fitX");
846  else if(plane.Atoi()==1) fit2[0] = (TF1*)h2[0]->FindObject("fitY");
847  else if(plane.Atoi()==2) fit2[0] = (TF1*)h2[0]->FindObject("fitMuX");
848  else if(plane.Atoi()==3) fit2[0] = (TF1*)h2[0]->FindObject("fitMuY");
849  if(fit2[0]) fit2[0]->SetLineColor(color2);
850  }//end of MC condition
851  else if(mcdata=="data" && !fb){
852  if(nddatamuoncatcher){
853  if(plane.Atoi()==2) fit2[0] = (TF1*)h2[0]->FindObject("fitMuX");
854  else if(plane.Atoi()==3) fit2[0] = (TF1*)h2[0]->FindObject("fitMuY");
855  if(fit2[0]) fit2[0]->SetLineColor(color2);
856  }//end of muon catcher condition
857  else if(plane.Atoi()%2==0) fit2[0] = (TF1*)h2[0]->FindObject("fitY");
858  else if(plane.Atoi()%2!=0) fit2[0] = (TF1*)h2[0]->FindObject("fitX");
859  if(fit2[0]) fit2[0]->SetLineColor(color2);
860  }//end of data condition
861 
862  new TCanvas;
863  gPad->SetGrid();
864  h1[0]->Draw();
865  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
866  if(!manytomany) h1[0]->SetLineColor(color1);
867  h1[0]->SetMarkerStyle(24);
868  h1[0]->SetMaximum(50);//old max: 120, then 60
869  if(fb) {
870  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
871  std::cout << "Drawing fiberbrightness = " << fiberbrightness << std::endl;
872  if(manytomany) {
873  h1[fiberbrightness]->Draw("same");
874  h1[fiberbrightness]->SetLineStyle(kDotted);
875  h1[fiberbrightness]->SetMarkerColor(fbcolors[fiberbrightness]);
876  h1[fiberbrightness]->SetMarkerStyle(25);
877  h1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
878  }
879  h2[fiberbrightness]->Draw("same");
880  h2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
881  h2[fiberbrightness]->SetMarkerColor(fbcolors[fiberbrightness]);
882  h2[fiberbrightness]->SetMarkerStyle(24);
883  }
884  } else {
885  h2[0]->Draw("same");
886  h2[0]->SetLineColor(color2);
887  h2[0]->SetMarkerColor(color2);
888  h2[0]->SetMarkerStyle(24);
889  }
890 
891  TLegend *legend = new TLegend(.12,.4,.64,.90);
892  legend->SetBorderSize(0);
893  legend->SetFillStyle(0);
894  //legend->SetTextSize(0.04);
895  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
896  if(!manytomany)legend->AddEntry(h1[0], analysis1);
897  else {
898  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
899  TString fbStr = Form(", FB %i", fiberbrightness);
900  legend->AddEntry(h1[fiberbrightness], analysis1+fbStr);
901  }
902  }
903  if(fb) {
904  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
905  TString fbStr = Form(", FB %i", fiberbrightness);
906  legend->AddEntry(h2[fiberbrightness], analysis2+fbStr);
907  }
908  } else {
909  legend->AddEntry(h2[0], analysis2);
910  }
911  legend->Draw();
912  gPad->Print("Images/Atten_Fit_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+plane+"_"+cell+".pdf");
913 
914  if(manytomany) {
915  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
916  new TCanvas();
917  h1[fiberbrightness]->Draw();
918  h2[fiberbrightness]->Draw("same");
919  TLegend * oneFBlegend = new TLegend(0.5, 0.7, 0.9, 0.9);
920  const TString viewStr = (plane.Atoi()%2==0) ? "Y" : "X";
921  oneFBlegend->AddEntry(h1[fiberbrightness], analysis1+" "+viewStr+" view "+fbStr, "l");
922  oneFBlegend->AddEntry(h2[fiberbrightness], analysis2+" "+viewStr+" view "+fbStr, "l");
923  oneFBlegend->Draw();
924  TString fbstr;
925  fbstr.Form("fb%i", fiberbrightness);
926  gPad->Print("Images/Atten_Fit_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_"+plane+"_"+cell+"_"+fbstr+".pdf");
927  }
928  }
929 
930  //ratio of data points
931  new TCanvas;
932  gPad->SetGrid();
933  const TString hrationame = Form("hratio%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
934  TH1D* hratio[9];
935  TString hrationamearr[9];
936  if(fb) {
937  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
938  hrationamearr[fiberbrightness] = Form("%s_fb%i", hrationame.Data(), fiberbrightness);
939  //hratio[fiberbrightness] = new TH1D(hrationamearr[fiberbrightness], "",h1->GetXaxis()->GetNbins(), h1->GetBinCenter(0), h1->GetBinCenter(h1->GetXaxis()->GetNbins()));
940  hratio[fiberbrightness] = new TH1D(hrationamearr[fiberbrightness], "",h1[0]->GetXaxis()->GetNbins(), -800, 800);
941  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
942  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
943  }
944  } else {
945  hratio[0] = = new TH1D(hrationame, "",h1[0]->GetXaxis()->GetNbins(), h1[0]->GetBinCenter(0), h1[0]->GetBinCenter(h1[0]->GetXaxis()->GetNbins()));
946  Ratiohist(h1[0], h2[0], hratio[0], kRed, 24, "W (cm)", "Ratio of Mean PE/cm data points", 0.3, 1.0);
947  }
948 
949  TLegend *legend = new TLegend(.12,.62,.64,.90);
950  legend->SetBorderSize(0);
951  legend->SetFillStyle(0);
952  legend->SetTextSize(0.04);
953  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
954  if(fb) {
955  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
956  TString fbStr = Form(" FB %i", fiberbrightness);
957  legend->AddEntry(hratio[fiberbrightness],analysis2+fbStr+" to "+analysis1);
958  }
959  } else {
960  legend->AddEntry(hratio[0],analysis2+" to "+analysis1);
961  }
962  legend->Draw();
963  gPad->Print("Images/Atten_Fit_ratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");
964 
965  //ratio of fitting curve
966  new TCanvas;
967  gPad->SetGrid();
968  const TString hratiofitname = Form("hratiofit%s_%s_%s_%s", det.Data(), mcdata.Data(), plane.Data(), cell.Data());
969  TString hratiofitnamearr[9];
970  float x_limit;
971  int x_bins;
972  if(det=="nd"){
973  x_limit = 250.;
974  x_bins = 2*x_limit;
975  }
976  else if(det=="fd"){
977  x_limit = 800.;
978  x_bins = 2*x_limit;
979  }
980 
981  TH1D* hfit = new TH1D("hfit", "", x_bins, -x_limit, +x_limit);
982  std::cout << "l592 about to use HistogramAttr1D" << std::endl;
983  HistogramAttr1D(hfit, "W (cm)", "Ratio of Mean PE/cm fitting curves", -x_limit, +x_limit, 0.9, 1.2, kWhite);//old limits: 0.55, 0.68
984  hfit->Draw();
985  TH1D* hratiofit[9];
986  if(fb){
987  std::cout << "running over " << nbrightnesses << "brightnesses" << std::endl;
988  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
989  hratiofitnamearr[fiberbrightness] = Form("%s_fb%i", hratiofitname.Data(), fiberbrightness);
990  hratiofit[fiberbrightness] = new TH1D(hratiofitnamearr[fiberbrightness], "", x_bins, -x_limit, +x_limit);
991  Ratiofit(fit1[fiberbrightness], 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
992  }
993  } else {
994  hratiofit[0] = new TH1D(hratiofitname, "", x_bins, -x_limit, +x_limit);
995  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);
996  }
997 
998  TLegend *legend = new TLegend(.12,.62,.64,.90);
999  legend->SetBorderSize(0);
1000  legend->SetFillStyle(0);
1001  legend->SetTextSize(0.04);
1002  legend->AddEntry((TObject*)0, det+" "+mcdata+" : Plane "+plane+" Cell "+cell," ");
1003 
1004  if(fb) {
1005  for(int fiberbrightness=0;fiberbrightness<nbrightnesses;fiberbrightness++){
1006  TString fbStr = Form(" FB %i", fiberbrightness);
1007  legend->AddEntry(hratiofit[fiberbrightness],analysis2+fbStr+" to "+analysis1);
1008  }
1009  } else {
1010  legend->AddEntry(hratiofit,analysis2+" to "+analysis1);
1011  }
1012  legend->Draw();
1013  gPad->Print("Images/Atten_Fit_fitratio_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+"_plane"+plane+"_cell"+cell+".pdf");
1014 }//end of corrected_PE()
1015 
1016 //---------------------------------------------------------------------------------------------------
1017 void corrected_PE_mean(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell){
1018  static const int nBrightnesses = 9;
1019  int ncells = 384;
1020  int ncells1 = ncells;
1021  int ncells2 = ncells;
1022 
1023  float x_limit;
1024  int x_bins;
1025  if(det=="nd"){
1026  x_limit = 250.;
1027  x_bins = 2*x_limit;
1028  }
1029  else if(det=="fd"){
1030  x_limit = 800.;
1031  x_bins = 2*x_limit;
1032  }
1033 
1034  TH1D* hfit1[9];
1035  TH1D* hfit2[9];
1036 
1037  for(int i = 0; i<nBrightnesses; ++i){
1038  TString fitname1 = Form("hfit1_FB%i", i);
1039  TString fitname2 = Form("hfit2_FB%i", i);
1040  hfit1[i] = new TH1D(fitname1, "", x_bins, -x_limit, +x_limit);
1041  hfit2[i] = new TH1D(fitname2, "", x_bins, -x_limit, +x_limit);
1042  }
1043 
1044  TH1F* hMean1[9];
1045  TH1F* hMean2[9];
1046 
1047  if(nddatamuoncatcher){
1048  if(plane.Atoi() > 3 || plane.Atoi() < 2){
1049  std::cout << "It is a ND muon catcher data file and has only 002, 003 planes. Plane " << plane << " is out of range. " << std::endl;
1050  abort();
1051  }
1052  if(cell.Atoi() > 95 || cell.Atoi() < 0){
1053  std::cout << "It is a Data file. Cell " << cell << " is out of range. " << std::endl;
1054  abort();
1055  }
1056  }
1057  if(det.Data()=="nd"){
1058  if(mcdata.Data()=="mc"){
1059  if(plane.Atoi() > 3 || plane.Atoi() < 0){
1060  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
1061  abort();
1062  }
1063  if(cell.Atoi() > 95 || cell.Atoi() < 0){
1064  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
1065  abort();
1066  }
1067  }//end of MC condition
1068  if(mcdata.Data()=="data"){
1069  if(plane.Atoi() > 191 || plane.Atoi() < 0){
1070  std::cout << "It is a ND data file and has only 191 main body planes. Plane " << plane << " out of range. " << std::endl;
1071  abort();
1072  }
1073  if(cell.Atoi() > 95 || cell.Atoi() < 0){
1074  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
1075  abort();
1076  }
1077  }//end of data condition
1078  }//end of the ND condition
1079 
1080  if(det.Data()=="fd"){
1081  if(mcdata.Data()=="mc"){
1082  if(plane.Atoi() > 1 || plane.Atoi() < 0){
1083  std::cout << "It is a Monte Carlo file. Plane " << plane << " out of range. " << std::endl;
1084  abort();
1085  }
1086  if(cell.Atoi() > 383 || cell.Atoi() < 0){
1087  std::cout << "It is a Monte Carlo file. Cell " << cell << " out of range. " << std::endl;
1088  abort();
1089  }
1090  }//end of MC condition
1091  else if(mcdata.Data()=="data"){
1092  if(plane.Atoi() > 895 || plane.Atoi() < 0){
1093  std::cout << "It is a Data file. Plane " << plane << " out of range. " << std::endl;
1094  abort();
1095  }
1096  if(cell.Atoi() > 383 || cell.Atoi() < 0){
1097  std::cout << "It is a Data file. Cell " << cell << " out of range. " << std::endl;
1098  abort();
1099  }
1100  }//end of data condition
1101  }//end of FD condition
1102 
1103  TFile* f1 = new TFile(c1, "read");
1104  TFile* f2 = new TFile(c2, "read");
1105  TF1* fit1[9];
1106  TF1* fit2[9];
1107 
1108  f1->cd();
1109  std::cout << "getting " << c1 << ": plane_" << plane << "/Atten_Fit_" << plane << "_" << cell << std::endl;
1110  TH1F* h1[9];
1111  if(!manytomany){
1112  h1[0] = (TH1F*)f1->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cell);
1113  h1[0]->ls();
1114  Color_t color1 = kBlack;
1115  if(mcdata=="mc"){
1116  h1[0]->GetListOfFunctions()->Print();
1117  if(plane.Atoi()==0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
1118  else if(plane.Atoi()==1) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
1119  else if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
1120  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
1121  if(fit1[0]) fit1[0]->SetLineColor(color1);
1122  }//end of MC condition
1123  else if(mcdata=="data"){
1124  if(nddatamuoncatcher){
1125  if(plane.Atoi()==2) fit1[0] = (TF1*)h1[0]->FindObject("fitMuX");
1126  else if(plane.Atoi()==3) fit1[0] = (TF1*)h1[0]->FindObject("fitMuY");
1127  if(fit1[0]) fit1[0]->SetLineColor(color1);
1128  }//end of muon catcher condition
1129  else if(plane.Atoi()%2==0) fit1[0] = (TF1*)h1[0]->FindObject("fitY");
1130  else if(plane.Atoi()%2!=0) fit1[0] = (TF1*)h1[0]->FindObject("fitX");
1131  if(fit1[0]) fit1[0]->SetLineColor(color1);
1132  }//end of data condition
1133  if(!fit1[0]) std::cout << "Can't find fit1[0]" << std::endl;
1134  } else {
1135  for(int fiberbrightness=0;fiberbrightness<nBrightnesses;fiberbrightness++){
1136  //std::cout << "making h1[" << fiberbrightness << "]" << std::endl;
1137  //std::cout << "Filling h1[" << fiberbrightness << "] with cell " << cellname << std::endl;
1138  TString fbStr = Form("fb%i_", fiberbrightness);
1139 
1140  for(int cellnumber = 0; cellnumber < ncells; ++cellnumber){
1141  TString cellname;
1142  if(cellnumber < 10) cellname.Form("00%i", cellnumber);
1143  else if(cellnumber < 100) cellname.Form("0%i", cellnumber);
1144  else cellname.Form("%i", cellnumber);
1145  //std::cout << "making h1[" << fiberbrightness << "]" << std::endl;
1146  std::cout << "Filling h1[" << fiberbrightness << "] with cell " << cellname << std::endl;
1147 
1148  h1[fiberbrightness] = (TH1F*)f1->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cellname);
1149  if(!h1[fiberbrightness]) std::cout << "Can't find h1 " << fiberbrightness << std::endl;
1150  if(cellnumber == 0) hMean1[fiberbrightness] = (TH1F*)h1[fiberbrightness]->Clone();
1151  else hMean1[fiberbrightness]->Add(h1[fiberbrightness]);
1152  //h1[fiberbrightness]->GetListOfFunctions()->Print();
1153 
1154  if(plane.Atoi()%2==0){
1155  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitY");
1156  std::cout << "using fitY" <<std::endl;
1157  } else {
1158  std::cout << "looking for fitX" << std::endl;
1159  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");
1160  if(fit1[fiberbrightness]) std::cout << "Got fitX" << std::endl;
1161  }
1162  if(!fit1[fiberbrightness]) fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitY");
1163  if(!fit1[fiberbrightness]) fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitX");
1164 
1165  if(!fit1[fiberbrightness]) fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitMuY");
1166  if(!fit1[fiberbrightness]) fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("fitMuX");
1167 
1168  if(fit1[fiberbrightness]){
1169  fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1170  fit1[fiberbrightness]->SetLineStyle(kDotted);
1171  for(int binnumber = 1; binnumber <= x_bins; ++binnumber){
1172  float w = hfit1[fiberbrightness]->GetXaxis()->GetBinCenter(binnumber);
1173  hfit1[fiberbrightness]->AddBinContent(binnumber, fit1[fiberbrightness]->Eval(w));
1174  }
1175  } else {
1176  std::cout << "can't find fit1[" << fiberbrightness << "]"<< std::endl;
1177  ncells1 -= 1;
1178  fit1[fiberbrightness] = (TF1*)h1[fiberbrightness]->FindObject("rollFitY");
1179  if(fit1[fiberbrightness]) fit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1180  }
1181  }
1182  hMean1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1183  hfit1[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1184  hfit1[fiberbrightness]->SetLineStyle(kDotted);
1185  }
1186  }
1187 
1188  f2->cd();
1189  TH1F* h2[9];
1190  //f2->Get("fb0_plane_000")->ls();
1191  if(fb){
1192  for(int fiberbrightness=0;fiberbrightness<nBrightnesses;fiberbrightness++){
1193 
1194  TString fbStr = Form("fb%i_", fiberbrightness);
1195  for(int cellnumber = 0; cellnumber < ncells; ++cellnumber){
1196  TString cellname;
1197  if(cellnumber < 10) cellname.Form("00%i", cellnumber);
1198  else if(cellnumber < 100) cellname.Form("0%i", cellnumber);
1199  else cellname.Form("%i", cellnumber);
1200  //std::cout << "making h1[" << fiberbrightness << "]" << std::endl;
1201  std::cout << "Filling h2[" << fiberbrightness << "] with cell " << cellname << std::endl;
1202 
1203  h2[fiberbrightness] = (TH1F*)f2->Get(fbStr+"plane_"+plane+"/Atten_Fit_"+fbStr+plane+"_"+cellname);
1204  if(!h2[fiberbrightness]) std::cout << "Can't find h2 " << fiberbrightness << std::endl;
1205  if(cellnumber == 0) hMean2[fiberbrightness] = (TH1F*)h2[fiberbrightness]->Clone();
1206  else hMean2[fiberbrightness]->Add(h2[fiberbrightness]);
1207 
1208  //h2[fiberbrightness]->GetListOfFunctions()->Print();
1209  //if(!h2[fiberbrightness]) std::cout << "Can't find h2 " << fiberbrightness << std::endl;
1210  //h2[fiberbrightness]->ls();
1211  if(plane.Atoi()%2==0){
1212  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitY");
1213  std::cout << "using fitY" <<std::endl;
1214  } else {
1215  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");
1216  std::cout << "using fitX" <<std::endl;
1217  }
1218  if(!fit2[fiberbrightness]) fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");
1219  if(!fit2[fiberbrightness]) fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitX");
1220  if(!fit2[fiberbrightness]) fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("fitMuX");
1221 
1222  if(fit2[fiberbrightness]) {
1223  for(int binnumber = 1; binnumber <= x_bins; ++binnumber){
1224  float w = hfit2[fiberbrightness]->GetXaxis()->GetBinCenter(binnumber);
1225  hfit2[fiberbrightness]->AddBinContent(binnumber, fit2[fiberbrightness]->Eval(w));
1226  }
1227  fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1228  } else {
1229  std::cout << "can't find fit2[" << fiberbrightness << "]"<< std::endl;
1230  ncells2 -= 1;
1231  fit2[fiberbrightness] = (TF1*)h2[fiberbrightness]->FindObject("rollFitY");
1232  if(fit2[fiberbrightness]) fit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1233  }
1234  }
1235  hMean2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1236  hfit2[fiberbrightness]->SetLineColor(fbcolors[fiberbrightness]);
1237  hfit2[fiberbrightness]->SetLineStyle(kDotted);
1238  }
1239  } else {
1240  h2[0] = (TH1F*)f2->Get("plane_"+plane+"/Atten_Fit_"+plane+"_"+cellname);
1241  }
1242 
1243  TCanvas * cMean = new TCanvas();
1244  TCanvas * cMeanRatio = new TCanvas();
1245  TCanvas * cFitMean = new TCanvas();
1246  TCanvas * cFitMeanRatio = new TCanvas();
1247 
1248  TLegend * legtemp = new TLegend(0.1, 0.75, 0.5, 0.9);
1249  legtemp->AddEntry(hfit2, analysis2, "l");
1250  legtemp->Draw();
1251 
1252  for(int brightness = 0; brightness < nBrightnesses; ++brightness){
1253  TString fbStr = Form("FB %i", brightness);
1254 
1255  std::cout << "About to scale (ncells1 = " << ncells1 << ", ncells2 = " << ncells2 << std::endl;
1256  hMean1[brightness]->Scale(1./(float)ncells);
1257  hMean2[brightness]->Scale(1./(float)ncells);
1258  hfit1[brightness]->Scale(1./(float)ncells1);
1259  hfit2[brightness]->Scale(1./(float)ncells2);
1260  std::cout << "Scaled" << std::endl;
1261 
1262  cMean->cd();
1263  hMean1[brightness]->SetLineStyle(kDotted);
1264  hMean1[brightness]->SetLineColor(fbcolors[fiberbrightness]);
1265  if(brightness == 0){
1266  hMean1[brightness]->SetTitle("PE/cm Mean;W(cm);PE/cm");
1267  hMean1[brightness]->Draw();
1268  } else
1269  hMean1[brightness]->Draw("same");
1270  hMean2[brightness]->SetLineColor(fbcolors[fiberbrightness]);
1271  hMean2[brightness]->Draw("same");
1272  cMean->SaveAs("Images/FB0Mean.pdf");
1273 
1274  legtemp->AddEntry(hMean1[brightness], analysis1+" "+fbStr, "l");
1275  legtemp->AddEntry(hMean2[brightness], analysis2+" "+fbStr, "l");
1276 
1277  if(brightness == nBrightnesses-1){
1278  legtemp->Draw();
1279  cMean->Print("Images/AttenFit_plane_mean.pdf");
1280  }
1281 
1282  // cFitMean->cd();
1283  // if(brightness == 0){
1284  // hfit1[brightness]->SetTitle("XView Mean Fits: FB 0;W (cm);PE/cm");
1285  // //hfit1[brightness]->SetLineColor(kBlack);
1286  // hfit1[brightness]->SetMaximum(60);
1287  // hfit1[brightness]->Draw();
1288  // } else
1289  // hfit1[brightness]->Draw("same");
1290  // // hfit2->SetLineColor(kRed);
1291  // hfit2[brightness]->Draw("same");
1292 
1293  // cFitMeanRatio->cd();
1294  // hfit1[brightness]->Divide(hfit2[brightness]);
1295  // hfit1[brightness]->SetTitle("Plane "+plane+" Ratio of mean fits;W (cm);");
1296  // hfit1[brightness]->SetLineStyle(kSolid);
1297  // hfit1[brightness]->SetMinimum(0.9);
1298  // hfit1[brightness]->SetMaximum(1.1);
1299  // if(brightness == 0)hfit1[brightness]->Draw();
1300  // hfit1[brightness]->Draw("same");
1301  }
1302  std::cout << "Go to canvas" << std::endl;
1303  cMean->cd();
1304  std::cout << "Draw legend" << std::endl;
1305  legtemp->Draw();
1306  std::cout << "Save canvas" << std::endl;
1307  // cMean->Print("Images/AttenFit_plane_mean.pdf");
1308  std::cout << "Saved";
1309 
1310  // cFitMean->cd();
1311  // legtemp->Draw();
1312  // cFitMean->SaveAs("AttenFit_plane"+plane+"_fit_mean.pdf");
1313 
1314  // cFitMeanRatio->cd();
1315  // legtemp->Draw();
1316  // TLine * oneline = new TLine(hfit1[0]->GetXaxis()->GetXmin(), 1., hfit1[0]->GetXaxis()->GetXmax(), 1.);
1317  // oneline->SetLineStyle(kDashed);
1318  // online->Draw();
1319  // cFitMeanRatio->SaveAs("FitRatio_xview_fb0.pdf");
1320 
1321 }//end corrected_PE_mean
1322 
1323 //---------------------------------------------------------------------------------------------------
1324 void calibrationconstants(char* c1, char* c2, TString det, TString mcdata, TString analysis1, TString analysis2){
1325  ifstream in_1;
1326  ifstream in_2;
1327 
1328  ofstream out;
1329 
1330  FILE* fout = fopen("uncalibrated_channels_"+det+"_"+mcdata+"_"+analysis1+"_vs_"+analysis2+".csv","w");
1331  fprintf(fout, "%s%50s%s\n", analysis1.Data(), " ", analysis2.Data());
1332  fprintf(fout, "channel,coeff_exp,atten_length,background,chisq, %9s, channel,coeff_exp,atten_length,background,chisq\n", " ");
1333 
1334  in_1.open(c1);
1335  in_2.open(c2);
1336 
1337 
1338  //x=coeff_exp, y=atten_length, z=background, a=edgx_low, b=edgx_high, c=coeff_low
1339  Float_t x_1, y_1, z_1, a_1, b_1, c_1, d_1, e_1, g_1;
1340  Float_t x_2, y_2, z_2, a_2, b_2, c_2, d_2, e_2, g_2;
1341  Float_t chisq = 0.2;
1342 
1343  Int_t channels;
1344  Int_t totalchannels1 = 0;
1345  Int_t totalchannels2 = 0;
1346  Int_t calibratedchannels1 = 0;
1347  Int_t calibratedchannels2 = 0;
1348  Int_t f_1;
1349  Int_t f_2;
1350  Int_t bins1;
1351 
1352  Double_t limit1D = 3.;
1353  Double_t limit2D = 1000.;
1354 
1355  if(det=="nd" && mcdata=="mc" ) { channels = 3063; }
1356  if(det=="nd" && mcdata=="data" ) { channels = 1003063; }
1357  if(det=="fd" && mcdata=="mc" ) { channels = 1383; }
1358  if(det=="fd" && mcdata=="data" ) { channels = 895383; }
1359 
1360 
1361  TH1F* hcoeffexp = new TH1F("hcoeffexp", "coeffexp difference distribution", 120, -limit1D, +limit1D );
1362  TH1F* hattenlength = new TH1F("hattenlength", "attenlength difference distribution", 600, -limit1D, +limit1D );
1363  TH1F* hbackground = new TH1F("hbackground", "background difference distribution", 120, -limit1D, +limit1D );
1364  TH1F* hchisq = new TH1F("hchisq", "chisq difference distribution", 600, -limit1D, +limit1D );
1365 
1366  TH2F* h2coeffexp = new TH2F("h2coeffexp", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1367  TH2F* h2attenlength = new TH2F("h2attenlength"," " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1368  TH2F* h2background = new TH2F("h2background", " " , 2*(Int_t)limit2D, -limit2D, limit2D, 2*(Int_t)limit2D, -limit2D, limit2D );
1369  TH2F* h2chisq = new TH2F("h2chisq"," ", 200,0.,2., 200,0.,2. );
1370 
1371  while(1){
1372  in_1 >> x_1 >> y_1 >> z_1 >> a_1 >> b_1 >> c_1 >> d_1 >> e_1 >> f_1 >> g_1;
1373  in_2 >> x_2 >> y_2 >> z_2 >> a_2 >> b_2 >> c_2 >> d_2 >> e_2 >> f_2 >> g_2;
1374  if(!in_1.good()) break;
1375  if(!in_2.good()) break;
1376 
1377  if(f_1<=channels && e_1 < chisq && e_2 < chisq){
1378  hcoeffexp->Fill((x_1 - x_2)/x_1);
1379  hattenlength->Fill((y_1 - y_2)/y_1);
1380  hbackground->Fill((z_1 - z_2)/z_1);
1381 
1382  h2coeffexp->Fill(x_2, x_1);
1383  h2attenlength->Fill(y_2, y_1);
1384  h2background->Fill(z_2, z_1);
1385  }//end of loop over channels and chisq conditions.
1386  if(f_1<=channels){
1387  totalchannels1++;
1388  totalchannels2++;
1389  if(e_1<chisq) calibratedchannels1++;
1390  if(e_2<chisq) calibratedchannels2++;
1391 
1392  hchisq->Fill((e_1 - e_2)/e_1);
1393  h2chisq->Fill(e_2, e_1);
1394 
1395  if(e_1>=chisq || e_2>=chisq){
1396  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);
1397  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;
1398  }
1399  }//end of only channels condition for chisq plots. For chisq plots we have to consider all channels.
1400  }//end of while condition.
1401 
1402  cout << " " << endl;
1403  cout << "Total channels in "+analysis1+" = " << totalchannels1 << " and calibrated channels = " << calibratedchannels1 << " = " << ((float)calibratedchannels1/(float)totalchannels1)*100. << " %." << endl;
1404  cout << "Total channels in "+analysis2+" = " << totalchannels2 << " and calibrated channels = " << calibratedchannels2 << " = " << ((float)calibratedchannels2/(float)totalchannels2)*100. << " %." << endl;
1405  cout << " " << endl;
1406 
1407  //coeffexp plots
1408  new TCanvas;
1409  gPad->SetGrid();
1410  std::cout << "l712 about to use HistogramAttr1D" << std::endl;
1411  HistogramAttr1D(hcoeffexp, "#Deltacoeffexp/coeffexp", "# of channels", 0.5*hcoeffexp->GetXaxis()->GetBinCenter(1), 0.5*hcoeffexp->GetXaxis()->GetBinCenter(hcoeffexp->GetXaxis()->GetNbins()), 0., 7e5, kBlack);
1412  hcoeffexp->Draw();
1413  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1414  legend->SetBorderSize(0);
1415  legend->SetFillStyle(0);
1416  legend->SetTextSize(0.05);
1417  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1418  legend->AddEntry(hcoeffexp,"#splitline{"+analysis1.Data()+"}{vs "+analysis2.Data()+"}","P");
1419  legend->Draw();
1420  gPad->Print("Images/Delta_coeffexp_"+analysis1.Data()+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1421 
1422  new TCanvas;
1423  gPad->SetGrid();
1424  gPad->SetLogz();
1425 
1426  HistogramAttr2D(h2coeffexp, "coeffexp_{"+analysis2+"}", "coeffexp_{"+analysis1+"}", " ", 0, 100, 0, 100);
1427  h2coeffexp->Draw("colz");
1428  gPad->Print("Images/coeffexp2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1429 
1430  //attenlength plots
1431  new TCanvas;
1432  gPad->SetGrid();
1433  HistogramAttr1D(hattenlength, "#Deltaattenlength/attenlength", "# of channels", 0.5*hattenlength->GetXaxis()->GetBinCenter(1), 0.5*hattenlength->GetXaxis()->GetBinCenter(hattenlength->GetXaxis()->GetNbins()), 0., 166500, kBlack);
1434  hattenlength->Draw();
1435  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1436  legend->SetBorderSize(0);
1437  legend->SetFillStyle(0);
1438  legend->SetTextSize(0.05);
1439  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1440  legend->AddEntry(hattenlength,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1441  legend->Draw();
1442  gPad->Print("Images/Delta_attenlength_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1443 
1444  new TCanvas;
1445  gPad->SetGrid();
1446  gPad->SetLogz();
1447  HistogramAttr2D(h2attenlength, "attenlength_{"+analysis2+"} (cm)", "attenlength_{"+analysis1+"} (cm)", " ", 0, 1000, 0, 1000);
1448  h2attenlength->Draw("colz");
1449  gPad->Print("Images/attenlength2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1450 
1451  //backgrounds plots
1452  new TCanvas;
1453  gPad->SetGrid();
1454  HistogramAttr1D(hbackground, "#Deltabackground/background", "# of channels", 0.5*hbackground->GetXaxis()->GetBinCenter(1), 0.5*hbackground->GetXaxis()->GetBinCenter(hbackground->GetXaxis()->GetNbins()), 0., 3e4, kBlack);
1455  hbackground->Draw();
1456  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1457  legend->SetBorderSize(0);
1458  legend->SetFillStyle(0);
1459  legend->SetTextSize(0.05);
1460  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1461  legend->AddEntry(hbackground,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1462  legend->Draw();
1463  gPad->Print("Images/Delta_background_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1464 
1465  new TCanvas;
1466  gPad->SetGrid();
1467  gPad->SetLogz();
1468  HistogramAttr2D(h2background, "background_{"+analysis2+"}", "background_{"+analysis1+"}", " ", -50, 50, -50, 50);
1469  h2background->Draw("colz");
1470  gPad->Print("Images/background2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1471 
1472  //chisq plots
1473  new TCanvas;
1474  gPad->SetGrid();
1475  HistogramAttr1D(hchisq, "#Deltachisq/chisq", "# of channels", 0.5*hchisq->GetXaxis()->GetBinCenter(1), 0.5*hchisq->GetXaxis()->GetBinCenter(hchisq->GetXaxis()->GetNbins()), 0., 4200., kBlack);
1476  hchisq->Draw();
1477  TLegend *legend = new TLegend(0.00, 0.60, 0.47, 0.93);
1478  legend->SetBorderSize(0);
1479  legend->SetFillStyle(0);
1480  legend->SetTextSize(0.05);
1481  legend->AddEntry((TObject*)0, det+" "+mcdata," ");
1482  legend->AddEntry(hchisq,"#splitline{"+analysis1+"}{vs "+analysis2+"}","P");
1483  legend->Draw();
1484  gPad->Print("Images/Delta_chisq_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1485 
1486  new TCanvas;
1487  gPad->SetGrid();
1488  gPad->SetLogz();
1489  TString chisq_tstring = "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1.Data+"}";
1490  char * chisq_char = chisq_tstring.Data();
1491  HistogramAttr2D(h2chisq, "#chi^{2}_{"+analysis2+"}", "#chi^{2}_{"+analysis1+"}", " ", 0, 0.5, 0, 0.5);
1492  h2chisq->Draw("colz");
1493  gPad->Print("Images/chisq2D_"+analysis1+"_vs_"+analysis2+"_"+det+"_"+mcdata+".pdf");
1494 
1495  in_1.close();
1496  in_2.close();
1497  out.close();
1498  fclose(fout);
1499 }//end of calibrationconstants()
1500 //---------------------------------------------------------------------------------------------------
1501 
1502 void zerosstudy(TFile * f0, TFile * f1, TFile * f2, TFile * f3, TFile * f4, TFile * f5, TFile * f6, TFile * f7, TFile * f8)
1503 {
1504  const int nbrightnesses = 9;
1505  //int nviews = 2;
1506  TFile * file[nbrightnesses] = {f0, f1, f2, f3, f4, f5, f6, f7, f8};
1507  TLegend * lzeros = new TLegend(0.5, 0.1, 0.88, 0.5);
1508  lzeros->SetBorderSize(0);
1509  lzeros->SetFillStyle(0);
1510  lzeros->SetTextSize(0.05);
1511 
1512  TH1F * hZerosVsWX[nbrightnesses];
1513  TH1F * hZerosVsWY[nbrightnesses];
1514  TH1F * hNCalibVsWX[nbrightnesses];
1515  TH1F * hNCalibVsWY[nbrightnesses];
1516  TH1F * hZerosVsCellX[nbrightnesses];
1517  TH1F * hZerosVsCellY[nbrightnesses];
1518  TH1F * hNCalibVsCellX[nbrightnesses];
1519  TH1F * hNCalibVsCellY[nbrightnesses];
1520 
1521  TCanvas * cZerosVsWX = new TCanvas;
1522  TCanvas * cZerosVsWY = new TCanvas;
1523  TCanvas * cCalibVsWX = new TCanvas;
1524  TCanvas * cCalibVsWY = new TCanvas;
1525  TCanvas * cZerosVsCellX = new TCanvas;
1526  TCanvas * cZerosVsCellY = new TCanvas;
1527  TCanvas * cCalibVsCellX = new TCanvas;
1528  TCanvas * cCalibVsCellY = new TCanvas;
1529 
1530  for(int fb = 0; fb < nbrightnesses; ++fb){
1531  TString fbstr;
1532  fbstr.Form("fb%i", fb);
1533  TString FBStr;
1534  FBStr.Form("FB %i", fb);
1535 
1536  std::cout << "Looking for nZerosVsW_xview_" << fbstr << std::endl;
1537  hZerosVsWX[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsW_xview_"+fbstr);
1538  hZerosVsWY[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsW_yview_"+fbstr);
1539  //hCalibVsWX[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsW_xview_"+fbstr);
1540  //hCalibVsWY[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsW_yview_"+fbstr);
1541  hZerosVsCellX[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsCell_xview_"+fbstr);
1542  hZerosVsCellY[fb] = (TH1F*)file[fb]->Get("thresh/nZerosVsCell_yview_"+fbstr);
1543  //hCalibVsCellX[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsCell_xview_"+fbstr);
1544  //hCalibVsCellY[fb] = (TH1F*)file[fb]->Get("thresh/nCalibVsCell_yview_"+fbstr);
1545 
1546  hZerosVsWX[fb]->SetTitle("XView");
1547  hZerosVsWY[fb]->SetTitle("YView");
1548 
1549  hZerosVsWX[fb]->SetLineColor(fbcolors[fb]);
1550  hZerosVsWY[fb]->SetLineColor(fbcolors[fb]);
1551  //hCalibVsWX[fb]->SetLineColor(fbcolors[fb]);
1552  //hCalibVsWY[fb]->SetLineColor(fbcolors[fb]);
1553  hZerosVsCellX[fb]->SetLineColor(fbcolors[fb]);
1554  hZerosVsCellY[fb]->SetLineColor(fbcolors[fb]);
1555  //hCalibVsCellX[fb]->SetLineColor(fbcolors[fb]);
1556  //hCalibVsCellY[fb]->SetLineColor(fbcolors[fb]);
1557 
1558  lzeros->AddEntry(hZerosVsWX[fb], FBStr.Data(), "l");
1559 
1560  cZerosVsWX->cd();
1561  hZerosVsWX[fb]->Draw("same");
1562  cZerosVsWY->cd();
1563  hZerosVsWY[fb]->Draw("same");
1564  cCalibVsWX->cd();
1565  //hCalibVsWX[fb]->Draw("same");
1566  cCalibVsWY->cd();
1567  //hCalibVsWY[fb]->Draw("same");
1568  cZerosVsCellX->cd();
1569  hZerosVsCellX[fb]->Draw("same");
1570  cZerosVsCellY->cd();
1571  hZerosVsCellY[fb]->Draw("same");
1572  cCalibVsCellX->cd();
1573  //hCalibVsCellX[fb]->Draw("same");
1574  cCalibVsCellY->cd();
1575  //hCalibVsCellY[fb]->Draw("same");
1576 
1577  }
1578 
1579  cZerosVsWX->cd();
1580  lzeros->Draw();
1581  cZerosVsWX->SaveAs("zeros_W_X.pdf");
1582 
1583  cZerosVsWY->cd();
1584  lzeros->Draw();
1585  cZerosVsWY->SaveAs("zeros_W_Y.pdf");
1586 
1587  cCalibVsWX->cd();
1588  lzeros->Draw();
1589  cCalibVsWX->SaveAs("calib_W_X.pdf");
1590 
1591  cCalibVsWY->cd();
1592  lzeros->Draw();
1593  cCalibVsWY->SaveAs("calib_W_Y.pdf");
1594 
1595  cZerosVsWX->cd();
1596  lzeros->Draw();
1597  cZerosVsWX->SaveAs("zeros_W_X.pdf");
1598 
1599  cZerosVsWY->cd();
1600  lzeros->Draw();
1601  cZerosVsWY->SaveAs("zeros_W_Y.pdf");
1602 
1603  cCalibVsWX->cd();
1604  lzeros->Draw();
1605  cCalibVsWX->SaveAs("calib_W_X.pdf");
1606 
1607  cCalibVsWY->cd();
1608  lzeros->Draw();
1609  cCalibVsWY->SaveAs("calib_W_Y.pdf");
1610 
1611 }
1612 
1613 // LocalWords: ProfileX
Float_t f4
TH1D * hist2
Definition: plotHisto.C:12
void thresholdshadowing(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString var)
void zerosstudy(TFile *f0, TFile *f1, TFile *f2, TFile *f3, TFile *f4, TFile *f5, TFile *f6, TFile *f7, TFile *f8)
void corrected_PE_mean(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell)
void calibrationconstants(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
Definition: AddMC.C:8
Float_t f2
static constexpr Double_t c_2
Definition: Munits.h:245
const XML_Char const XML_Char * data
Definition: expat.h:268
c2
Definition: demo5.py:33
fclose(fg1)
Float_t f3
std::map< ToFCounter, std::vector< unsigned int > > channels
correl_xv GetXaxis() -> SetDecimals()
TH1D * hist1
Definition: plotHisto.C:9
void uncorrected_PE_tricellvstrajectory(char *c1, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
void Ratiohist(TH1 *h1, TH1 *h2, TH1 *hratio, Color_t color, Style_t mstyle, TString labelx, TString labely, Double_t lowy, Double_t highy)
Float_t f1
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
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)
void corrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString plane, TString cell)
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)
uncorrected_PE(char *c1, char *c2, TString det, TString mcdata, TString analysis1, TString analysis2, TString hist1, TString hist2)
TFile * file
Definition: cellShifts.C:17
int ncells
Definition: geom.C:124
c1
Definition: demo5.py:24
HistogramAttr1D(TH1 *h, char *title_x, char *title_y, Double_t binlowx, Double_t binhighx, Double_t binlowy, Double_t binhighy, Color_t lcolor)
Float_t w
Definition: plot.C:20