20 const char * inputpattern=
"/sphenix/u/skurdi/CMCalibration/cmDistHitsTree_Event*.root";
23 TFileCollection *filelist=
new TFileCollection();
24 filelist->Add(inputpattern);
25 TString sourcefilename;
29 nEvents=filelist->GetNFiles();
30 }
else if(nMaxEvents<filelist->GetNFiles()){
33 nEvents= filelist->GetNFiles();
37 TCanvas *canvas1=
new TCanvas(
"canvas1",
"CMDistortionReco1",1200,800);
41 TCanvas *canvas=
new TCanvas(
"canvas",
"CMDistortionReco2",400,400);
44 position =
new TVector3(1.,1.,1.);
45 newposition =
new TVector3(1.,1.,1.);
48 TH1F *hTimePerEvent =
new TH1F(
"hTimePerEvent",
"Time Per Event; time (ms)",20,0,10000);
50 for (
int ifile=0;ifile <
nEvents;ifile++){
54 unsigned long before = now;
57 sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
59 char const *treename=
"cmDistHitsTree";
60 TFile *
input=TFile::Open(sourcefilename,
"READ");
61 TTree *inTree=(TTree*)input->Get(
"tree");
63 inTree->SetBranchAddress(
"position",&position);
64 inTree->SetBranchAddress(
"newposition",&newposition);
68 double lowphi = -0.078539819;
69 double highphi = 6.3617253;
76 TH2F *hStripesPerBin =
new TH2F(
"hStripesPerBin",
"CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
79 TH2F *hStripesPerBinPhiR =
new TH2F(
"hStripesPerBinPhiR",
"CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
81 TH2F *hCartesianForward[3];
82 hCartesianForward[0] =
new TH2F(
"hForwardX",
"X Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
83 hCartesianForward[1] =
new TH2F(
"hForwardY",
"Y Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
84 hCartesianForward[2] =
new TH2F(
"hForwardZ",
"Z Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
87 TH2F *hCartesianForwardPhiR[3];
88 hCartesianForwardPhiR[0] =
new TH2F(
"hForwardX_PhiR",
"X Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
89 hCartesianForwardPhiR[1] =
new TH2F(
"hForwardY_PhiR",
"Y Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
90 hCartesianForwardPhiR[2] =
new TH2F(
"hForwardZ_PhiR",
"Z Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
92 TH2F *hCylindricalForwardPhiR[2];
93 hCylindricalForwardPhiR[0] =
new TH2F(
"hForwardR_PhiR",
"Radial Shift Forward of Stripe Centers, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
94 hCylindricalForwardPhiR[1] =
new TH2F(
"hForwardPhi_PhiR",
"Phi Shift Forward of Stripe Centers, Phi,R binning (rad); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
96 for (
int i=0;
i<inTree->GetEntries();
i++){
99 hStripesPerBin->Fill(position->X(),position->Y(),1);
102 double r = position->Perp();
103 double phi = position->Phi();
105 if(position->Phi() < 0.0){
106 phi = position->Phi() + 2.0*TMath::Pi();
109 hStripesPerBinPhiR->Fill(phi,r,1);
111 deltaX = (newposition->X() - position->X())*(1e4);
112 deltaY = (newposition->Y() - position->Y())*(1e4);
113 deltaZ = (newposition->Z() - position->Z())*(1e4);
115 deltaR = (newposition->Perp() - position->Perp())*(1e4);
116 deltaPhi = newposition->DeltaPhi(*position);
118 hCartesianForward[0]->Fill(position->X(),position->Y(),deltaX);
119 hCartesianForward[1]->Fill(position->X(),position->Y(),deltaY);
120 hCartesianForward[2]->Fill(position->X(),position->Y(),deltaZ);
123 hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
124 hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
125 hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
127 hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
128 hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
131 TH2F *hCartesianAveShift[3];
132 hCartesianAveShift[0] =
new TH2F(
"AveShiftX",
"Average of CM Model X over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
133 hCartesianAveShift[1] =
new TH2F(
"AveShiftY",
"Average of CM Model Y over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
134 hCartesianAveShift[2] =
new TH2F(
"AveShiftZ",
"Average of CM Model Z over Stripes per Bin (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
136 TH2F *hCylindricalAveShift[2];
137 hCylindricalAveShift[0] =
new TH2F(
"AveShiftRCart",
"Average of CM Model R over Stripes per Bin from Cartesian (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
138 hCylindricalAveShift[1] =
new TH2F(
"AveShiftPhiCart",
"Average of CM Model Phi over Stripes per Bin from Cartesian (rad); x (cm); y (cm)",nbins,low,high,nbins,low,high);
141 TH2F *hCartesianAveShiftPhiR[3];
142 hCartesianAveShiftPhiR[0] =
new TH2F(
"AveShiftX_PhiR",
"Average of CM Model X over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
143 hCartesianAveShiftPhiR[1] =
new TH2F(
"AveShiftY_PhiR",
"Average of CM Model Y over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
144 hCartesianAveShiftPhiR[2] =
new TH2F(
"AveShiftZ_PhiR",
"Average of CM Model Z over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
146 TH2F *hCylindricalAveShiftPhiR[2];
147 hCylindricalAveShiftPhiR[0] =
new TH2F(
"AveShiftR_PhiR",
"Average of CM Model R over Stripes per Bin, Phi,R binning (#mum); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
148 hCylindricalAveShiftPhiR[1] =
new TH2F(
"AveShiftPhi_PhiR",
"Average of CM Model Phi over Stripes per Bin, Phi,R binning (rad); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
150 for (
int i = 0;
i < 3;
i ++){
151 hCartesianAveShift[
i]->Divide(hCartesianForward[
i],hStripesPerBin);
152 hCartesianAveShiftPhiR[
i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
155 hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
156 hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
160 double x = low + ((high - low)/(1.0*nbins))*(
i+0.5);
163 double y = low + ((high - low)/(1.0*nbins))*(
j+0.5);
165 int xbin = hCartesianAveShift[0]->FindBin(x,y);
166 int ybin = hCartesianAveShift[1]->FindBin(x,y);
167 double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1
e-4);
168 double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1
e-4);
170 TVector3 shifted, original;
173 shifted.SetX(x+xaveshift);
174 shifted.SetY(y+yaveshift);
178 double raveshift = (shifted.Perp() - original.Perp())*(1e4);
179 double phiaveshift = shifted.DeltaPhi(original);
182 hCylindricalAveShift[0]->Fill(x,y,raveshift);
183 hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
193 double minphi = -0.078539819;
194 double minr = 18.884615;
196 double minz = -1.3187500;
198 double maxphi = 6.3617253;
199 double maxr = 79.115387;
200 double maxz = 106.81875;
202 TH3F *hCartesianCMModel[3];
203 hCartesianCMModel[0]=
new TH3F(
"hCMModelX",
"CM Model: X Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
204 hCartesianCMModel[1]=
new TH3F(
"hCMModelY",
"CM Model: Y Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
205 hCartesianCMModel[2]=
new TH3F(
"hCMModelZ",
"CM Model: Z Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
207 TH3F *hCylindricalCMModel[2];
208 hCylindricalCMModel[0]=
new TH3F(
"hCMModelRCart",
"CM Model: Radial Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
209 hCylindricalCMModel[1]=
new TH3F(
"hCMModelPhiCart",
"CM Model: Phi Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
212 TH3F *hCartesianCMModelPhiR[3];
213 hCartesianCMModelPhiR[0]=
new TH3F(
"hCMModelX_PhiR",
"CM Model: X Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
214 hCartesianCMModelPhiR[1]=
new TH3F(
"hCMModelY_PhiR",
"CM Model: Y Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
215 hCartesianCMModelPhiR[2]=
new TH3F(
"hCMModelZ_PhiR",
"CM Model: Z Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
217 TH3F *hCylindricalCMModelPhiR[2];
218 hCylindricalCMModelPhiR[0]=
new TH3F(
"hCMModelR_PhiR",
"CM Model: Radial Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
219 hCylindricalCMModelPhiR[1]=
new TH3F(
"hCMModelPhi_PhiR",
"CM Model: Phi Shift Forward of Stripe Centers, Phi,R binning", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
221 double xshift, yshift, zshift, rshiftcart, phishiftcart;
222 double xshiftPhiR, yshiftPhiR, zshiftPhiR, rshiftcartPhiR, phishiftcartPhiR;
225 double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(
i+0.5);
227 for(
int j = 0;
j < nr;
j++){
228 double r = minr + ((maxr - minr)/(1.0*nr))*(
j+0.5);
230 double x = r*cos(phi);
231 double y = r*sin(phi);
233 for(
int k = 0;
k <
nz;
k++){
234 double z = minz + ((maxz - minz)/(1.0*nz))*(
k+0.5);
236 xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1
e-4);
237 yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1
e-4);
238 zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1
e-4);
240 rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1
e-4);
241 phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
243 hCartesianCMModel[0]->Fill(phi,r,z,xshift*(1-z/105.5));
244 hCartesianCMModel[1]->Fill(phi,r,z,yshift*(1-z/105.5));
245 hCartesianCMModel[2]->Fill(phi,r,z,zshift*(1-z/105.5));
247 hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
248 hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
251 xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1
e-4);
252 yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1
e-4);
253 zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1
e-4);
255 rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1
e-4);
256 phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
258 hCartesianCMModelPhiR[0]->Fill(phi,r,z,xshiftPhiR*(1-z/105.5));
259 hCartesianCMModelPhiR[1]->Fill(phi,r,z,yshiftPhiR*(1-z/105.5));
260 hCartesianCMModelPhiR[2]->Fill(phi,r,z,zshiftPhiR*(1-z/105.5));
262 hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
263 hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
270 plots=TFile::Open(Form(
"CMModels_Event%d.root",ifile),
"RECREATE");
271 hStripesPerBin->Write();
273 for(
int i = 0;
i < 3;
i++){
274 hCartesianForward[
i]->Write();
275 hCartesianAveShift[
i]->Write();
276 hCartesianCMModel[
i]->Write();
278 hCartesianForwardPhiR[
i]->Write();
279 hCartesianAveShiftPhiR[
i]->Write();
280 hCartesianCMModelPhiR[
i]->Write();
283 for(
int i = 0;
i < 2;
i++){
284 hCylindricalAveShift[
i]->Write();
285 hCylindricalCMModel[
i]->Write();
287 hCylindricalForwardPhiR[
i]->Write();
288 hCylindricalAveShiftPhiR[
i]->Write();
289 hCylindricalCMModelPhiR[
i]->Write();
297 unsigned long after = now;
299 hTimePerEvent->Fill(after-before);
302 for (
int i = 0;
i < 3;
i++){
303 hCartesianForward[
i]->SetStats(0);
304 hCartesianForwardPhiR[
i]->SetStats(0);
307 hStripesPerBin->SetStats(0);
309 hStripesPerBinPhiR->Draw(
"colz");
311 hCartesianForward[1]->Draw(
"colz");
313 hCartesianForward[2]->Draw(
"colz");
315 hCartesianForwardPhiR[0]->Draw(
"colz");
317 hCartesianForwardPhiR[1]->Draw(
"colz");
319 hCartesianForwardPhiR[2]->Draw(
"colz");
322 canvas1->Print(
"CMDistortionReco1.pdf(",
"pdf");
323 }
else if (ifile == nEvents - 1){
324 canvas1->Print(
"CMDistortionReco1.pdf)",
"pdf");
326 canvas1->Print(
"CMDistortionReco1.pdf",
"pdf");
332 hTimePerEvent->Draw();
333 canvas->Print(
"CMDistortionReco2.pdf",
"pdf");