Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CMDistortionAnalysisCart.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CMDistortionAnalysisCart.C
1 //step 3 in cart coords
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 CMDistortionAnalysisCart(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","CMDistortionAnalysisCart",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("CMModelsCart_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  TH1F *hCartesianShiftDifference[3];
147  hCartesianShiftDifference[0] = new TH1F("hShiftDifferenceX", "Difference between CM Model X and True (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
148  hCartesianShiftDifference[1] = new TH1F("hShiftDifferenceY", "Difference between CM Model Y and True (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
149  hCartesianShiftDifference[2] = new TH1F("hShiftDifferenceZ", "Difference between CM Model Z and True (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
150 
151  TH1F *hCylindricalShiftDifference[2];
152  hCylindricalShiftDifference[0] = new TH1F("hShiftDifferenceRCart", "Difference between CM Model R from Cartesian and True (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
153  hCylindricalShiftDifference[1] = new TH1F("hShiftDifferencePhiCart", "Difference between CM Model Phi from Cartesian and True (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
154 
155  TH1F *hRShiftTrue = new TH1F("hRShiftTrue", "True R Distortion Model (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
156  TH1F *hPhiShiftTrue = new TH1F("hPhiShiftTrue", "True Phi Distortion Model (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
157 
158  TH2F *hCartesianDiff[6];
159  hCartesianDiff[0] = new TH2F("hDiffXYX", "Difference in XY for CM Model X; x (cm); y (cm)",nbins,low,high,nbins,low,high);
160  hCartesianDiff[1] = new TH2F("hDiffRZX", "Difference in RZ for CM Model X; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
161  hCartesianDiff[2] = new TH2F("hDiffXYY", "Difference in XY for CM Model Y; x (cm); y (cm)",nbins,low,high,nbins,low,high);
162  hCartesianDiff[3] = new TH2F("hDiffRZY", "Difference in RZ for CM Model Y; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
163  hCartesianDiff[4] = new TH2F("hDiffXYZ", "Difference in XY for CM Model Z; x (cm); y (cm)",nbins,low,high,nbins,low,high);
164  hCartesianDiff[5] = new TH2F("hDiffRZZ", "Difference in RZ for CM Model Z; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
165 
166  TH2F *hCylindricalDiff[4];
167  hCylindricalDiff[0] = new TH2F("hDiffXYRCart", "Difference in XY for CM Model R from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
168  hCylindricalDiff[1] = new TH2F("hDiffRZRCart", "Difference in RZ for CM Model R from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
169  hCylindricalDiff[2] = new TH2F("hDiffXYPhiCart", "Difference in XY for CM Model Phi from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
170  hCylindricalDiff[3] = new TH2F("hDiffRZPhiCart", "Difference in RZ for CM Model Phi from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
171 
172  TH2F *hCartesianAveDiff[6];
173  hCartesianAveDiff[0] = new TH2F("hAveDiffXYX", "X Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
174  hCartesianAveDiff[1] = new TH2F("hAveDiffRZX", "X Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
175  hCartesianAveDiff[2] = new TH2F("hAveDiffXYY", "Y Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
176  hCartesianAveDiff[3] = new TH2F("hAveDiffRZY", "Y Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
177  hCartesianAveDiff[4] = new TH2F("hAveDiffXYZ", "Z Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
178  hCartesianAveDiff[5] = new TH2F("hAveDiffRZZ", "Z Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
179 
180  TH2F *hCylindricalAveDiff[4];
181  hCylindricalAveDiff[0] = new TH2F("hAveDiffXYRCart", "R Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
182  hCylindricalAveDiff[1] = new TH2F("hAveDiffRZRCart", "R Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
183  hCylindricalAveDiff[2] = new TH2F("hAveDiffXYPhiCart", "Phi Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
184  hCylindricalAveDiff[3] = new TH2F("hAveDiffRZPhiCart", "Phi Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
185 
186  TH2F *hSamplePerBinXY = new TH2F("hSamplePerBinXY", "Filling each xy bin; x (cm); y (cm)",nbins,low,high,nbins,low,high);
187  TH2F *hSamplePerBinRZ = new TH2F("hSamplePerBinRZ", "Filling each rz bin; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
188 
189  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);
190  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);
191 
192  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);
193  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);
194  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);
195 
196  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);
197  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);
198  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);
199 
200  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);
201  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);
202  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);
203 
204  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);
205  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);
206  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);
207 
208  for(int i = 1; i < nphi - 1; i++){
209  double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
210  for(int j = 1; j < nr - 1; j++){
211  double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
212  for(int k = 1; k < nz - 1; k++){
213  double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
214 
215  double shiftrecoCart[3];
216  double shifttrueCart[3];
217  double differenceCart[3];
218 
219  double shiftrecoCyl[2];
220  double shifttrueCyl[2];
221  double differenceCyl[2];
222 
223  double differenceR, differencePhi;
224 
225  int bin = hCartCMModel[0]->FindBin(phi,r,z); //same for all
226 
227  if((r > 30.0) && (r < 76.0)){
228  //x y and z
229  shifttrueCart[0] = (shifter->hX->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
230  shifttrueCart[1] = (shifter->hY->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
231  shifttrueCart[2] = (shifter->hZ->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
232  for(int l = 0; l < 3; l ++){
233  shiftrecoCart[l] = (hCartCMModel[l]->GetBinContent(bin))*(1e4);
234 
235  differenceCart[l] = shiftrecoCart[l] - shifttrueCart[l];
236 
237  hCartesianShiftDifference[l]->Fill(differenceCart[l]);
238  }
239 
240  //r from cart
241  shiftrecoCyl[0] = (hCylCMModel[0]->GetBinContent(bin))*(1e4);
242  shifttrueCyl[0] = (shifter->hR->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
243  differenceCyl[0] = shiftrecoCyl[0] - shifttrueCyl[0];
244  hCylindricalShiftDifference[0]->Fill(differenceCyl[0]);
245 
246  //phi from cart
247  shiftrecoCyl[1] = r*(1e4)*(hCylCMModel[1]->GetBinContent(bin));
248  shifttrueCyl[1] = (shifter->hPhi->Interpolate(phi,r,z))*(1e4);
249  differenceCyl[1] = (shiftrecoCyl[1] - shifttrueCyl[1]);
250  hCylindricalShiftDifference[1]->Fill(differenceCyl[1]);
251 
252  hRShiftTrue->Fill(shifttrueCyl[0]);
253  hPhiShiftTrue->Fill(shifttrueCyl[1]);
254 
255  double x = r*cos(phi);
256  double y = r*sin(phi);
257 
258  //x
259  hCartesianDiff[0]->Fill(x,y, differenceCart[0]);
260  hCartesianDiff[1]->Fill(z,r, differenceCart[0]);
261  //y
262  hCartesianDiff[2]->Fill(x,y, differenceCart[1]);
263  hCartesianDiff[3]->Fill(z,r, differenceCart[1]);
264  //z
265  hCartesianDiff[4]->Fill(x,y, differenceCart[2]);
266  hCartesianDiff[5]->Fill(z,r, differenceCart[2]);
267 
268  //r cart
269  hCylindricalDiff[0]->Fill(x,y, differenceCyl[0]);
270  hCylindricalDiff[1]->Fill(z,r, differenceCyl[0]);
271  //phi cart
272  hCylindricalDiff[2]->Fill(x,y, differenceCyl[1]);
273  hCylindricalDiff[3]->Fill(z,r, differenceCyl[1]);
274 
275  hCompareRTrue->Fill(shiftrecoCyl[0],shifttrueCyl[0]);
276  hComparePhiTrue->Fill(shiftrecoCyl[1],shifttrueCyl[1]);
277 
278  hRDiffvR->Fill(r,differenceCyl[0],1);
279  hRDiffvPhi->Fill(phi,differenceCyl[0],1);
280  hRDiffvZ->Fill(z,differenceCyl[0],1);
281 
282  hPhiDiffvR->Fill(r,differenceCyl[1],1);
283  hPhiDiffvPhi->Fill(phi,differenceCyl[1],1);
284  hPhiDiffvZ->Fill(z,differenceCyl[1],1);
285 
286  hSamplePerBinXY->Fill(x,y,1);
287 
288  hSamplePerBinRZ->Fill(z,r,1);
289  }
290  }
291  }
292  }
293 
294  //average over z
295  for (int m = 0; m < 6; m = m+2){
296  hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinXY);
297  }
298  for (int m = 0; m < 4; m = m+2){
299  hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinXY);
300  }
301 
302  //average over phi
303  for (int m = 1; m < 6; m = m+2){
304  hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinRZ);
305  }
306  for (int m = 1; m < 4; m = m+2){
307  hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinRZ);
308  }
309 
310  //summary plots
311  hDifferenceMeanR->Fill(hCylindricalShiftDifference[0]->GetMean(1));
312  hDifferenceStdDevR->Fill(hCylindricalShiftDifference[0]->GetStdDev(1));
313 
314  hTrueMeanR->Fill(hRShiftTrue->GetMean(1));
315  hTrueStdDevR->Fill(hRShiftTrue->GetStdDev(1));
316 
317  hDifferenceMeanPhi->Fill(hCylindricalShiftDifference[1]->GetMean(1));
318  hDifferenceStdDevPhi->Fill(hCylindricalShiftDifference[1]->GetStdDev(1));
319 
320  hTrueMeanPhi->Fill(hPhiShiftTrue->GetMean(1));
321  hTrueStdDevPhi->Fill(hPhiShiftTrue->GetStdDev(1));
322 
323  for (int m = 0; m < 6; m++){
324  hCartesianAveDiff[m]->SetStats(0);
325  }
326  for (int m = 0; m < 4; m++){
327  hCylindricalAveDiff[m]->SetStats(0);
328  }
329 
330  hCompareRTrue->SetStats(0);
331  hComparePhiTrue->SetStats(0);
332 
333  hRDiffvR->SetStats(0);
334  hRDiffvZ->SetStats(0);
335  hRDiffvPhi->SetStats(0);
336 
337  hPhiDiffvR->SetStats(0);
338  hPhiDiffvZ->SetStats(0);
339  hPhiDiffvPhi->SetStats(0);
340 
341  TPad *c1=new TPad("c1","",0.0,0.8,1.0,0.93); //can i do an array of pads?
342  TPad *c2=new TPad("c2","",0.0,0.64,1.0,0.77);
343  TPad *c3=new TPad("c3","",0.0,0.48,1.0,0.61);
344  TPad *c4=new TPad("c4","",0.0,0.32,1.0,0.45);
345  TPad *c5=new TPad("c5","",0.0,0.16,1.0,0.29);
346  TPad *c6=new TPad("c6","",0.0,0.0,1.0,0.13);
347 
348  TPad *titlepad=new TPad("titlepad","",0.0,0.96,1.0,1.0);
349 
350  TPad *stitlepad1=new TPad("stitlepad1","",0.0,0.93,1.0,0.96);
351  TPad *stitlepad2=new TPad("stitlepad2","",0.0,0.77,1.0,0.8);
352  TPad *stitlepad3=new TPad("stitlepad3","",0.0,0.61,1.0,0.64);
353  TPad *stitlepad4=new TPad("stitlepad4","",0.0,0.45,1.0,0.48);
354  TPad *stitlepad5=new TPad("stitlepad5","",0.0,0.29,1.0,0.32);
355  TPad *stitlepad6=new TPad("stitlepad6","",0.0,0.13,1.0,0.16);
356 
357  TLatex * title = new TLatex(0.0,0.0,"");
358 
359  TLatex * stitle1 = new TLatex(0.0,0.0,""); //array?
360  TLatex * stitle2 = new TLatex(0.0,0.0,"");
361  TLatex * stitle3 = new TLatex(0.0,0.0,"");
362  TLatex * stitle4 = new TLatex(0.0,0.0,"");
363  TLatex * stitle5 = new TLatex(0.0,0.0,"");
364  TLatex * stitle6 = new TLatex(0.0,0.0,"");
365 
366  title->SetNDC();
367  stitle1->SetNDC();
368  stitle2->SetNDC();
369  stitle3->SetNDC();
370  stitle4->SetNDC();
371  stitle5->SetNDC();
372  stitle6->SetNDC();
373 
374  title->SetTextSize(0.32);
375  stitle1->SetTextSize(0.35);
376  stitle2->SetTextSize(0.35);
377  stitle3->SetTextSize(0.35);
378  stitle4->SetTextSize(0.35);
379  stitle5->SetTextSize(0.35);
380  stitle6->SetTextSize(0.35);
381 
382  canvas->cd();
383  c1->Draw();
384  stitlepad1->Draw();
385  c2->Draw();
386  stitlepad2->Draw();
387  c3->Draw();
388  stitlepad3->Draw();
389  c4->Draw();
390  stitlepad4->Draw();
391  c5->Draw();
392  stitlepad5->Draw();
393  c6->Draw();
394  stitlepad6->Draw();
395  titlepad->Draw();
396 
397  //x plots
398  c1->Divide(4,1);
399  c1->cd(1);
400  hCartesianAveDiff[0]->Draw("colz");
401  c1->cd(2);
402  hCartesianAveDiff[1]->Draw("colz");
403  c1->cd(3);
404  hCartesianShiftDifference[0]->Draw();
405  //c1->cd(4)->Clear();
406  c1->cd(4);
407  //hCMmodelSliceRvTrue->Draw("colz");
408  hSamplePerBinRZ->Draw("colz");
409 
410  //y plots
411  c2->Divide(4,1);
412  c2->cd(1);
413  hCartesianAveDiff[2]->Draw("colz");
414  c2->cd(2);
415  hCartesianAveDiff[3]->Draw("colz");
416  c2->cd(3);
417  hCartesianShiftDifference[1]->Draw();
418  //c2->cd(4)->Clear();
419  c2->cd(4);
420  //hStripesPerBin->Draw("colz");
421  hSamplePerBinXY->Draw("colz");
422 
423  //r cart
424  c3->Divide(4,1);
425  c3->cd(1);
426  hCylindricalAveDiff[0]->Draw("colz");
427  c3->cd(2);
428  hCylindricalAveDiff[1]->Draw("colz");
429  c3->cd(3);
430  hCylindricalShiftDifference[0]->Draw();
431  c3->cd(4);
432  hRShiftTrue->Draw();
433 
434  //phi cart
435  c4->Divide(4,1);
436  c4->cd(1);
437  hCylindricalAveDiff[2]->Draw("colz");
438  c4->cd(2);
439  hCylindricalAveDiff[3]->Draw("colz");
440  c4->cd(3);
441  hCylindricalShiftDifference[1]->Draw();
442  c4->cd(4);
443  hPhiShiftTrue->Draw();
444 
445  //r to true comparison
446  c5->Divide(4,1);
447  c5->cd(1);
448  hCompareRTrue->Draw("colz");
449  c5->cd(2);
450  hRDiffvR->Draw("colz");
451  c5->cd(3);
452  hRDiffvZ->Draw("colz");
453  c5->cd(4);
454  hRDiffvPhi->Draw("colz");
455 
456  //phi to true comparison
457  c6->Divide(4,1);
458  c6->cd(1);
459  hComparePhiTrue->Draw("colz");
460  c6->cd(2);
461  hPhiDiffvR->Draw("colz");
462  c6->cd(3);
463  hPhiDiffvZ->Draw("colz");
464  c6->cd(4);
465  hPhiDiffvPhi->Draw("colz");
466 
467  titlepad->cd();
468  titlepad->Clear();
469  title->DrawLatex(0.01,0.4,Form("Event %d; %s", ifile, sourcefilename.Data()));
470  title->Draw();
471 
472  stitlepad1->cd();
473  stitlepad1->Clear();
474  stitle1->DrawLatex(0.45,0.2,"X Model");
475  stitle1->Draw();
476 
477  stitlepad2->cd();
478  stitlepad2->Clear();
479  stitle2->DrawLatex(0.45,0.2,"Y Model");
480  stitle2->Draw();
481 
482  stitlepad3->cd();
483  stitlepad3->Clear();
484  stitle3->DrawLatex(0.45,0.2,"R Model");
485  stitle3->Draw();
486 
487  stitlepad4->cd();
488  stitlepad4->Clear();
489  stitle4->DrawLatex(0.45,0.2,"Phi Model");
490  stitle4->Draw();
491 
492  stitlepad5->cd();
493  stitlepad5->Clear();
494  stitle5->DrawLatex(0.4,0.2,"Comparing R Model to True");
495  stitle5->Draw();
496 
497  stitlepad6->cd();
498  stitlepad6->Clear();
499  stitle6->DrawLatex(0.4,0.2,"Comparing Phi Model to True");
500  stitle6->Draw();
501 
502  if(ifile == 0){
503  //if(ifile == 1){
504  canvas->Print("CMDistortionAnalysisCart.pdf(","pdf");
505  } else if((ifile == 1) || (ifile == nEvents - 1)){
506  canvas->Print("CMDistortionAnalysisCart.pdf","pdf");
507  }
508  }
509 
510  TCanvas *summary = new TCanvas("summary","ShiftPlotsSummary",2000,3000);
511 
512  TPad *sumtitlepad = new TPad("sumtitlepad","",0.0,0.96,1.0,1.0);
513  TPad *sumplots = new TPad("sumplotspad","",0.0,0.0,1.0,0.96);
514 
515  TLatex *sumtitle = new TLatex(0.0,0.0,"");
516 
517  sumtitle->SetNDC();
518  sumtitle->SetTextSize(0.4);
519 
520  summary->cd();
521  sumplots->Draw();
522  sumtitlepad->Draw();
523 
524  sumplots->Divide(4,6);
525  sumplots->cd(1);
526  hDifferenceMeanR->Draw();
527  sumplots->cd(2);
528  hDifferenceStdDevR->Draw();
529  sumplots->cd(3);
530  hTrueMeanR->Draw();
531  sumplots->cd(4);
532  hTrueStdDevR->Draw();
533  sumplots->cd(5);
534  hDifferenceMeanPhi->Draw();
535  sumplots->cd(6);
536  hDifferenceStdDevPhi->Draw();
537  sumplots->cd(7);
538  hTrueMeanPhi->Draw();
539  sumplots->cd(8);
540  hTrueStdDevPhi->Draw();
541  sumplots->cd(9);
542  sumplots->cd(10)->Clear();
543  sumplots->cd(11)->Clear();
544  sumplots->cd(12)->Clear();
545  sumplots->cd(13)->Clear();
546  sumplots->cd(14)->Clear();
547  sumplots->cd(15)->Clear();
548  sumplots->cd(16)->Clear();
549  sumplots->cd(17)->Clear();
550  sumplots->cd(18)->Clear();
551  sumplots->cd(19)->Clear();
552  sumplots->cd(20)->Clear();
553  sumplots->cd(21)->Clear();
554  sumplots->cd(22)->Clear();
555  sumplots->cd(23)->Clear();
556  sumplots->cd(24)->Clear();
557 
558  sumtitlepad->cd();
559  sumtitlepad->Clear();
560  sumtitle->DrawLatex(0.4,0.4,"Summary of Events");
561  summary->Print("CMDistortionAnalysisCart.pdf)","pdf");
562 
563  return 0;
564 }