tpc/cm/CMDistortionAnalysisCart.C

0001 //step 3 in cart coords
0002 #include <iostream>
0003 #include <cmath>
0004 #include <vector>
0005 #include "TMath.h"
0006 #include "TVector3.h"
0007 #include "TTree.h"
0008
0009 using namespace std;
0010
0011 class Shifter {
0012 public:
0013 Shifter(TString sourcefilename);
0014   TFile *forward, *average;
0015   TH3F *hX, *hY, *hZ, *hR, *hPhi, *hXave, *hYave, *hZave, *hRave, *hPhiave;
0016 };
0017
0018 Shifter::Shifter(TString sourcefilename){
0019   //single event distortion file
0020   forward=TFile::Open(sourcefilename,"READ");
0021
0022   hX=(TH3F*)forward->Get("hIntDistortionPosX");
0023   hY=(TH3F*)forward->Get("hIntDistortionPosY");
0024   hZ=(TH3F*)forward->Get("hIntDistortionPosZ");
0025
0026   hR=(TH3F*)forward->Get("hIntDistortionPosR");
0027   hPhi=(TH3F*)forward->Get("hIntDistortionPosP");
0028
0029   //average distortion file
0030   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");
0031
0032   hXave=(TH3F*)average->Get("hIntDistortionPosX");
0033   hYave=(TH3F*)average->Get("hIntDistortionPosY");
0034   hZave=(TH3F*)average->Get("hIntDistortionPosZ");
0035
0036   hRave=(TH3F*)average->Get("hIntDistortionPosR");
0037   hPhiave=(TH3F*)average->Get("hIntDistortionPosP");
0038
0039   //subtract average from total distortions to study fluctuations
0040   hX->Add(hXave,-1);
0041   hY->Add(hYave,-1);
0042   hZ->Add(hZave,-1);
0043
0044   hR->Add(hRave,-1);
0045   hPhi->Add(hPhiave,-1);
0046 }
0047
0048 int CMDistortionAnalysisCart(int nMaxEvents = -1) {
0049   Shifter *shifter;
0050   int nbins = 35;
0051   double x, y, z;
0052   double low = -80.0;
0053   double high = 80.0;
0054   double deltaX, deltaY, deltaZ, deltaR, deltaPhi;
0055   int nEvents;
0056
0057   TCanvas *canvas=new TCanvas("canvas","CMDistortionAnalysisCart",2000,3000);
0058
0059   int nsumbins = 20;
0060   int minsum = -10;
0061   int maxsum = 10;
0062
0063   //set up summary plots
0064   TH1F *hDifferenceMeanR = new TH1F("hDifferenceMeanR", "Average Difference between R Model and True of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
0065     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);
0066
0067     TH1F *hTrueMeanR = new TH1F("hTrueMeanR", "Mean True R Distortion Model of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
0068     TH1F *hTrueStdDevR = new TH1F("hTrueStdDevR", "Std Dev of True R Distortion Model of All Events (R > 30); #Delta R (#mum)", nsumbins, minsum, maxsum);
0069
0070     TH1F *hDifferenceMeanPhi = new TH1F("hDifferenceMeanPhi", "Average Difference between Phi Model and True of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
0071     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);
0072
0073     TH1F *hTrueMeanPhi = new TH1F("hTrueMeanPhi", "Mean True Phi Distortion Model of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
0074     TH1F *hTrueStdDevPhi = new TH1F("hTrueStdDevPhi", "Std Dev of True Phi Distortion Model of All Events (R > 30); #Delta Phi (#mum)", nsumbins, minsum, maxsum);
0075
0076     const char * inputpattern="/sphenix/user/rcorliss/distortion_maps/2021.04/*h_Charge_*.root"; //updated
0077
0078   //find all files that match the input string (includes wildcards)
0079   TFileCollection *filelist=new TFileCollection();
0080   filelist->Add(inputpattern);
0081   TString sourcefilename;
0082   /*
0083    //how many events
0084   if (nMaxEvents<0){
0085     nEvents=filelist->GetNFiles();
0086   } else if(nMaxEvents<filelist->GetNFiles()){
0087     nEvents=nMaxEvents;
0088   } else {
0089     nEvents= filelist->GetNFiles();
0090     }
0091   */
0092   nEvents = 2;
0093
0094   for (int ifile=0;ifile < nEvents;ifile++){
0095     //for each file, find all histograms in that file.
0096     sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
0097
0098     //create shifter
0099     shifter = new Shifter(sourcefilename);
0100
0101     TFile *plots;
0102
0103     plots=TFile::Open(Form("CMModelsCart_Event%d.root",ifile),"READ");
0104
0105     TH3F *hCartCMModel[3];
0106     hCartCMModel[0]=(TH3F*)plots->Get("hCMModelX");
0107     hCartCMModel[1]=(TH3F*)plots->Get("hCMModelY");
0108     hCartCMModel[2]=(TH3F*)plots->Get("hCMModelZ");
0109
0110     TH3F *hCylCMModel[2];
0111     hCylCMModel[0]=(TH3F*)plots->Get("hCMModelRCart");
0112     hCylCMModel[1]=(TH3F*)plots->Get("hCMModelPhiCart");
0113
0114     //phi,r binning
0115     TH3F *hCartCMModelPhiR[3];
0116     hCartCMModelPhiR[0]=(TH3F*)plots->Get("hCMModelX_PhiR");
0117     hCartCMModelPhiR[1]=(TH3F*)plots->Get("hCMModelY_PhiR");
0118     hCartCMModelPhiR[2]=(TH3F*)plots->Get("hCMModelZ_PhiR");
0119
0120     TH3F *hCylCMModelPhiR[2];
0121     hCylCMModelPhiR[0]=(TH3F*)plots->Get("hCMModelR_PhiR");
0122     hCylCMModelPhiR[1]=(TH3F*)plots->Get("hCMModelPhi_PhiR");
0123
0124     //for forward only
0125
0126     //same range and bins for each coordinate, binned in cm
0127     //hardcoded numbers from average distortion file's hIntDistortionPosX
0128     int nphi = 82;
0129     int nr = 54;
0130     int nz = 82;
0131
0132     double minphi = -0.078539819;
0133     double minr = 18.884615;
0134     double minz = -1.3187500;
0135
0136     double maxphi = 6.3617253;
0137     double maxr = 79.115387;
0138     double maxz = 106.81875;
0139
0140     double rshiftcart, phishiftcart;
0141
0142     int ndiff = 300;
0143     int mindiff = -20;
0144     int maxdiff = 20;
0145
0146     TH1F *hCartesianShiftDifference[3];
0147     hCartesianShiftDifference[0] = new TH1F("hShiftDifferenceX", "Difference between CM Model X and True (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
0148     hCartesianShiftDifference[1] = new TH1F("hShiftDifferenceY", "Difference between CM Model Y and True (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
0149     hCartesianShiftDifference[2] = new TH1F("hShiftDifferenceZ", "Difference between CM Model Z and True (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
0150
0151     TH1F *hCylindricalShiftDifference[2];
0152     hCylindricalShiftDifference[0] = new TH1F("hShiftDifferenceRCart", "Difference between CM Model R from Cartesian and True (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
0153     hCylindricalShiftDifference[1] = new TH1F("hShiftDifferencePhiCart", "Difference between CM Model Phi from Cartesian and True (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
0154
0155     TH1F *hRShiftTrue = new TH1F("hRShiftTrue", "True R Distortion Model (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
0156     TH1F *hPhiShiftTrue = new TH1F("hPhiShiftTrue", "True Phi Distortion Model (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
0157
0158     TH2F *hCartesianDiff[6];
0159     hCartesianDiff[0] = new TH2F("hDiffXYX", "Difference in XY for CM Model X; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0160     hCartesianDiff[1] = new TH2F("hDiffRZX", "Difference in RZ for CM Model X; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0161     hCartesianDiff[2] = new TH2F("hDiffXYY", "Difference in XY for CM Model Y; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0162     hCartesianDiff[3] = new TH2F("hDiffRZY", "Difference in RZ for CM Model Y; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0163     hCartesianDiff[4] = new TH2F("hDiffXYZ", "Difference in XY for CM Model Z; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0164     hCartesianDiff[5] = new TH2F("hDiffRZZ", "Difference in RZ for CM Model Z; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0165
0166     TH2F *hCylindricalDiff[4];
0167     hCylindricalDiff[0] = new TH2F("hDiffXYRCart", "Difference in XY for CM Model R from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0168     hCylindricalDiff[1] = new TH2F("hDiffRZRCart", "Difference in RZ for CM Model R from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
0169     hCylindricalDiff[2] = new TH2F("hDiffXYPhiCart", "Difference in XY for CM Model Phi from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0170     hCylindricalDiff[3] = new TH2F("hDiffRZPhiCart", "Difference in RZ for CM Model Phi from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
0171
0172     TH2F *hCartesianAveDiff[6];
0173     hCartesianAveDiff[0] = new TH2F("hAveDiffXYX", "X Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0174     hCartesianAveDiff[1] = new TH2F("hAveDiffRZX", "X Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0175     hCartesianAveDiff[2] = new TH2F("hAveDiffXYY", "Y Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0176     hCartesianAveDiff[3] = new TH2F("hAveDiffRZY", "Y Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0177     hCartesianAveDiff[4] = new TH2F("hAveDiffXYZ", "Z Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0178     hCartesianAveDiff[5] = new TH2F("hAveDiffRZZ", "Z Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0179
0180     TH2F *hCylindricalAveDiff[4];
0181     hCylindricalAveDiff[0] = new TH2F("hAveDiffXYRCart", "R Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0182     hCylindricalAveDiff[1] = new TH2F("hAveDiffRZRCart", "R Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0183     hCylindricalAveDiff[2] = new TH2F("hAveDiffXYPhiCart", "Phi Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0184     hCylindricalAveDiff[3] = new TH2F("hAveDiffRZPhiCart", "Phi Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0185
0186     TH2F *hSamplePerBinXY = new TH2F("hSamplePerBinXY", "Filling each xy bin; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0187     TH2F *hSamplePerBinRZ = new TH2F("hSamplePerBinRZ", "Filling each rz bin; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0188
0189     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);
0190     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);
0191
0192     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);
0193     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);
0194     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);
0195
0196     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);
0197     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);
0198     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);
0199
0200     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);
0201     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);
0202     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);
0203
0204     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);
0205     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);
0206     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);
0207
0208     for(int i = 1; i < nphi - 1; i++){
0209       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0210       for(int j = 1; j < nr - 1; j++){
0211     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0212     for(int k = 1; k < nz - 1; k++){
0213       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0214
0215       double shiftrecoCart[3];
0216       double shifttrueCart[3];
0217       double differenceCart[3];
0218
0219       double shiftrecoCyl[2];
0220       double shifttrueCyl[2];
0221       double differenceCyl[2];
0222
0223       double differenceR, differencePhi;
0224
0225       int bin = hCartCMModel[0]->FindBin(phi,r,z); //same for all
0226
0227       if((r > 30.0) && (r < 76.0)){
0228         //x y and z
0229         shifttrueCart[0] = (shifter->hX->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0230         shifttrueCart[1] = (shifter->hY->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0231         shifttrueCart[2] = (shifter->hZ->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0232         for(int l = 0; l < 3; l ++){
0233           shiftrecoCart[l] =  (hCartCMModel[l]->GetBinContent(bin))*(1e4);
0234
0235           differenceCart[l] = shiftrecoCart[l] - shifttrueCart[l];
0236
0237           hCartesianShiftDifference[l]->Fill(differenceCart[l]);
0238         }
0239
0240         //r from cart
0241         shiftrecoCyl[0] =  (hCylCMModel[0]->GetBinContent(bin))*(1e4);
0242         shifttrueCyl[0] = (shifter->hR->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0243         differenceCyl[0] = shiftrecoCyl[0] - shifttrueCyl[0];
0244         hCylindricalShiftDifference[0]->Fill(differenceCyl[0]);
0245
0246         //phi from cart
0247         shiftrecoCyl[1] = r*(1e4)*(hCylCMModel[1]->GetBinContent(bin));
0248         shifttrueCyl[1] = (shifter->hPhi->Interpolate(phi,r,z))*(1e4);
0249         differenceCyl[1] = (shiftrecoCyl[1] - shifttrueCyl[1]);
0250         hCylindricalShiftDifference[1]->Fill(differenceCyl[1]);
0251
0252         hRShiftTrue->Fill(shifttrueCyl[0]);
0253         hPhiShiftTrue->Fill(shifttrueCyl[1]);
0254
0255         double x = r*cos(phi);
0256         double y = r*sin(phi);
0257
0258         //x
0259         hCartesianDiff[0]->Fill(x,y, differenceCart[0]);
0260         hCartesianDiff[1]->Fill(z,r, differenceCart[0]);
0261         //y
0262         hCartesianDiff[2]->Fill(x,y, differenceCart[1]);
0263         hCartesianDiff[3]->Fill(z,r, differenceCart[1]);
0264         //z
0265         hCartesianDiff[4]->Fill(x,y, differenceCart[2]);
0266         hCartesianDiff[5]->Fill(z,r, differenceCart[2]);
0267
0268         //r cart
0269         hCylindricalDiff[0]->Fill(x,y, differenceCyl[0]);
0270         hCylindricalDiff[1]->Fill(z,r, differenceCyl[0]);
0271         //phi cart
0272         hCylindricalDiff[2]->Fill(x,y, differenceCyl[1]);
0273         hCylindricalDiff[3]->Fill(z,r, differenceCyl[1]);
0274
0275         hCompareRTrue->Fill(shiftrecoCyl[0],shifttrueCyl[0]);
0276         hComparePhiTrue->Fill(shiftrecoCyl[1],shifttrueCyl[1]);
0277
0278         hRDiffvR->Fill(r,differenceCyl[0],1);
0279         hRDiffvPhi->Fill(phi,differenceCyl[0],1);
0280         hRDiffvZ->Fill(z,differenceCyl[0],1);
0281
0282         hPhiDiffvR->Fill(r,differenceCyl[1],1);
0283         hPhiDiffvPhi->Fill(phi,differenceCyl[1],1);
0284         hPhiDiffvZ->Fill(z,differenceCyl[1],1);
0285
0286         hSamplePerBinXY->Fill(x,y,1);
0287
0288         hSamplePerBinRZ->Fill(z,r,1);
0289       }
0290     }
0291       }
0292     }
0293
0294     //average over z
0295     for (int m = 0; m < 6; m = m+2){
0296       hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinXY);
0297     }
0298     for (int m = 0; m < 4; m = m+2){
0299       hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinXY);
0300     }
0301
0302     //average over phi
0303     for (int m = 1; m < 6; m = m+2){
0304       hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinRZ);
0305     }
0306     for (int m = 1; m < 4; m = m+2){
0307       hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinRZ);
0308     }
0309
0310     //summary plots
0311     hDifferenceMeanR->Fill(hCylindricalShiftDifference[0]->GetMean(1));
0312     hDifferenceStdDevR->Fill(hCylindricalShiftDifference[0]->GetStdDev(1));
0313
0314     hTrueMeanR->Fill(hRShiftTrue->GetMean(1));
0315     hTrueStdDevR->Fill(hRShiftTrue->GetStdDev(1));
0316
0317     hDifferenceMeanPhi->Fill(hCylindricalShiftDifference[1]->GetMean(1));
0318     hDifferenceStdDevPhi->Fill(hCylindricalShiftDifference[1]->GetStdDev(1));
0319
0320     hTrueMeanPhi->Fill(hPhiShiftTrue->GetMean(1));
0321     hTrueStdDevPhi->Fill(hPhiShiftTrue->GetStdDev(1));
0322
0323     for (int m = 0; m < 6; m++){
0324       hCartesianAveDiff[m]->SetStats(0);
0325     }
0326     for (int m = 0; m < 4; m++){
0327       hCylindricalAveDiff[m]->SetStats(0);
0328     }
0329
0330     hCompareRTrue->SetStats(0);
0331     hComparePhiTrue->SetStats(0);
0332
0333     hRDiffvR->SetStats(0);
0334     hRDiffvZ->SetStats(0);
0335     hRDiffvPhi->SetStats(0);
0336
0337     hPhiDiffvR->SetStats(0);
0338     hPhiDiffvZ->SetStats(0);
0339     hPhiDiffvPhi->SetStats(0);
0340
0341     TPad *c1=new TPad("c1","",0.0,0.8,1.0,0.93); //can i do an array of pads?
0342     TPad *c2=new TPad("c2","",0.0,0.64,1.0,0.77);
0343     TPad *c3=new TPad("c3","",0.0,0.48,1.0,0.61);
0344     TPad *c4=new TPad("c4","",0.0,0.32,1.0,0.45);
0345     TPad *c5=new TPad("c5","",0.0,0.16,1.0,0.29);
0346     TPad *c6=new TPad("c6","",0.0,0.0,1.0,0.13);
0347
0348     TPad *titlepad=new TPad("titlepad","",0.0,0.96,1.0,1.0);
0349
0350     TPad *stitlepad1=new TPad("stitlepad1","",0.0,0.93,1.0,0.96);
0351     TPad *stitlepad2=new TPad("stitlepad2","",0.0,0.77,1.0,0.8);
0352     TPad *stitlepad3=new TPad("stitlepad3","",0.0,0.61,1.0,0.64);
0353     TPad *stitlepad4=new TPad("stitlepad4","",0.0,0.45,1.0,0.48);
0354     TPad *stitlepad5=new TPad("stitlepad5","",0.0,0.29,1.0,0.32);
0355     TPad *stitlepad6=new TPad("stitlepad6","",0.0,0.13,1.0,0.16);
0356
0357     TLatex * title = new TLatex(0.0,0.0,"");
0358
0359     TLatex * stitle1 = new TLatex(0.0,0.0,""); //array?
0360     TLatex * stitle2 = new TLatex(0.0,0.0,"");
0361     TLatex * stitle3 = new TLatex(0.0,0.0,"");
0362     TLatex * stitle4 = new TLatex(0.0,0.0,"");
0363     TLatex * stitle5 = new TLatex(0.0,0.0,"");
0364     TLatex * stitle6 = new TLatex(0.0,0.0,"");
0365
0366     title->SetNDC();
0367     stitle1->SetNDC();
0368     stitle2->SetNDC();
0369     stitle3->SetNDC();
0370     stitle4->SetNDC();
0371     stitle5->SetNDC();
0372     stitle6->SetNDC();
0373
0374     title->SetTextSize(0.32);
0375     stitle1->SetTextSize(0.35);
0376     stitle2->SetTextSize(0.35);
0377     stitle3->SetTextSize(0.35);
0378     stitle4->SetTextSize(0.35);
0379     stitle5->SetTextSize(0.35);
0380     stitle6->SetTextSize(0.35);
0381
0382     canvas->cd();
0383     c1->Draw();
0384     stitlepad1->Draw();
0385     c2->Draw();
0386     stitlepad2->Draw();
0387     c3->Draw();
0388     stitlepad3->Draw();
0389     c4->Draw();
0390     stitlepad4->Draw();
0391     c5->Draw();
0392     stitlepad5->Draw();
0393     c6->Draw();
0394     stitlepad6->Draw();
0395     titlepad->Draw();
0396
0397     //x plots
0398     c1->Divide(4,1);
0399     c1->cd(1);
0400     hCartesianAveDiff[0]->Draw("colz");
0401     c1->cd(2);
0402     hCartesianAveDiff[1]->Draw("colz");
0403     c1->cd(3);
0404     hCartesianShiftDifference[0]->Draw();
0405     //c1->cd(4)->Clear();
0406     c1->cd(4);
0407     //hCMmodelSliceRvTrue->Draw("colz");
0408     hSamplePerBinRZ->Draw("colz");
0409
0410     //y plots
0411     c2->Divide(4,1);
0412     c2->cd(1);
0413     hCartesianAveDiff[2]->Draw("colz");
0414     c2->cd(2);
0415     hCartesianAveDiff[3]->Draw("colz");
0416     c2->cd(3);
0417     hCartesianShiftDifference[1]->Draw();
0418     //c2->cd(4)->Clear();
0419     c2->cd(4);
0420     //hStripesPerBin->Draw("colz");
0421     hSamplePerBinXY->Draw("colz");
0422
0423     //r cart
0424     c3->Divide(4,1);
0425     c3->cd(1);
0426     hCylindricalAveDiff[0]->Draw("colz");
0427     c3->cd(2);
0428     hCylindricalAveDiff[1]->Draw("colz");
0429     c3->cd(3);
0430     hCylindricalShiftDifference[0]->Draw();
0431     c3->cd(4);
0432     hRShiftTrue->Draw();
0433
0434     //phi cart
0435     c4->Divide(4,1);
0436     c4->cd(1);
0437     hCylindricalAveDiff[2]->Draw("colz");
0438     c4->cd(2);
0439     hCylindricalAveDiff[3]->Draw("colz");
0440     c4->cd(3);
0441     hCylindricalShiftDifference[1]->Draw();
0442     c4->cd(4);
0443     hPhiShiftTrue->Draw();
0444
0445     //r to true comparison
0446     c5->Divide(4,1);
0447     c5->cd(1);
0448     hCompareRTrue->Draw("colz");
0449     c5->cd(2);
0450     hRDiffvR->Draw("colz");
0451     c5->cd(3);
0452     hRDiffvZ->Draw("colz");
0453     c5->cd(4);
0454     hRDiffvPhi->Draw("colz");
0455
0456     //phi to true comparison
0457     c6->Divide(4,1);
0458     c6->cd(1);
0459     hComparePhiTrue->Draw("colz");
0460     c6->cd(2);
0461     hPhiDiffvR->Draw("colz");
0462     c6->cd(3);
0463     hPhiDiffvZ->Draw("colz");
0464     c6->cd(4);
0465     hPhiDiffvPhi->Draw("colz");
0466
0467     titlepad->cd();
0468     titlepad->Clear();
0469     title->DrawLatex(0.01,0.4,Form("Event %d; %s", ifile, sourcefilename.Data()));
0470     title->Draw();
0471
0472     stitlepad1->cd();
0473     stitlepad1->Clear();
0474     stitle1->DrawLatex(0.45,0.2,"X Model");
0475     stitle1->Draw();
0476
0477     stitlepad2->cd();
0478     stitlepad2->Clear();
0479     stitle2->DrawLatex(0.45,0.2,"Y Model");
0480     stitle2->Draw();
0481
0482     stitlepad3->cd();
0483     stitlepad3->Clear();
0484     stitle3->DrawLatex(0.45,0.2,"R Model");
0485     stitle3->Draw();
0486
0487     stitlepad4->cd();
0488     stitlepad4->Clear();
0489     stitle4->DrawLatex(0.45,0.2,"Phi Model");
0490     stitle4->Draw();
0491
0492     stitlepad5->cd();
0493     stitlepad5->Clear();
0494     stitle5->DrawLatex(0.4,0.2,"Comparing R Model to True");
0495     stitle5->Draw();
0496
0497     stitlepad6->cd();
0498     stitlepad6->Clear();
0499     stitle6->DrawLatex(0.4,0.2,"Comparing Phi Model to True");
0500     stitle6->Draw();
0501
0502     if(ifile == 0){
0503       //if(ifile == 1){
0504       canvas->Print("CMDistortionAnalysisCart.pdf(","pdf");
0505     } else if((ifile == 1) || (ifile == nEvents - 1)){
0506       canvas->Print("CMDistortionAnalysisCart.pdf","pdf");
0507     }
0508   }
0509
0510   TCanvas *summary = new TCanvas("summary","ShiftPlotsSummary",2000,3000);
0511
0512   TPad *sumtitlepad = new TPad("sumtitlepad","",0.0,0.96,1.0,1.0);
0513   TPad *sumplots = new TPad("sumplotspad","",0.0,0.0,1.0,0.96);
0514
0515   TLatex *sumtitle = new TLatex(0.0,0.0,"");
0516
0517   sumtitle->SetNDC();
0518   sumtitle->SetTextSize(0.4);
0519
0520   summary->cd();
0521   sumplots->Draw();
0522   sumtitlepad->Draw();
0523
0524   sumplots->Divide(4,6);
0525   sumplots->cd(1);
0526   hDifferenceMeanR->Draw();
0527   sumplots->cd(2);
0528   hDifferenceStdDevR->Draw();
0529   sumplots->cd(3);
0530   hTrueMeanR->Draw();
0531   sumplots->cd(4);
0532   hTrueStdDevR->Draw();
0533   sumplots->cd(5);
0534   hDifferenceMeanPhi->Draw();
0535   sumplots->cd(6);
0536   hDifferenceStdDevPhi->Draw();
0537   sumplots->cd(7);
0538   hTrueMeanPhi->Draw();
0539   sumplots->cd(8);
0540   hTrueStdDevPhi->Draw();
0541   sumplots->cd(9);
0542   sumplots->cd(10)->Clear();
0543   sumplots->cd(11)->Clear();
0544   sumplots->cd(12)->Clear();
0545   sumplots->cd(13)->Clear();
0546   sumplots->cd(14)->Clear();
0547   sumplots->cd(15)->Clear();
0548   sumplots->cd(16)->Clear();
0549   sumplots->cd(17)->Clear();
0550   sumplots->cd(18)->Clear();
0551   sumplots->cd(19)->Clear();
0552   sumplots->cd(20)->Clear();
0553   sumplots->cd(21)->Clear();
0554   sumplots->cd(22)->Clear();
0555   sumplots->cd(23)->Clear();
0556   sumplots->cd(24)->Clear();
0557
0558   sumtitlepad->cd();
0559   sumtitlepad->Clear();
0560   sumtitle->DrawLatex(0.4,0.4,"Summary of Events");
0561   summary->Print("CMDistortionAnalysisCart.pdf)","pdf");
0562
0563   return 0;
0564 }