ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CMDistortionReco.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CMDistortionReco.C
1 // step 2 with phi,r coords
2 #include <iostream>
3 #include <cmath>
4 #include <vector>
5 #include "TMath.h"
6 #include "TVector3.h"
7 #include "TTree.h"
8 #include <TTime.h>
9 
10 using namespace std;
11 
12 int CMDistortionReco(int nMaxEvents = -1) {
13  int nbins = 35;
14  double low = -80.0;
15  double high = 80.0;
16  double deltaX, deltaY, deltaZ, deltaR, deltaPhi;
17  int nEvents;
18 
19  //take in events
20  const char * inputpattern="/sphenix/u/skurdi/CMCalibration/cmDistHitsTree_Event*.root";
21 
22  //find all files that match the input string (includes wildcards)
23  TFileCollection *filelist=new TFileCollection();
24  filelist->Add(inputpattern);
25  TString sourcefilename;
26 
27  //how many events
28  if (nMaxEvents<0){
29  nEvents=filelist->GetNFiles();
30  } else if(nMaxEvents<filelist->GetNFiles()){
31  nEvents=nMaxEvents;
32  } else {
33  nEvents= filelist->GetNFiles();
34  }
35  // nEvents = 2;
36 
37  TCanvas *canvas1=new TCanvas("canvas1","CMDistortionReco1",1200,800);
38  canvas1->Divide(3,2);
39 
40  //canvas for time plot
41  TCanvas *canvas=new TCanvas("canvas","CMDistortionReco2",400,400);
42 
43  TVector3 *position, *newposition;
44  position = new TVector3(1.,1.,1.);
45  newposition = new TVector3(1.,1.,1.);
46 
47  //histogram to compare times
48  TH1F *hTimePerEvent = new TH1F("hTimePerEvent","Time Per Event; time (ms)",20,0,10000);
49 
50  for (int ifile=0;ifile < nEvents;ifile++){
51  //call to TTime before opening ttree
52  TTime now;
53  now=gSystem->Now();
54  unsigned long before = now;
55 
56  //get data from ttree
57  sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
58 
59  char const *treename="cmDistHitsTree";
60  TFile *input=TFile::Open(sourcefilename, "READ");
61  TTree *inTree=(TTree*)input->Get("tree");
62 
63  inTree->SetBranchAddress("position",&position);
64  inTree->SetBranchAddress("newposition",&newposition);
65 
66  //hardcoded numbers from average distortion file's hIntDistortionPosX
67  int nbinsphi = 30; //when using 35, blank spots at around r = 22 cm, phi just above n below pi
68  double lowphi = -0.078539819;
69  double highphi = 6.3617253;
70  int nbinsr = 30; // when using 35, blank stripe around r = 58 cm
71  double lowr = 0.0;
72  double highr = 90.0;
73 
74  //for forward only
75 
76  TH2F *hStripesPerBin = new TH2F("hStripesPerBin","CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
77 
78  //phi,r binning
79  TH2F *hStripesPerBinPhiR = new TH2F("hStripesPerBinPhiR","CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
80 
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);
85 
86  //phi,r binning
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);
91 
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);
95 
96  for (int i=0;i<inTree->GetEntries();i++){
97  inTree->GetEntry(i);
98 
99  hStripesPerBin->Fill(position->X(),position->Y(),1);
100 
101  //for phi,r binning
102  double r = position->Perp();
103  double phi = position->Phi();
104 
105  if(position->Phi() < 0.0){
106  phi = position->Phi() + 2.0*TMath::Pi();
107  }
108 
109  hStripesPerBinPhiR->Fill(phi,r,1);
110 
111  deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
112  deltaY = (newposition->Y() - position->Y())*(1e4);
113  deltaZ = (newposition->Z() - position->Z())*(1e4);
114 
115  deltaR = (newposition->Perp() - position->Perp())*(1e4);
116  deltaPhi = newposition->DeltaPhi(*position);
117 
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);
121 
122  // phi,r binning
123  hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
124  hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
125  hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
126 
127  hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
128  hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
129  }
130 
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);
135 
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);
139 
140  // phi,r binning
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);
145 
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);
149 
150  for (int i = 0; i < 3; i ++){
151  hCartesianAveShift[i]->Divide(hCartesianForward[i],hStripesPerBin);
152  hCartesianAveShiftPhiR[i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
153  }
154 
155  hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
156  hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
157 
158  //cyl models from cart coords
159  for(int i = 0; i < nbins; i++){
160  double x = low + ((high - low)/(1.0*nbins))*(i+0.5); //center of bin
161 
162  for(int j = 0; j < nbins; j++){
163  double y = low + ((high - low)/(1.0*nbins))*(j+0.5); //center of bin
164 
165  int xbin = hCartesianAveShift[0]->FindBin(x,y);
166  int ybin = hCartesianAveShift[1]->FindBin(x,y);
167  double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1e-4); // converts microns to cm
168  double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1e-4);
169 
170  TVector3 shifted, original;
171  original.SetX(x);
172  original.SetY(y);
173  shifted.SetX(x+xaveshift);
174  shifted.SetY(y+yaveshift);
175  //x n y above for orig
176  //shifted is orig + ave shift
177 
178  double raveshift = (shifted.Perp() - original.Perp())*(1e4);
179  double phiaveshift = shifted.DeltaPhi(original);
180 
181  //fill with r from x n y
182  hCylindricalAveShift[0]->Fill(x,y,raveshift);
183  hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
184  }
185  }
186 
187  //same range and bins for each coordinate, binned in cm
188  //hardcoded numbers from average distortion file's hIntDistortionPosX
189  int nphi = 82;
190  int nr = 54;
191  int nz = 82;
192 
193  double minphi = -0.078539819;
194  double minr = 18.884615;
195  //double minz = 5.0;
196  double minz = -1.3187500;
197 
198  double maxphi = 6.3617253;
199  double maxr = 79.115387;
200  double maxz = 106.81875;
201 
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); //rad, cm, cm
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);
206 
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);
210 
211  //phi,r binning
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); //rad, cm, cm
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);
216 
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);
220 
221  double xshift, yshift, zshift, rshiftcart, phishiftcart;
222  double xshiftPhiR, yshiftPhiR, zshiftPhiR, rshiftcartPhiR, phishiftcartPhiR;
223 
224  for(int i = 0; i < nphi; i++){
225  double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
226 
227  for(int j = 0; j < nr; j++){
228  double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
229 
230  double x = r*cos(phi); //cm
231  double y = r*sin(phi);
232 
233  for(int k = 0; k < nz; k++){
234  double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
235 
236  xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1e-4);//coordinate of your stripe
237  yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1e-4);//convert micron to cm
238  zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1e-4);
239 
240  rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1e-4);
241  phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
242 
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));
246 
247  hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
248  hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
249 
250  //phi,r binning
251  xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);//coordinate of your stripe
252  yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1e-4);//convert micron to cm
253  zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1e-4);
254 
255  rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);
256  phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
257 
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));
261 
262  hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
263  hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
264  }
265  }
266  }
267 
268  TFile *plots;
269 
270  plots=TFile::Open(Form("CMModels_Event%d.root",ifile),"RECREATE");
271  hStripesPerBin->Write();
272 
273  for(int i = 0; i < 3; i++){
274  hCartesianForward[i]->Write();
275  hCartesianAveShift[i]->Write();
276  hCartesianCMModel[i]->Write();
277 
278  hCartesianForwardPhiR[i]->Write();
279  hCartesianAveShiftPhiR[i]->Write();
280  hCartesianCMModelPhiR[i]->Write();
281  }
282 
283  for(int i = 0; i < 2; i++){
284  hCylindricalAveShift[i]->Write();
285  hCylindricalCMModel[i]->Write();
286 
287  hCylindricalForwardPhiR[i]->Write();
288  hCylindricalAveShiftPhiR[i]->Write();
289  hCylindricalCMModelPhiR[i]->Write();
290  }
291 
292  plots->Close();
293 
294 
295  //call to TTime after outputting TH3Fs
296  now=gSystem->Now();
297  unsigned long after = now;
298 
299  hTimePerEvent->Fill(after-before);
300 
301  //to check histograms
302  for (int i = 0; i < 3; i++){
303  hCartesianForward[i]->SetStats(0);
304  hCartesianForwardPhiR[i]->SetStats(0);
305  }
306 
307  hStripesPerBin->SetStats(0);
308  canvas1->cd(1);
309  hStripesPerBinPhiR->Draw("colz");
310  canvas1->cd(2);
311  hCartesianForward[1]->Draw("colz");
312  canvas1->cd(3);
313  hCartesianForward[2]->Draw("colz");
314  canvas1->cd(4);
315  hCartesianForwardPhiR[0]->Draw("colz");
316  canvas1->cd(5);
317  hCartesianForwardPhiR[1]->Draw("colz");
318  canvas1->cd(6);
319  hCartesianForwardPhiR[2]->Draw("colz");
320 
321  if(ifile == 0){
322  canvas1->Print("CMDistortionReco1.pdf(","pdf");
323  } else if (ifile == nEvents - 1){
324  canvas1->Print("CMDistortionReco1.pdf)","pdf");
325  } else {
326  canvas1->Print("CMDistortionReco1.pdf","pdf");
327  }
328 
329  }
330 
331  canvas->cd();
332  hTimePerEvent->Draw();
333  canvas->Print("CMDistortionReco2.pdf","pdf");
334 
335  return 0;
336 }