tpc/cm/CMDistortionAnalysis.C

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