tpc/cm/CMDistortionAnalysisCart.C

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