tpc/cm/CMDistortionAnalysis.C

0001 //step 3
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 CMDistortionAnalysis(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","CMDistortionAnalysisPhiR",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("CMModels_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
0147     TH1F *hCartesianShiftDifference[3];
0148     hCartesianShiftDifference[0] = new TH1F("hShiftDifferenceX", "Difference between CM Model X and True (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
0149     hCartesianShiftDifference[1] = new TH1F("hShiftDifferenceY", "Difference between CM Model Y and True (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
0150     hCartesianShiftDifference[2] = new TH1F("hShiftDifferenceZ", "Difference between CM Model Z and True (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
0151
0152     //phi,r binning
0153      TH1F *hCartesianShiftDifferencePhiR[3];
0154     hCartesianShiftDifferencePhiR[0] = new TH1F("hShiftDifferenceX_PhiR", "Difference between CM Model X and True, Phi,R binning (R > 30); #Delta X (#mum)", ndiff, mindiff, maxdiff);
0155     hCartesianShiftDifferencePhiR[1] = new TH1F("hShiftDifferenceY_PhiR", "Difference between CM Model Y and True, Phi,R binning (R > 30); #Delta Y (#mum)", ndiff, mindiff, maxdiff);
0156     hCartesianShiftDifferencePhiR[2] = new TH1F("hShiftDifferenceZ_PhiR", "Difference between CM Model Z and True, Phi,R binning (R > 30); #Delta Z (#mum)", ndiff, mindiff, maxdiff);
0157
0158     TH1F *hCylindricalShiftDifference[2];
0159     hCylindricalShiftDifference[0] = new TH1F("hShiftDifferenceRCart", "Difference between CM Model R from Cartesian and True (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
0160     hCylindricalShiftDifference[1] = new TH1F("hShiftDifferencePhiCart", "Difference between CM Model Phi from Cartesian and True (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
0161
0162     //phi,r binning
0163     TH1F *hCylindricalShiftDifferencePhiR[2];
0164     hCylindricalShiftDifferencePhiR[0] = new TH1F("hShiftDifferenceR_PhiR", "Difference between CM Model R and True, Phi,R binning (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
0165     hCylindricalShiftDifferencePhiR[1] = new TH1F("hShiftDifferencePhi_PhiR", "Difference between CM Model Phi and True, Phi,R binning (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
0166
0167     TH1F *hRShiftTrue = new TH1F("hRShiftTrue", "True R Distortion Model (R > 30); #Delta R (#mum)", ndiff, mindiff, maxdiff);
0168     TH1F *hPhiShiftTrue = new TH1F("hPhiShiftTrue", "True Phi Distortion Model (R > 30); #Delta Phi (#mum)", ndiff, mindiff, maxdiff);
0169
0170     TH2F *hCartesianDiff[6];
0171     hCartesianDiff[0] = new TH2F("hDiffXYX", "Difference in XY for CM Model X; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0172     hCartesianDiff[1] = new TH2F("hDiffRZX", "Difference in RZ for CM Model X; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0173     hCartesianDiff[2] = new TH2F("hDiffXYY", "Difference in XY for CM Model Y; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0174     hCartesianDiff[3] = new TH2F("hDiffRZY", "Difference in RZ for CM Model Y; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0175     hCartesianDiff[4] = new TH2F("hDiffXYZ", "Difference in XY for CM Model Z; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0176     hCartesianDiff[5] = new TH2F("hDiffRZZ", "Difference in RZ for CM Model Z; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0177
0178     //phi,r binning
0179      TH2F *hCartesianDiffPhiR[6];
0180     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);
0181     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);
0182     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);
0183     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);
0184     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);
0185     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);
0186
0187     TH2F *hCylindricalDiff[4];
0188     hCylindricalDiff[0] = new TH2F("hDiffXYRCart", "Difference in XY for CM Model R from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0189     hCylindricalDiff[1] = new TH2F("hDiffRZRCart", "Difference in RZ for CM Model R from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
0190     hCylindricalDiff[2] = new TH2F("hDiffXYPhiCart", "Difference in XY for CM Model Phi from Cartesian; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0191     hCylindricalDiff[3] = new TH2F("hDiffRZPhiCart", "Difference in RZ for CM Model Phi from Cartesian; z (cm); r (cm)",nz,minz,maxz,nr,minr,maxr);
0192
0193     //phi,r binning
0194     TH2F *hCylindricalDiffPhiR[4];
0195     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);
0196     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);
0197     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);
0198     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);
0199
0200     TH2F *hCartesianAveDiff[6];
0201     hCartesianAveDiff[0] = new TH2F("hAveDiffXYX", "X Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0202     hCartesianAveDiff[1] = new TH2F("hAveDiffRZX", "X Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0203     hCartesianAveDiff[2] = new TH2F("hAveDiffXYY", "Y Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0204     hCartesianAveDiff[3] = new TH2F("hAveDiffRZY", "Y Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0205     hCartesianAveDiff[4] = new TH2F("hAveDiffXYZ", "Z Model - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0206     hCartesianAveDiff[5] = new TH2F("hAveDiffRZZ", "Z Model - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0207
0208     //phi,r binning
0209     TH2F *hCartesianAveDiffPhiR[6];
0210     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);
0211     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);
0212     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);
0213     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);
0214     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);
0215     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);
0216
0217     TH2F *hCylindricalAveDiff[4];
0218     hCylindricalAveDiff[0] = new TH2F("hAveDiffXYRCart", "R Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0219     hCylindricalAveDiff[1] = new TH2F("hAveDiffRZRCart", "R Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0220     hCylindricalAveDiff[2] = new TH2F("hAveDiffXYPhiCart", "Phi Model from Cartesian - Truth Averaged Over z (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
0221     hCylindricalAveDiff[3] = new TH2F("hAveDiffRZPhiCart", "Phi Model from Cartesian - Truth Averaged Over phi (#mum); z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0222
0223     //phi,r binning
0224      TH2F *hCylindricalAveDiffPhiR[4];
0225     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);
0226     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);
0227     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);
0228     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);
0229
0230     TH2F *hSamplePerBinXY = new TH2F("hSamplePerBinXY", "Filling each xy bin; x (cm); y (cm)",nbins,low,high,nbins,low,high);
0231     TH2F *hSamplePerBinRZ = new TH2F("hSamplePerBinRZ", "Filling each rz bin; z (cm); r (cm)", nz,minz,maxz,nr,minr,maxr);
0232
0233     //phi,r binning
0234     TH2F *hSamplePerBinPhiR = new TH2F("hSamplePerBinPhiR", "Filling each PhiR bin; phi (rad); r (cm)",nphi,minphi,maxphi,nr,minr,maxr);
0235
0236     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);
0237     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);
0238
0239     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);
0240     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);
0241     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);
0242
0243     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);
0244     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);
0245     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);
0246
0247     //phi,r binning
0248     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);
0249     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);
0250
0251     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);
0252     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);
0253     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);
0254
0255     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);
0256     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);
0257     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);
0258
0259     for(int i = 1; i < nphi - 1; i++){
0260       double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
0261       for(int j = 1; j < nr - 1; j++){
0262     double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
0263     for(int k = 1; k < nz - 1; k++){
0264       double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
0265
0266       double shiftrecoCart[3];
0267       double shifttrueCart[3];
0268       double differenceCart[3];
0269
0270       double shiftrecoCyl[2];
0271       double shifttrueCyl[2];
0272       double differenceCyl[2];
0273
0274       double differenceR, differencePhi;
0275
0276       int bin = hCartCMModel[0]->FindBin(phi,r,z); //same for all
0277
0278       //phi,r binning
0279       double shiftrecoCartPhiR[3];
0280       double differenceCartPhiR[3];
0281       double shiftrecoCylPhiR[2];
0282       double differenceCylPhiR[2];
0283       double differenceR_PhiR, differencePhi_PhiR;
0284       int binPhiR = hCartCMModelPhiR[0]->FindBin(phi,r,z);
0285
0286       if((r > 30.0) && (r < 76.0)){
0287         //x y and z
0288         shifttrueCart[0] = (shifter->hX->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0289         shifttrueCart[1] = (shifter->hY->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0290         shifttrueCart[2] = (shifter->hZ->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0291         for(int l = 0; l < 3; l ++){
0292           shiftrecoCart[l] =  (hCartCMModel[l]->GetBinContent(bin))*(1e4);
0293
0294           differenceCart[l] = shiftrecoCart[l] - shifttrueCart[l];
0295
0296           hCartesianShiftDifference[l]->Fill(differenceCart[l]);
0297         }
0298
0299         //r from cart
0300         shiftrecoCyl[0] =  (hCylCMModel[0]->GetBinContent(bin))*(1e4);
0301         shifttrueCyl[0] = (shifter->hR->Interpolate(phi,r,z))*(1e4); //convert from cm to micron
0302         differenceCyl[0] = shiftrecoCyl[0] - shifttrueCyl[0];
0303         hCylindricalShiftDifference[0]->Fill(differenceCyl[0]);
0304
0305         //phi from cart
0306         shiftrecoCyl[1] = r*(1e4)*(hCylCMModel[1]->GetBinContent(bin));
0307         shifttrueCyl[1] = (shifter->hPhi->Interpolate(phi,r,z))*(1e4);
0308         differenceCyl[1] = (shiftrecoCyl[1] - shifttrueCyl[1]);
0309         hCylindricalShiftDifference[1]->Fill(differenceCyl[1]);
0310
0311         hRShiftTrue->Fill(shifttrueCyl[0]);
0312         hPhiShiftTrue->Fill(shifttrueCyl[1]);
0313
0314         double x = r*cos(phi);
0315         double y = r*sin(phi);
0316
0317         //x
0318         hCartesianDiff[0]->Fill(x,y, differenceCart[0]);
0319         hCartesianDiff[1]->Fill(z,r, differenceCart[0]);
0320         //y
0321         hCartesianDiff[2]->Fill(x,y, differenceCart[1]);
0322         hCartesianDiff[3]->Fill(z,r, differenceCart[1]);
0323         //z
0324         hCartesianDiff[4]->Fill(x,y, differenceCart[2]);
0325         hCartesianDiff[5]->Fill(z,r, differenceCart[2]);
0326
0327         //r cart
0328         hCylindricalDiff[0]->Fill(x,y, differenceCyl[0]);
0329         hCylindricalDiff[1]->Fill(z,r, differenceCyl[0]);
0330         //phi cart
0331         hCylindricalDiff[2]->Fill(x,y, differenceCyl[1]);
0332         hCylindricalDiff[3]->Fill(z,r, differenceCyl[1]);
0333
0334         hCompareRTrue->Fill(shiftrecoCyl[0],shifttrueCyl[0]);
0335         hComparePhiTrue->Fill(shiftrecoCyl[1],shifttrueCyl[1]);
0336
0337         hRDiffvR->Fill(r,differenceCyl[0],1);
0338         hRDiffvPhi->Fill(phi,differenceCyl[0],1);
0339         hRDiffvZ->Fill(z,differenceCyl[0],1);
0340
0341         hPhiDiffvR->Fill(r,differenceCyl[1],1);
0342         hPhiDiffvPhi->Fill(phi,differenceCyl[1],1);
0343         hPhiDiffvZ->Fill(z,differenceCyl[1],1);
0344
0345         hSamplePerBinXY->Fill(x,y,1);
0346
0347         hSamplePerBinRZ->Fill(z,r,1);
0348
0349         //phi,r binning
0350         //x y and z
0351         for(int l = 0; l < 3; l ++){
0352           shiftrecoCartPhiR[l] =  (hCartCMModelPhiR[l]->GetBinContent(binPhiR))*(1e4);
0353
0354           differenceCartPhiR[l] = shiftrecoCartPhiR[l] - shifttrueCart[l];
0355
0356           hCartesianShiftDifferencePhiR[l]->Fill(differenceCartPhiR[l]);
0357         }
0358
0359         //r
0360         shiftrecoCylPhiR[0] =  (hCylCMModelPhiR[0]->GetBinContent(binPhiR))*(1e4);
0361         differenceCylPhiR[0] = shiftrecoCylPhiR[0] - shifttrueCyl[0];
0362         hCylindricalShiftDifferencePhiR[0]->Fill(differenceCylPhiR[0]);
0363
0364         //phi
0365         shiftrecoCylPhiR[1] = r*(1e4)*(hCylCMModelPhiR[1]->GetBinContent(binPhiR));
0366         differenceCylPhiR[1] = (shiftrecoCylPhiR[1] - shifttrueCyl[1]);
0367         hCylindricalShiftDifferencePhiR[1]->Fill(differenceCylPhiR[1]);
0368
0369         //x
0370         hCartesianDiffPhiR[0]->Fill(phi,r, differenceCartPhiR[0]);
0371         hCartesianDiffPhiR[1]->Fill(z,r, differenceCartPhiR[0]);
0372         //y
0373         hCartesianDiffPhiR[2]->Fill(phi,r, differenceCartPhiR[1]);
0374         hCartesianDiffPhiR[3]->Fill(z,r, differenceCartPhiR[1]);
0375         //z
0376         hCartesianDiffPhiR[4]->Fill(phi,r, differenceCartPhiR[2]);
0377         hCartesianDiffPhiR[5]->Fill(z,r, differenceCartPhiR[2]);
0378
0379         //r
0380         hCylindricalDiffPhiR[0]->Fill(phi,r, differenceCylPhiR[0]);
0381         hCylindricalDiffPhiR[1]->Fill(z,r, differenceCylPhiR[0]);
0382         //phi
0383         hCylindricalDiffPhiR[2]->Fill(phi,r, differenceCylPhiR[1]);
0384         hCylindricalDiffPhiR[3]->Fill(z,r, differenceCylPhiR[1]);
0385
0386         hCompareRTrue_PhiR->Fill(shiftrecoCylPhiR[0],shifttrueCyl[0]);
0387         hComparePhiTrue_PhiR->Fill(shiftrecoCylPhiR[1],shifttrueCyl[1]);
0388
0389         hRDiffvR_PhiR->Fill(r,differenceCylPhiR[0],1);
0390         hRDiffvPhi_PhiR->Fill(phi,differenceCylPhiR[0],1);
0391         hRDiffvZ_PhiR->Fill(z,differenceCylPhiR[0],1);
0392
0393         hPhiDiffvR_PhiR->Fill(r,differenceCylPhiR[1],1);
0394         hPhiDiffvPhi_PhiR->Fill(phi,differenceCylPhiR[1],1);
0395         hPhiDiffvZ_PhiR->Fill(z,differenceCylPhiR[1],1);
0396
0397         hSamplePerBinPhiR->Fill(phi,r,1);
0398       }
0399     }
0400       }
0401     }
0402
0403     //average over z
0404     for (int m = 0; m < 6; m = m+2){
0405       hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinXY);
0406     }
0407     for (int m = 0; m < 4; m = m+2){
0408       hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinXY);
0409     }
0410
0411     //average over phi
0412     for (int m = 1; m < 6; m = m+2){
0413       hCartesianAveDiff[m]->Divide(hCartesianDiff[m],hSamplePerBinRZ);
0414     }
0415     for (int m = 1; m < 4; m = m+2){
0416       hCylindricalAveDiff[m]->Divide(hCylindricalDiff[m],hSamplePerBinRZ);
0417     }
0418
0419     //phi,r binning
0420     //average over z
0421     for (int m = 0; m < 6; m = m+2){
0422       hCartesianAveDiffPhiR[m]->Divide(hCartesianDiffPhiR[m],hSamplePerBinPhiR);
0423     }
0424     for (int m = 0; m < 4; m = m+2){
0425       hCylindricalAveDiffPhiR[m]->Divide(hCylindricalDiffPhiR[m],hSamplePerBinPhiR);
0426     }
0427
0428     //average over phi
0429     for (int m = 1; m < 6; m = m+2){
0430       hCartesianAveDiffPhiR[m]->Divide(hCartesianDiffPhiR[m],hSamplePerBinRZ);
0431     }
0432     for (int m = 1; m < 4; m = m+2){
0433       hCylindricalAveDiffPhiR[m]->Divide(hCylindricalDiffPhiR[m],hSamplePerBinRZ);
0434     }
0435
0436     //summary plots
0437     hDifferenceMeanR->Fill(hCylindricalShiftDifferencePhiR[0]->GetMean(1));
0438     hDifferenceStdDevR->Fill(hCylindricalShiftDifferencePhiR[0]->GetStdDev(1));
0439
0440     hTrueMeanR->Fill(hRShiftTrue->GetMean(1));
0441     hTrueStdDevR->Fill(hRShiftTrue->GetStdDev(1));
0442
0443     hDifferenceMeanPhi->Fill(hCylindricalShiftDifferencePhiR[1]->GetMean(1));
0444     hDifferenceStdDevPhi->Fill(hCylindricalShiftDifferencePhiR[1]->GetStdDev(1));
0445
0446     hTrueMeanPhi->Fill(hPhiShiftTrue->GetMean(1));
0447     hTrueStdDevPhi->Fill(hPhiShiftTrue->GetStdDev(1));
0448
0449     for (int m = 0; m < 6; m++){
0450       hCartesianAveDiff[m]->SetStats(0);
0451       hCartesianAveDiffPhiR[m]->SetStats(0);
0452     }
0453     for (int m = 0; m < 4; m++){
0454       hCylindricalAveDiff[m]->SetStats(0);
0455       hCylindricalAveDiffPhiR[m]->SetStats(0);
0456     }
0457
0458     hCompareRTrue->SetStats(0);
0459     hComparePhiTrue->SetStats(0);
0460
0461     hRDiffvR->SetStats(0);
0462     hRDiffvZ->SetStats(0);
0463     hRDiffvPhi->SetStats(0);
0464
0465     hPhiDiffvR->SetStats(0);
0466     hPhiDiffvZ->SetStats(0);
0467     hPhiDiffvPhi->SetStats(0);
0468
0469     hCompareRTrue_PhiR->SetStats(0);
0470     hComparePhiTrue_PhiR->SetStats(0);
0471
0472     hRDiffvR_PhiR->SetStats(0);
0473     hRDiffvZ_PhiR->SetStats(0);
0474     hRDiffvPhi_PhiR->SetStats(0);
0475
0476     hPhiDiffvR_PhiR->SetStats(0);
0477     hPhiDiffvZ_PhiR->SetStats(0);
0478     hPhiDiffvPhi_PhiR->SetStats(0);
0479
0480     TPad *c1=new TPad("c1","",0.0,0.8,1.0,0.93); //can i do an array of pads?
0481     TPad *c2=new TPad("c2","",0.0,0.64,1.0,0.77);
0482     TPad *c3=new TPad("c3","",0.0,0.48,1.0,0.61);
0483     TPad *c4=new TPad("c4","",0.0,0.32,1.0,0.45);
0484     TPad *c5=new TPad("c5","",0.0,0.16,1.0,0.29);
0485     TPad *c6=new TPad("c6","",0.0,0.0,1.0,0.13);
0486
0487     TPad *titlepad=new TPad("titlepad","",0.0,0.96,1.0,1.0);
0488
0489     TPad *stitlepad1=new TPad("stitlepad1","",0.0,0.93,1.0,0.96);
0490     TPad *stitlepad2=new TPad("stitlepad2","",0.0,0.77,1.0,0.8);
0491     TPad *stitlepad3=new TPad("stitlepad3","",0.0,0.61,1.0,0.64);
0492     TPad *stitlepad4=new TPad("stitlepad4","",0.0,0.45,1.0,0.48);
0493     TPad *stitlepad5=new TPad("stitlepad5","",0.0,0.29,1.0,0.32);
0494     TPad *stitlepad6=new TPad("stitlepad6","",0.0,0.13,1.0,0.16);
0495
0496     TLatex * title = new TLatex(0.0,0.0,"");
0497
0498     TLatex * stitle1 = new TLatex(0.0,0.0,""); //array?
0499     TLatex * stitle2 = new TLatex(0.0,0.0,"");
0500     TLatex * stitle3 = new TLatex(0.0,0.0,"");
0501     TLatex * stitle4 = new TLatex(0.0,0.0,"");
0502     TLatex * stitle5 = new TLatex(0.0,0.0,"");
0503     TLatex * stitle6 = new TLatex(0.0,0.0,"");
0504
0505     title->SetNDC();
0506     stitle1->SetNDC();
0507     stitle2->SetNDC();
0508     stitle3->SetNDC();
0509     stitle4->SetNDC();
0510     stitle5->SetNDC();
0511     stitle6->SetNDC();
0512
0513     title->SetTextSize(0.32);
0514     stitle1->SetTextSize(0.35);
0515     stitle2->SetTextSize(0.35);
0516     stitle3->SetTextSize(0.35);
0517     stitle4->SetTextSize(0.35);
0518     stitle5->SetTextSize(0.35);
0519     stitle6->SetTextSize(0.35);
0520
0521     canvas->cd();
0522     c1->Draw();
0523     stitlepad1->Draw();
0524     c2->Draw();
0525     stitlepad2->Draw();
0526     c3->Draw();
0527     stitlepad3->Draw();
0528     c4->Draw();
0529     stitlepad4->Draw();
0530     c5->Draw();
0531     stitlepad5->Draw();
0532     c6->Draw();
0533     stitlepad6->Draw();
0534     titlepad->Draw();
0535
0536     //x plots
0537     c1->Divide(4,1);
0538     c1->cd(1);
0539     hCartesianAveDiffPhiR[0]->Draw("colz");
0540     c1->cd(2);
0541     hCartesianAveDiffPhiR[1]->Draw("colz");
0542     c1->cd(3);
0543     hCartesianShiftDifferencePhiR[0]->Draw();
0544     //c1->cd(4)->Clear();
0545     c1->cd(4);
0546     //hCMmodelSliceRvTrue->Draw("colz");
0547     hSamplePerBinRZ->Draw("colz");
0548
0549     //y plots
0550     c2->Divide(4,1);
0551     c2->cd(1);
0552     hCartesianAveDiffPhiR[2]->Draw("colz");
0553     c2->cd(2);
0554     hCartesianAveDiffPhiR[3]->Draw("colz");
0555     c2->cd(3);
0556     hCartesianShiftDifferencePhiR[1]->Draw();
0557     //c2->cd(4)->Clear();
0558     c2->cd(4);
0559     //hStripesPerBin->Draw("colz");
0560     hSamplePerBinPhiR->Draw("colz");
0561
0562     //r cart
0563     c3->Divide(4,1);
0564     c3->cd(1);
0565     hCylindricalAveDiffPhiR[0]->Draw("colz");
0566     c3->cd(2);
0567     hCylindricalAveDiffPhiR[1]->Draw("colz");
0568     c3->cd(3);
0569     hCylindricalShiftDifferencePhiR[0]->Draw();
0570     c3->cd(4);
0571     hRShiftTrue->Draw();
0572
0573     //phi cart
0574     c4->Divide(4,1);
0575     c4->cd(1);
0576     hCylindricalAveDiffPhiR[2]->Draw("colz");
0577     c4->cd(2);
0578     hCylindricalAveDiffPhiR[3]->Draw("colz");
0579     c4->cd(3);
0580     hCylindricalShiftDifferencePhiR[1]->Draw();
0581     c4->cd(4);
0582     hPhiShiftTrue->Draw();
0583
0584     //r to true comparison
0585     c5->Divide(4,1);
0586     c5->cd(1);
0587     hCompareRTrue_PhiR->Draw("colz");
0588     c5->cd(2);
0589     hRDiffvR_PhiR->Draw("colz");
0590     c5->cd(3);
0591     hRDiffvZ_PhiR->Draw("colz");
0592     c5->cd(4);
0593     hRDiffvPhi_PhiR->Draw("colz");
0594
0595     //phi to true comparison
0596     c6->Divide(4,1);
0597     c6->cd(1);
0598     hComparePhiTrue_PhiR->Draw("colz");
0599     c6->cd(2);
0600     hPhiDiffvR_PhiR->Draw("colz");
0601     c6->cd(3);
0602     hPhiDiffvZ_PhiR->Draw("colz");
0603     c6->cd(4);
0604     hPhiDiffvPhi_PhiR->Draw("colz");
0605
0606     titlepad->cd();
0607     titlepad->Clear();
0608     title->DrawLatex(0.01,0.4,Form("Event %d; %s", ifile, sourcefilename.Data()));
0609     title->Draw();
0610
0611     stitlepad1->cd();
0612     stitlepad1->Clear();
0613     stitle1->DrawLatex(0.45,0.2,"X Model");
0614     stitle1->Draw();
0615
0616     stitlepad2->cd();
0617     stitlepad2->Clear();
0618     stitle2->DrawLatex(0.45,0.2,"Y Model");
0619     stitle2->Draw();
0620
0621     stitlepad3->cd();
0622     stitlepad3->Clear();
0623     stitle3->DrawLatex(0.45,0.2,"R Model");
0624     stitle3->Draw();
0625
0626     stitlepad4->cd();
0627     stitlepad4->Clear();
0628     stitle4->DrawLatex(0.45,0.2,"Phi Model");
0629     stitle4->Draw();
0630
0631     stitlepad5->cd();
0632     stitlepad5->Clear();
0633     stitle5->DrawLatex(0.4,0.2,"Comparing R Model to True");
0634     stitle5->Draw();
0635
0636     stitlepad6->cd();
0637     stitlepad6->Clear();
0638     stitle6->DrawLatex(0.4,0.2,"Comparing Phi Model to True");
0639     stitle6->Draw();
0640
0641     if(ifile == 0){
0642       //if(ifile == 1){
0643       canvas->Print("CMDistortionAnalysisPhiR.pdf(","pdf");
0644     } else if((ifile == 1) || (ifile == nEvents - 1)){
0645       canvas->Print("CMDistortionAnalysisPhiR.pdf","pdf");
0646     }
0647   }
0648
0649   TCanvas *summary = new TCanvas("summary","ShiftPlotsSummary",2000,3000);
0650
0651   TPad *sumtitlepad = new TPad("sumtitlepad","",0.0,0.96,1.0,1.0);
0652   TPad *sumplots = new TPad("sumplotspad","",0.0,0.0,1.0,0.96);
0653
0654   TLatex *sumtitle = new TLatex(0.0,0.0,"");
0655
0656   sumtitle->SetNDC();
0657   sumtitle->SetTextSize(0.4);
0658
0659   summary->cd();
0660   sumplots->Draw();
0661   sumtitlepad->Draw();
0662
0663   sumplots->Divide(4,6);
0664   sumplots->cd(1);
0665   hDifferenceMeanR->Draw();
0666   sumplots->cd(2);
0667   hDifferenceStdDevR->Draw();
0668   sumplots->cd(3);
0669   hTrueMeanR->Draw();
0670   sumplots->cd(4);
0671   hTrueStdDevR->Draw();
0672   sumplots->cd(5);
0673   hDifferenceMeanPhi->Draw();
0674   sumplots->cd(6);
0675   hDifferenceStdDevPhi->Draw();
0676   sumplots->cd(7);
0677   hTrueMeanPhi->Draw();
0678   sumplots->cd(8);
0679   hTrueStdDevPhi->Draw();
0680   sumplots->cd(9);
0681   sumplots->cd(10)->Clear();
0682   sumplots->cd(11)->Clear();
0683   sumplots->cd(12)->Clear();
0684   sumplots->cd(13)->Clear();
0685   sumplots->cd(14)->Clear();
0686   sumplots->cd(15)->Clear();
0687   sumplots->cd(16)->Clear();
0688   sumplots->cd(17)->Clear();
0689   sumplots->cd(18)->Clear();
0690   sumplots->cd(19)->Clear();
0691   sumplots->cd(20)->Clear();
0692   sumplots->cd(21)->Clear();
0693   sumplots->cd(22)->Clear();
0694   sumplots->cd(23)->Clear();
0695   sumplots->cd(24)->Clear();
0696
0697   sumtitlepad->cd();
0698   sumtitlepad->Clear();
0699   sumtitle->DrawLatex(0.4,0.4,"Summary of Events");
0700   summary->Print("CMDistortionAnalysisPhiR.pdf)","pdf");
0701
0702   return 0;
0703 }