ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CMDistortionAnalysis.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CMDistortionAnalysis.C
1 //step 3
2 #include <iostream>
3 #include <cmath>
4 #include <vector>
5 #include "TMath.h"
6 #include "TVector3.h"
7 #include "TTree.h"
8 
9 using namespace std;
10 
11 class Shifter {
12 public:
13 Shifter(TString sourcefilename);
14  TFile *forward, *average;
15  TH3F *hX, *hY, *hZ, *hR, *hPhi, *hXave, *hYave, *hZave, *hRave, *hPhiave;
16 };
17 
18 Shifter::Shifter(TString sourcefilename){
19  //single event distortion file
20  forward=TFile::Open(sourcefilename,"READ");
21 
22  hX=(TH3F*)forward->Get("hIntDistortionPosX");
23  hY=(TH3F*)forward->Get("hIntDistortionPosY");
24  hZ=(TH3F*)forward->Get("hIntDistortionPosZ");
25 
26  hR=(TH3F*)forward->Get("hIntDistortionPosR");
27  hPhi=(TH3F*)forward->Get("hIntDistortionPosP");
28 
29  //average distortion file
30  average=TFile::Open("/sphenix/user/rcorliss/distortion_maps/2021.04/apr07.average.real_B1.4_E-400.0.ross_phi1_sphenix_phislice_lookup_r26xp40xz40.distortion_map.hist.root","READ");
31 
32  hXave=(TH3F*)average->Get("hIntDistortionPosX");
33  hYave=(TH3F*)average->Get("hIntDistortionPosY");
34  hZave=(TH3F*)average->Get("hIntDistortionPosZ");
35 
36  hRave=(TH3F*)average->Get("hIntDistortionPosR");
37  hPhiave=(TH3F*)average->Get("hIntDistortionPosP");
38 
39  //subtract average from total distortions to study fluctuations
40  hX->Add(hXave,-1);
41  hY->Add(hYave,-1);
42  hZ->Add(hZave,-1);
43 
44  hR->Add(hRave,-1);
45  hPhi->Add(hPhiave,-1);
46 }
47 
48 int CMDistortionAnalysis(int nMaxEvents = -1) {
49  Shifter *shifter;
50  int nbins = 35;
51  double x, y, z;
52  double low = -80.0;
53  double high = 80.0;
54  double deltaX, deltaY, deltaZ, deltaR, deltaPhi;
55  int nEvents;
56 
57  TCanvas *canvas=new TCanvas("canvas","CMDistortionAnalysisPhiR",2000,3000);
58 
59  int nsumbins = 20;
60  int minsum = -10;
61  int maxsum = 10;
62 
63  //set up summary plots
64  TH1F *hDifferenceMeanR = new TH1F("hDifferenceMeanR", "Average Difference between R Model and True of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
65  TH1F *hDifferenceStdDevR = new TH1F("hDifferenceStdDevR", "Std Dev of Difference between R Model and True of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
66 
67  TH1F *hTrueMeanR = new TH1F("hTrueMeanR", "Mean True R Distortion Model of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
68  TH1F *hTrueStdDevR = new TH1F("hTrueStdDevR", "Std Dev of True R Distortion Model of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
69 
70  TH1F *hDifferenceMeanPhi = new TH1F("hDifferenceMeanPhi", "Average Difference between Phi Model and True of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
71  TH1F *hDifferenceStdDevPhi = new TH1F("hDifferenceStdDevPhi", "Std Dev of Difference between Phi Model and True of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
72 
73  TH1F *hTrueMeanPhi = new TH1F("hTrueMeanPhi", "Mean True Phi Distortion Model of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
74  TH1F *hTrueStdDevPhi = new TH1F("hTrueStdDevPhi", "Std Dev of True Phi Distortion Model of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
75 
76  const char * inputpattern="/sphenix/user/rcorliss/distortion_maps/2021.04/*h_Charge_*.root"; //updated
77 
78  //find all files that match the input string (includes wildcards)
79  TFileCollection *filelist=new TFileCollection();
80  filelist->Add(inputpattern);
81  TString sourcefilename;
82  /*
83  //how many events
84  if (nMaxEvents<0){
85  nEvents=filelist->GetNFiles();
86  } else if(nMaxEvents<filelist->GetNFiles()){
87  nEvents=nMaxEvents;
88  } else {
89  nEvents= filelist->GetNFiles();
90  }
91  */
92  nEvents = 2;
93 
94  for (int ifile=0;ifile < nEvents;ifile++){
95  //for each file, find all histograms in that file.
96  sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
97 
98  //create shifter
99  shifter = new Shifter(sourcefilename);
100 
101  TFile *plots;
102 
103  plots=TFile::Open(Form("CMModels_Event%d.root",ifile),"READ");
104 
105  TH3F *hCartCMModel[3];
106  hCartCMModel[0]=(TH3F*)plots->Get("hCMModelX");
107  hCartCMModel[1]=(TH3F*)plots->Get("hCMModelY");
108  hCartCMModel[2]=(TH3F*)plots->Get("hCMModelZ");
109 
110  TH3F *hCylCMModel[2];
111  hCylCMModel[0]=(TH3F*)plots->Get("hCMModelRCart");
112  hCylCMModel[1]=(TH3F*)plots->Get("hCMModelPhiCart");
113 
114  //phi,r binning
115  TH3F *hCartCMModelPhiR[3];
116  hCartCMModelPhiR[0]=(TH3F*)plots->Get("hCMModelX_PhiR");
117  hCartCMModelPhiR[1]=(TH3F*)plots->Get("hCMModelY_PhiR");
118  hCartCMModelPhiR[2]=(TH3F*)plots->Get("hCMModelZ_PhiR");
119 
120  TH3F *hCylCMModelPhiR[2];
121  hCylCMModelPhiR[0]=(TH3F*)plots->Get("hCMModelR_PhiR");
122  hCylCMModelPhiR[1]=(TH3F*)plots->Get("hCMModelPhi_PhiR");
123 
124  //for forward only
125 
126  //same range and bins for each coordinate, binned in cm
127  //hardcoded numbers from average distortion file's hIntDistortionPosX
128  int nphi = 82;
129  int nr = 54;
130  int nz = 82;
131 
132  double minphi = -0.078539819;
133  double minr = 18.884615;
134  double minz = -1.3187500;
135 
136  double maxphi = 6.3617253;
137  double maxr = 79.115387;
138  double maxz = 106.81875;
139 
140  double rshiftcart, phishiftcart;
141 
142  int ndiff = 300;
143  int mindiff = -20;
144  int maxdiff = 20;
145 
146 
147  TH1F *hCartesianShiftDifference[3];
148  hCartesianShiftDifference[0] = new TH1F("hShiftDifferenceX", "Difference between CM Model X and True (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
149  hCartesianShiftDifference[1] = new TH1F("hShiftDifferenceY", "Difference between CM Model Y and True (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
150  hCartesianShiftDifference[2] = new TH1F("hShiftDifferenceZ", "Difference between CM Model Z and True (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
151 
152  //phi,r binning
153  TH1F *hCartesianShiftDifferencePhiR[3];
154  hCartesianShiftDifferencePhiR[0] = new TH1F("hShiftDifferenceX_PhiR", "Difference between CM Model X and True, Phi,R binning (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
155  hCartesianShiftDifferencePhiR[1] = new TH1F("hShiftDifferenceY_PhiR", "Difference between CM Model Y and True, Phi,R binning (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
156  hCartesianShiftDifferencePhiR[2] = new TH1F("hShiftDifferenceZ_PhiR", "Difference between CM Model Z and True, Phi,R binning (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
157 
158  TH1F *hCylindricalShiftDifference[2];
159  hCylindricalShiftDifference[0] = new TH1F("hShiftDifferenceRCart", "Difference between CM Model R from Cartesian and True (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
160  hCylindricalShiftDifference[1] = new TH1F("hShiftDifferencePhiCart", "Difference between CM Model Phi from Cartesian and True (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
161 
162  //phi,r binning
163  TH1F *hCylindricalShiftDifferencePhiR[2];
164  hCylindricalShiftDifferencePhiR[0] = new TH1F("hShiftDifferenceR_PhiR", "Difference between CM Model R and True, Phi,R binning (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
165  hCylindricalShiftDifferencePhiR[1] = new TH1F("hShiftDifferencePhi_PhiR", "Difference between CM Model Phi and True, Phi,R binning (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
166 
167  TH1F *hRShiftTrue = new TH1F("hRShiftTrue", "True R Distortion Model (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
168  TH1F *hPhiShiftTrue = new TH1F("hPhiShiftTrue", "True Phi Distortion Model (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
169 
170  TH2F *hCartesianDiff[6];
171  hCartesianDiff[0] = new TH2F("hDiffXYX", "Difference in XY for CM Model X; x (cm); y (cm)",nbins,low,high,nbins,low,high);
172  hCartesianDiff[1] = new TH2F("hDiffRZX", "Difference in RZ for CM Model X; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
173  hCartesianDiff[2] = new TH2F("hDiffXYY", "Difference in XY for CM Model Y; x (cm); y (cm)",nbins,low,high,nbins,low,high);
174  hCartesianDiff[3] = new TH2F("hDiffRZY", "Difference in RZ for CM Model Y; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
175  hCartesianDiff[4] = new TH2F("hDiffXYZ", "Difference in XY for CM Model Z; x (cm); y (cm)",nbins,low,high,nbins,low,high);
176  hCartesianDiff[5] = new TH2F("hDiffRZZ", "Difference in RZ for CM Model Z; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
177 
178  //phi,r binning
179  TH2F *hCartesianDiffPhiR[6];
180  hCartesianDiffPhiR[0] = new TH2F("hDiffXYX_PhiR", "Difference in PhiR for CM Model X, Phi,R binning; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
181  hCartesianDiffPhiR[1] = new TH2F("hDiffRZX_PhiR", "Difference in RZ for CM Model X, Phi,R binning; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
182  hCartesianDiffPhiR[2] = new TH2F("hDiffXYY_PhiR", "Difference in PhiR for CM Model Y, Phi,R binning; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
183  hCartesianDiffPhiR[3] = new TH2F("hDiffRZY_PhiR", "Difference in RZ for CM Model Y, Phi,R binning; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
184  hCartesianDiffPhiR[4] = new TH2F("hDiffXYZ_PhiR", "Difference in PhiR for CM Model Z, Phi,R binning; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
185  hCartesianDiffPhiR[5] = new TH2F("hDiffRZZ_PhiR", "Difference in RZ for CM Model Z, Phi,R binning; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
186 
187  TH2F *hCylindricalDiff[4];
188  hCylindricalDiff[0] = new TH2F("hDiffXYRCart", "Difference in XY for CM Model R from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
189  hCylindricalDiff[1] = new TH2F("hDiffRZRCart", "Difference in RZ for CM Model R from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
190  hCylindricalDiff[2] = new TH2F("hDiffXYPhiCart", "Difference in XY for CM Model Phi from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
191  hCylindricalDiff[3] = new TH2F("hDiffRZPhiCart", "Difference in RZ for CM Model Phi from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
192 
193  //phi,r binning
194  TH2F *hCylindricalDiffPhiR[4];
195  hCylindricalDiffPhiR[0] = new TH2F("hDiffXYR_PhiR", "Difference in PhiR for CM Model R, Phi,R binning; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
196  hCylindricalDiffPhiR[1] = new TH2F("hDiffRZR_PhiR", "Difference in RZ for CM Model R, Phi,R binning; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
197  hCylindricalDiffPhiR[2] = new TH2F("hDiffXYPhi_PhiR", "Difference in PhiR for CM Model Phi, Phi,R binning; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
198  hCylindricalDiffPhiR[3] = new TH2F("hDiffRZPhi_PhiR", "Difference in RZ for CM Model Phi, Phi,R binning; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
199 
200  TH2F *hCartesianAveDiff[6];
201  hCartesianAveDiff[0] = new TH2F("hAveDiffXYX", "X Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
202  hCartesianAveDiff[1] = new TH2F("hAveDiffRZX", "X Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
203  hCartesianAveDiff[2] = new TH2F("hAveDiffXYY", "Y Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
204  hCartesianAveDiff[3] = new TH2F("hAveDiffRZY", "Y Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
205  hCartesianAveDiff[4] = new TH2F("hAveDiffXYZ", "Z Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
206  hCartesianAveDiff[5] = new TH2F("hAveDiffRZZ", "Z Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
207 
208  //phi,r binning
209  TH2F *hCartesianAveDiffPhiR[6];
210  hCartesianAveDiffPhiR[0] = new TH2F("hAveDiffXYX_PhiR", "X Model - Truth Averaged Over z, Phi,R binning (#mum); phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
211  hCartesianAveDiffPhiR[1] = new TH2F("hAveDiffRZX_PhiR", "X Model - Truth Averaged Over phi, Phi,R binning (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
212  hCartesianAveDiffPhiR[2] = new TH2F("hAveDiffXYY_PhiR", "Y Model - Truth Averaged Over z, Phi,R binning (#mum); phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
213  hCartesianAveDiffPhiR[3] = new TH2F("hAveDiffRZY_PhiR", "Y Model - Truth Averaged Over phi, Phi,R binning (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
214  hCartesianAveDiffPhiR[4] = new TH2F("hAveDiffXYZ_PhiR", "Z Model - Truth Averaged Over z, Phi,R binning (#mum); phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
215  hCartesianAveDiffPhiR[5] = new TH2F("hAveDiffRZZ_PhiR", "Z Model - Truth Averaged Over phi, Phi,R binning (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
216 
217  TH2F *hCylindricalAveDiff[4];
218  hCylindricalAveDiff[0] = new TH2F("hAveDiffXYRCart", "R Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
219  hCylindricalAveDiff[1] = new TH2F("hAveDiffRZRCart", "R Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
220  hCylindricalAveDiff[2] = new TH2F("hAveDiffXYPhiCart", "Phi Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
221  hCylindricalAveDiff[3] = new TH2F("hAveDiffRZPhiCart", "Phi Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
222 
223  //phi,r binning
224  TH2F *hCylindricalAveDiffPhiR[4];
225  hCylindricalAveDiffPhiR[0] = new TH2F("hAveDiffXYR_PhiR", "R Model - Truth Averaged Over z, Phi,R binning (#mum); phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
226  hCylindricalAveDiffPhiR[1] = new TH2F("hAveDiffRZR_PhiR", "R Model - Truth Averaged Over phi, Phi,R binning (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
227  hCylindricalAveDiffPhiR[2] = new TH2F("hAveDiffXYPhi_PhiR", "Phi Model - Truth Averaged Over z, Phi,R binning (#mum); phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
228  hCylindricalAveDiffPhiR[3] = new TH2F("hAveDiffRZPhi_PhiR", "Phi Model - Truth Averaged Over phi, Phi,R binning (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
229 
230  TH2F *hSamplePerBinXY = new TH2F("hSamplePerBinXY", "Filling each xy bin; x (cm); y (cm)",nbins,low,high,nbins,low,high);
231  TH2F *hSamplePerBinRZ = new TH2F("hSamplePerBinRZ", "Filling each rz bin; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
232 
233  //phi,r binning
234  TH2F *hSamplePerBinPhiR = new TH2F("hSamplePerBinPhiR", "Filling each PhiR bin; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
235 
236  TH2F *hCompareRTrue = new TH2F("hCompareRTrue", "Compare Difference from R Model and True (R > 30, 10 < z < 90); reco shift (#mum); true shift (#mum)",nbins,-550,550,nbins,-550,550);
237  TH2F *hComparePhiTrue = new TH2F("hComparePhiTrue", "Compare Difference from Phi Model and True (R > 30, 10 < z < 90); reco shift (#mum); true shift (#mum)",nbins,-550,550,nbins,-550,550);
238 
239  TH2F *hRDiffvR = new TH2F("hRDiffvR", "Difference between R Model and True vs. r (R > 30, 10 < z < 90); r (cm); shift difference (#mum)",nr,minr,maxr,ndiff,mindiff,maxdiff);
240  TH2F *hRDiffvZ = new TH2F("hRDiffvZ", "Difference between R Model and True vs. z (R > 30); z (cm); shift difference (#mum)",nz,minz,maxz,ndiff,mindiff,maxdiff);
241  TH2F *hRDiffvPhi = new TH2F("hRDiffvPhi", "Difference between R Model and True vs. phi (R > 30, 10 < z < 90); phi (rad); shift difference (#mum)",nphi,minphi,maxphi,ndiff,mindiff,maxdiff);
242 
243  TH2F *hPhiDiffvR = new TH2F("hPhiDiffvR", "Difference between Phi Model and True vs. r (R > 30, 10 < z < 90); r (cm); shift difference (#mum)",nr,minr,maxr,ndiff,mindiff,maxdiff);
244  TH2F *hPhiDiffvZ = new TH2F("hPhiDiffvZ", "Difference between Phi Model and True vs. z (R > 30); z (cm); shift difference (#mum)",nz,minz,maxz,ndiff,mindiff,maxdiff);
245  TH2F *hPhiDiffvPhi = new TH2F("hPhiDiffvPhi", "Difference between Phi Model and True vs. phi (R > 30, 10 < z < 90); phi (rad); shift difference (#mum)",nphi,minphi,maxphi,ndiff,mindiff,maxdiff);
246 
247  //phi,r binning
248  TH2F *hCompareRTrue_PhiR = new TH2F("hCompareRTrue_PhiR", "Compare Difference from R Model and True, Phi,R binning (R > 30, 10 < z < 90); reco shift (#mum); true shift (#mum)",nbins,-550,550,nbins,-550,550);
249  TH2F *hComparePhiTrue_PhiR = new TH2F("hComparePhiTrue_PhiR", "Compare Difference from Phi Model and True, Phi,R binning (R > 30, 10 < z < 90); reco shift (#mum); true shift (#mum)",nbins,-550,550,nbins,-550,550);
250 
251  TH2F *hRDiffvR_PhiR = new TH2F("hRDiffvR_PhiR", "Difference between R Model and True vs. r, Phi,R binning (R > 30, 10 < z < 90); r (cm); shift difference (#mum)",nr,minr,maxr,ndiff,mindiff,maxdiff);
252  TH2F *hRDiffvZ_PhiR = new TH2F("hRDiffvZ_PhiR", "Difference between R Model and True vs. z, Phi,R binning (R > 30); z (cm); shift difference (#mum)",nz,minz,maxz,ndiff,mindiff,maxdiff);
253  TH2F *hRDiffvPhi_PhiR = new TH2F("hRDiffvPhi_PhiR", "Difference between R Model and True vs. phi, Phi,R binning (R > 30, 10 < z < 90); phi (rad); shift difference (#mum)",nphi,minphi,maxphi,ndiff,mindiff,maxdiff);
254 
255  TH2F *hPhiDiffvR_PhiR = new TH2F("hPhiDiffvR_PhiR", "Difference between Phi Model and True vs. r, Phi,R binning (R > 30, 10 < z < 90); r (cm); shift difference (#mum)",nr,minr,maxr,ndiff,mindiff,maxdiff);
256  TH2F *hPhiDiffvZ_PhiR = new TH2F("hPhiDiffvZ_PhiR", "Difference between Phi Model and True vs. z, Phi,R binning (R > 30); z (cm); shift difference (#mum)",nz,minz,maxz,ndiff,mindiff,maxdiff);
257  TH2F *hPhiDiffvPhi_PhiR = new TH2F("hPhiDiffvPhi_PhiR", "Difference between Phi Model and True vs. phi, Phi,R binning (R > 30, 10 < z < 90); phi (rad); shift difference (#mum)",nphi,minphi,maxphi,ndiff,mindiff,maxdiff);
258 
259  for(int i = 1; i < nphi - 1; i++){
260  double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
261  for(int j = 1; j < nr - 1; j++){
262  double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
263  for(int k = 1; k < nz - 1; k++){
264  double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
265 
266  double shiftrecoCart[3];
267  double shifttrueCart[3];
268  double differenceCart[3];
269 
270  double shiftrecoCyl[2];
271  double shifttrueCyl[2];
272  double differenceCyl[2];
273 
274  double differenceR, differencePhi;
275 
276  int bin = hCartCMModel[0]->FindBin(phi,r,z); //same for all
277 
278  //phi,r binning
279  double shiftrecoCartPhiR[3];
280  double differenceCartPhiR[3];
281  double shiftrecoCylPhiR[2];
282  double differenceCylPhiR[2];
283  double differenceR_PhiR, differencePhi_PhiR;
284  int binPhiR = hCartCMModelPhiR[0]->FindBin(phi,r,z);
285 
286  if((r > 30.0) && (r < 76.0)){
287  //x y and z
288  shifttrueCart[0] = (shifter->hX->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
289  shifttrueCart[1] = (shifter->hY->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
290  shifttrueCart[2] = (shifter->hZ->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
291  for(int l = 0; l < 3; l ++){
292  shiftrecoCart[l] = (hCartCMModel[l]->GetBinContent(bin))*(1e4);
293 
294  differenceCart[l] = shiftrecoCart[l] - shifttrueCart[l];
295 
296  hCartesianShiftDifference[l]->Fill(differenceCart[l]);
297  }
298 
299  //r from cart
300  shiftrecoCyl[0] = (hCylCMModel[0]->GetBinContent(bin))*(1e4);
301  shifttrueCyl[0] = (shifter->hR->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
302  differenceCyl[0] = shiftrecoCyl[0] - shifttrueCyl[0];
303  hCylindricalShiftDifference[0]->Fill(differenceCyl[0]);
304 
305  //phi from cart
306  shiftrecoCyl[1] = r*(1e4)*(hCylCMModel[1]->GetBinContent(bin));
307  shifttrueCyl[1] = (shifter->hPhi->Interpolate(phi,r,z))*(1e4);
308  differenceCyl[1] = (shiftrecoCyl[1] - shifttrueCyl[1]);
309  hCylindricalShiftDifference[1]->Fill(differenceCyl[1]);
310 
311  hRShiftTrue->Fill(shifttrueCyl[0]);
312  hPhiShiftTrue->Fill(shifttrueCyl[1]);
313 
314  double x = r*cos(phi);
315  double y = r*sin(phi);
316 
317  //x
318  hCartesianDiff[0]->Fill(x,y, differenceCart[0]);
319  hCartesianDiff[1]->Fill(z,r, differenceCart[0]);
320  //y
321  hCartesianDiff[2]->Fill(x,y, differenceCart[1]);
322  hCartesianDiff[3]->Fill(z,r, differenceCart[1]);
323  //z
324  hCartesianDiff[4]->Fill(x,y, differenceCart[2]);
325  hCartesianDiff[5]->Fill(z,r, differenceCart[2]);
326 
327  //r cart
328  hCylindricalDiff[0]->Fill(x,y, differenceCyl[0]);
329  hCylindricalDiff[1]->Fill(z,r, differenceCyl[0]);
330  //phi cart
331  hCylindricalDiff[2]->Fill(x,y, differenceCyl[1]);
332  hCylindricalDiff[3]->Fill(z,r, differenceCyl[1]);
333 
334  hCompareRTrue->Fill(shiftrecoCyl[0],shifttrueCyl[0]);
335  hComparePhiTrue->Fill(shiftrecoCyl[1],shifttrueCyl[1]);
336 
337  hRDiffvR->Fill(r,differenceCyl[0],1);
338  hRDiffvPhi->Fill(phi,differenceCyl[0],1);
339  hRDiffvZ->Fill(z,differenceCyl[0],1);
340 
341  hPhiDiffvR->Fill(r,differenceCyl[1],1);
342  hPhiDiffvPhi->Fill(phi,differenceCyl[1],1);
343  hPhiDiffvZ->Fill(z,differenceCyl[1],1);
344 
345  hSamplePerBinXY->Fill(x,y,1);
346 
347  hSamplePerBinRZ->Fill(z,r,1);
348 
349  //phi,r binning
350  //x y and z
351  for(int l = 0; l < 3; l ++){
352  shiftrecoCartPhiR[l] = (hCartCMModelPhiR[l]->GetBinContent(binPhiR))*(1e4);
353 
354  differenceCartPhiR[l] = shiftrecoCartPhiR[l] - shifttrueCart[l];
355 
356  hCartesianShiftDifferencePhiR[l]->Fill(differenceCartPhiR[l]);
357  }
358 
359  //r
360  shiftrecoCylPhiR[0] = (hCylCMModelPhiR[0]->GetBinContent(binPhiR))*(1e4);
361  differenceCylPhiR[0] = shiftrecoCylPhiR[0] - shifttrueCyl[0];
362  hCylindricalShiftDifferencePhiR[0]->Fill(differenceCylPhiR[0]);
363 
364  //phi
365  shiftrecoCylPhiR[1] = r*(1e4)*(hCylCMModelPhiR[1]->GetBinContent(binPhiR));
366  differenceCylPhiR[1] = (shiftrecoCylPhiR[1] - shifttrueCyl[1]);
367  hCylindricalShiftDifferencePhiR[1]->Fill(differenceCylPhiR[1]);
368 
369  //x
370  hCartesianDiffPhiR[0]->Fill(phi,r, differenceCartPhiR[0]);
371  hCartesianDiffPhiR[1]->Fill(z,r, differenceCartPhiR[0]);
372  //y
373  hCartesianDiffPhiR[2]->Fill(phi,r, differenceCartPhiR[1]);
374  hCartesianDiffPhiR[3]->Fill(z,r, differenceCartPhiR[1]);
375  //z
376  hCartesianDiffPhiR[4]->Fill(phi,r, differenceCartPhiR[2]);
377  hCartesianDiffPhiR[5]->Fill(z,r, differenceCartPhiR[2]);
378 
379  //r
380  hCylindricalDiffPhiR[0]->Fill(phi,r, differenceCylPhiR[0]);
381  hCylindricalDiffPhiR[1]->Fill(z,r, differenceCylPhiR[0]);
382  //phi
383  hCylindricalDiffPhiR[2]->Fill(phi,r, differenceCylPhiR[1]);
384  hCylindricalDiffPhiR[3]->Fill(z,r, differenceCylPhiR[1]);
385 
386  hCompareRTrue_PhiR->Fill(shiftrecoCylPhiR[0],shifttrueCyl[0]);
387  hComparePhiTrue_PhiR->Fill(shiftrecoCylPhiR[1],shifttrueCyl[1]);
388 
389  hRDiffvR_PhiR->Fill(r,differenceCylPhiR[0],1);
390  hRDiffvPhi_PhiR->Fill(phi,differenceCylPhiR[0],1);
391  hRDiffvZ_PhiR->Fill(z,differenceCylPhiR[0],1);
392 
393  hPhiDiffvR_PhiR->Fill(r,differenceCylPhiR[1],1);
394  hPhiDiffvPhi_PhiR->Fill(phi,differenceCylPhiR[1],1);
395  hPhiDiffvZ_PhiR->Fill(z,differenceCylPhiR[1],1);
396 
397  hSamplePerBinPhiR->Fill(phi,r,1);
398  }
399  }
400  }
401  }
402 
403  //average over z
404  for (int m = 0; m < 6; m = m+2){
405  hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinXY);
406  }
407  for (int m = 0; m < 4; m = m+2){
408  hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinXY);
409  }
410 
411  //average over phi
412  for (int m = 1; m < 6; m = m+2){
413  hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinRZ);
414  }
415  for (int m = 1; m < 4; m = m+2){
416  hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinRZ);
417  }
418 
419  //phi,r binning
420  //average over z
421  for (int m = 0; m < 6; m = m+2){
422  hCartesianAveDiffPhiR[m]->Divide(hCartesianDiffPhiR[m],hSamplePerBinPhiR);
423  }
424  for (int m = 0; m < 4; m = m+2){
425  hCylindricalAveDiffPhiR[m]->Divide(hCylindricalDiffPhiR[m],hSamplePerBinPhiR);
426  }
427 
428  //average over phi
429  for (int m = 1; m < 6; m = m+2){
430  hCartesianAveDiffPhiR[m]->Divide(hCartesianDiffPhiR[m],hSamplePerBinRZ);
431  }
432  for (int m = 1; m < 4; m = m+2){
433  hCylindricalAveDiffPhiR[m]->Divide(hCylindricalDiffPhiR[m],hSamplePerBinRZ);
434  }
435 
436  //summary plots
437  hDifferenceMeanR->Fill(hCylindricalShiftDifferencePhiR[0]->GetMean(1));
438  hDifferenceStdDevR->Fill(hCylindricalShiftDifferencePhiR[0]->GetStdDev(1));
439 
440  hTrueMeanR->Fill(hRShiftTrue->GetMean(1));
441  hTrueStdDevR->Fill(hRShiftTrue->GetStdDev(1));
442 
443  hDifferenceMeanPhi->Fill(hCylindricalShiftDifferencePhiR[1]->GetMean(1));
444  hDifferenceStdDevPhi->Fill(hCylindricalShiftDifferencePhiR[1]->GetStdDev(1));
445 
446  hTrueMeanPhi->Fill(hPhiShiftTrue->GetMean(1));
447  hTrueStdDevPhi->Fill(hPhiShiftTrue->GetStdDev(1));
448 
449  for (int m = 0; m < 6; m++){
450  hCartesianAveDiff[m]->SetStats(0);
451  hCartesianAveDiffPhiR[m]->SetStats(0);
452  }
453  for (int m = 0; m < 4; m++){
454  hCylindricalAveDiff[m]->SetStats(0);
455  hCylindricalAveDiffPhiR[m]->SetStats(0);
456  }
457 
458  hCompareRTrue->SetStats(0);
459  hComparePhiTrue->SetStats(0);
460 
461  hRDiffvR->SetStats(0);
462  hRDiffvZ->SetStats(0);
463  hRDiffvPhi->SetStats(0);
464 
465  hPhiDiffvR->SetStats(0);
466  hPhiDiffvZ->SetStats(0);
467  hPhiDiffvPhi->SetStats(0);
468 
469  hCompareRTrue_PhiR->SetStats(0);
470  hComparePhiTrue_PhiR->SetStats(0);
471 
472  hRDiffvR_PhiR->SetStats(0);
473  hRDiffvZ_PhiR->SetStats(0);
474  hRDiffvPhi_PhiR->SetStats(0);
475 
476  hPhiDiffvR_PhiR->SetStats(0);
477  hPhiDiffvZ_PhiR->SetStats(0);
478  hPhiDiffvPhi_PhiR->SetStats(0);
479 
480  TPad *c1=new TPad("c1","",0.0,0.8,1.0,0.93); //can i do an array of pads?
481  TPad *c2=new TPad("c2","",0.0,0.64,1.0,0.77);
482  TPad *c3=new TPad("c3","",0.0,0.48,1.0,0.61);
483  TPad *c4=new TPad("c4","",0.0,0.32,1.0,0.45);
484  TPad *c5=new TPad("c5","",0.0,0.16,1.0,0.29);
485  TPad *c6=new TPad("c6","",0.0,0.0,1.0,0.13);
486 
487  TPad *titlepad=new TPad("titlepad","",0.0,0.96,1.0,1.0);
488 
489  TPad *stitlepad1=new TPad("stitlepad1","",0.0,0.93,1.0,0.96);
490  TPad *stitlepad2=new TPad("stitlepad2","",0.0,0.77,1.0,0.8);
491  TPad *stitlepad3=new TPad("stitlepad3","",0.0,0.61,1.0,0.64);
492  TPad *stitlepad4=new TPad("stitlepad4","",0.0,0.45,1.0,0.48);
493  TPad *stitlepad5=new TPad("stitlepad5","",0.0,0.29,1.0,0.32);
494  TPad *stitlepad6=new TPad("stitlepad6","",0.0,0.13,1.0,0.16);
495 
496  TLatex * title = new TLatex(0.0,0.0,"");
497 
498  TLatex * stitle1 = new TLatex(0.0,0.0,""); //array?
499  TLatex * stitle2 = new TLatex(0.0,0.0,"");
500  TLatex * stitle3 = new TLatex(0.0,0.0,"");
501  TLatex * stitle4 = new TLatex(0.0,0.0,"");
502  TLatex * stitle5 = new TLatex(0.0,0.0,"");
503  TLatex * stitle6 = new TLatex(0.0,0.0,"");
504 
505  title->SetNDC();
506  stitle1->SetNDC();
507  stitle2->SetNDC();
508  stitle3->SetNDC();
509  stitle4->SetNDC();
510  stitle5->SetNDC();
511  stitle6->SetNDC();
512 
513  title->SetTextSize(0.32);
514  stitle1->SetTextSize(0.35);
515  stitle2->SetTextSize(0.35);
516  stitle3->SetTextSize(0.35);
517  stitle4->SetTextSize(0.35);
518  stitle5->SetTextSize(0.35);
519  stitle6->SetTextSize(0.35);
520 
521  canvas->cd();
522  c1->Draw();
523  stitlepad1->Draw();
524  c2->Draw();
525  stitlepad2->Draw();
526  c3->Draw();
527  stitlepad3->Draw();
528  c4->Draw();
529  stitlepad4->Draw();
530  c5->Draw();
531  stitlepad5->Draw();
532  c6->Draw();
533  stitlepad6->Draw();
534  titlepad->Draw();
535 
536  //x plots
537  c1->Divide(4,1);
538  c1->cd(1);
539  hCartesianAveDiffPhiR[0]->Draw("colz");
540  c1->cd(2);
541  hCartesianAveDiffPhiR[1]->Draw("colz");
542  c1->cd(3);
543  hCartesianShiftDifferencePhiR[0]->Draw();
544  //c1->cd(4)->Clear();
545  c1->cd(4);
546  //hCMmodelSliceRvTrue->Draw("colz");
547  hSamplePerBinRZ->Draw("colz");
548 
549  //y plots
550  c2->Divide(4,1);
551  c2->cd(1);
552  hCartesianAveDiffPhiR[2]->Draw("colz");
553  c2->cd(2);
554  hCartesianAveDiffPhiR[3]->Draw("colz");
555  c2->cd(3);
556  hCartesianShiftDifferencePhiR[1]->Draw();
557  //c2->cd(4)->Clear();
558  c2->cd(4);
559  //hStripesPerBin->Draw("colz");
560  hSamplePerBinPhiR->Draw("colz");
561 
562  //r cart
563  c3->Divide(4,1);
564  c3->cd(1);
565  hCylindricalAveDiffPhiR[0]->Draw("colz");
566  c3->cd(2);
567  hCylindricalAveDiffPhiR[1]->Draw("colz");
568  c3->cd(3);
569  hCylindricalShiftDifferencePhiR[0]->Draw();
570  c3->cd(4);
571  hRShiftTrue->Draw();
572 
573  //phi cart
574  c4->Divide(4,1);
575  c4->cd(1);
576  hCylindricalAveDiffPhiR[2]->Draw("colz");
577  c4->cd(2);
578  hCylindricalAveDiffPhiR[3]->Draw("colz");
579  c4->cd(3);
580  hCylindricalShiftDifferencePhiR[1]->Draw();
581  c4->cd(4);
582  hPhiShiftTrue->Draw();
583 
584  //r to true comparison
585  c5->Divide(4,1);
586  c5->cd(1);
587  hCompareRTrue_PhiR->Draw("colz");
588  c5->cd(2);
589  hRDiffvR_PhiR->Draw("colz");
590  c5->cd(3);
591  hRDiffvZ_PhiR->Draw("colz");
592  c5->cd(4);
593  hRDiffvPhi_PhiR->Draw("colz");
594 
595  //phi to true comparison
596  c6->Divide(4,1);
597  c6->cd(1);
598  hComparePhiTrue_PhiR->Draw("colz");
599  c6->cd(2);
600  hPhiDiffvR_PhiR->Draw("colz");
601  c6->cd(3);
602  hPhiDiffvZ_PhiR->Draw("colz");
603  c6->cd(4);
604  hPhiDiffvPhi_PhiR->Draw("colz");
605 
606  titlepad->cd();
607  titlepad->Clear();
608  title->DrawLatex(0.01,0.4,Form("Event %d; %s", ifile, sourcefilename.Data()));
609  title->Draw();
610 
611  stitlepad1->cd();
612  stitlepad1->Clear();
613  stitle1->DrawLatex(0.45,0.2,"X Model");
614  stitle1->Draw();
615 
616  stitlepad2->cd();
617  stitlepad2->Clear();
618  stitle2->DrawLatex(0.45,0.2,"Y Model");
619  stitle2->Draw();
620 
621  stitlepad3->cd();
622  stitlepad3->Clear();
623  stitle3->DrawLatex(0.45,0.2,"R Model");
624  stitle3->Draw();
625 
626  stitlepad4->cd();
627  stitlepad4->Clear();
628  stitle4->DrawLatex(0.45,0.2,"Phi Model");
629  stitle4->Draw();
630 
631  stitlepad5->cd();
632  stitlepad5->Clear();
633  stitle5->DrawLatex(0.4,0.2,"Comparing R Model to True");
634  stitle5->Draw();
635 
636  stitlepad6->cd();
637  stitlepad6->Clear();
638  stitle6->DrawLatex(0.4,0.2,"Comparing Phi Model to True");
639  stitle6->Draw();
640 
641  if(ifile == 0){
642  //if(ifile == 1){
643  canvas->Print("CMDistortionAnalysisPhiR.pdf(","pdf");
644  } else if((ifile == 1) || (ifile == nEvents - 1)){
645  canvas->Print("CMDistortionAnalysisPhiR.pdf","pdf");
646  }
647  }
648 
649  TCanvas *summary = new TCanvas("summary","ShiftPlotsSummary",2000,3000);
650 
651  TPad *sumtitlepad = new TPad("sumtitlepad","",0.0,0.96,1.0,1.0);
652  TPad *sumplots = new TPad("sumplotspad","",0.0,0.0,1.0,0.96);
653 
654  TLatex *sumtitle = new TLatex(0.0,0.0,"");
655 
656  sumtitle->SetNDC();
657  sumtitle->SetTextSize(0.4);
658 
659  summary->cd();
660  sumplots->Draw();
661  sumtitlepad->Draw();
662 
663  sumplots->Divide(4,6);
664  sumplots->cd(1);
665  hDifferenceMeanR->Draw();
666  sumplots->cd(2);
667  hDifferenceStdDevR->Draw();
668  sumplots->cd(3);
669  hTrueMeanR->Draw();
670  sumplots->cd(4);
671  hTrueStdDevR->Draw();
672  sumplots->cd(5);
673  hDifferenceMeanPhi->Draw();
674  sumplots->cd(6);
675  hDifferenceStdDevPhi->Draw();
676  sumplots->cd(7);
677  hTrueMeanPhi->Draw();
678  sumplots->cd(8);
679  hTrueStdDevPhi->Draw();
680  sumplots->cd(9);
681  sumplots->cd(10)->Clear();
682  sumplots->cd(11)->Clear();
683  sumplots->cd(12)->Clear();
684  sumplots->cd(13)->Clear();
685  sumplots->cd(14)->Clear();
686  sumplots->cd(15)->Clear();
687  sumplots->cd(16)->Clear();
688  sumplots->cd(17)->Clear();
689  sumplots->cd(18)->Clear();
690  sumplots->cd(19)->Clear();
691  sumplots->cd(20)->Clear();
692  sumplots->cd(21)->Clear();
693  sumplots->cd(22)->Clear();
694  sumplots->cd(23)->Clear();
695  sumplots->cd(24)->Clear();
696 
697  sumtitlepad->cd();
698  sumtitlepad->Clear();
699  sumtitle->DrawLatex(0.4,0.4,"Summary of Events");
700  summary->Print("CMDistortionAnalysisPhiR.pdf)","pdf");
701 
702  return 0;
703 }