ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CMDistortionRecoCart.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CMDistortionRecoCart.C
1 // step 2
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 CMDistortionRecoCart(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 
36  //canvas for checking data
37  //TCanvas *canvas1=new TCanvas("canvas1","CMDistortionRecoCart1",1200,800);
38  //canvas1->Divide(3,2);
39 
40  //canvas for time plot
41  TCanvas *canvas=new TCanvas("canvas","CMDistortionRecoCart2",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  //for forward only
67 
68  TH2F *hStripesPerBin = new TH2F("hStripesPerBin","CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
69 
70  TH2F *hCartesianForward[3];
71  hCartesianForward[0] = new TH2F("hForwardX","X Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
72  hCartesianForward[1] = new TH2F("hForwardY","Y Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
73  hCartesianForward[2] = new TH2F("hForwardZ","Z Shift Forward of Stripe Centers (#mum); x (cm); y (cm)",nbins,low,high,nbins,low,high);
74 
75  for (int i=0;i<inTree->GetEntries();i++){
76  inTree->GetEntry(i);
77 
78  hStripesPerBin->Fill(position->X(),position->Y(),1);
79 
80  deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
81  deltaY = (newposition->Y() - position->Y())*(1e4);
82  deltaZ = (newposition->Z() - position->Z())*(1e4);
83 
84  deltaR = (newposition->Perp() - position->Perp())*(1e4);
85  deltaPhi = newposition->DeltaPhi(*position);
86 
87  hCartesianForward[0]->Fill(position->X(),position->Y(),deltaX);
88  hCartesianForward[1]->Fill(position->X(),position->Y(),deltaY);
89  hCartesianForward[2]->Fill(position->X(),position->Y(),deltaZ);
90  }
91 
92  TH2F *hCartesianAveShift[3];
93  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);
94  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);
95  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);
96 
97  TH2F *hCylindricalAveShift[2];
98  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);
99  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);
100 
101  for (int i = 0; i < 3; i ++){
102  hCartesianAveShift[i]->Divide(hCartesianForward[i],hStripesPerBin);
103  }
104 
105  //cyl models from cart coords
106  for(int i = 0; i < nbins; i++){
107  double x = low + ((high - low)/(1.0*nbins))*(i+0.5); //center of bin
108 
109  for(int j = 0; j < nbins; j++){
110  double y = low + ((high - low)/(1.0*nbins))*(j+0.5); //center of bin
111 
112  int xbin = hCartesianAveShift[0]->FindBin(x,y);
113  int ybin = hCartesianAveShift[1]->FindBin(x,y);
114  double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1e-4); // converts microns to cm
115  double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1e-4);
116 
117  TVector3 shifted, original;
118  original.SetX(x);
119  original.SetY(y);
120  shifted.SetX(x+xaveshift);
121  shifted.SetY(y+yaveshift);
122  //x n y above for orig
123  //shifted is orig + ave shift
124 
125  double raveshift = (shifted.Perp() - original.Perp())*(1e4);
126  double phiaveshift = shifted.DeltaPhi(original);
127 
128  //fill with r from x n y
129  hCylindricalAveShift[0]->Fill(x,y,raveshift);
130  hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
131  }
132  }
133 
134  //same range and bins for each coordinate, binned in cm
135  //hardcoded numbers from average distortion file's hIntDistortionPosX
136  int nphi = 82;
137  int nr = 54;
138  int nz = 82;
139 
140  double minphi = -0.078539819;
141  double minr = 18.884615;
142  //double minz = 5.0;
143  double minz = -1.3187500;
144 
145  double maxphi = 6.3617253;
146  double maxr = 79.115387;
147  double maxz = 106.81875;
148 
149  TH3F *hCartesianCMModel[3];
150  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
151  hCartesianCMModel[1]=new TH3F("hCMModelY", "CM Model: Y Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
152  hCartesianCMModel[2]=new TH3F("hCMModelZ", "CM Model: Z Shift Forward of Stripe Centers", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
153 
154  TH3F *hCylindricalCMModel[2];
155  hCylindricalCMModel[0]=new TH3F("hCMModelRCart", "CM Model: Radial Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
156  hCylindricalCMModel[1]=new TH3F("hCMModelPhiCart", "CM Model: Phi Shift Forward of Stripe Centers from Cartesian", nphi,minphi,maxphi, nr,minr,maxr, nz,minz,maxz);
157 
158  double xshift, yshift, zshift, rshiftcart, phishiftcart;
159 
160  for(int i = 0; i < nphi; i++){
161  double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
162 
163  for(int j = 0; j < nr; j++){
164  double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
165 
166  double x = r*cos(phi); //cm
167  double y = r*sin(phi);
168 
169  for(int k = 0; k < nz; k++){
170  double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
171 
172  xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1e-4);//coordinate of your stripe
173  yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1e-4);//convert micron to cm
174  zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1e-4);
175 
176  rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1e-4);
177  phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
178 
179  hCartesianCMModel[0]->Fill(phi,r,z,xshift*(1-z/105.5));
180  hCartesianCMModel[1]->Fill(phi,r,z,yshift*(1-z/105.5));
181  hCartesianCMModel[2]->Fill(phi,r,z,zshift*(1-z/105.5));
182 
183  hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
184  hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
185  }
186  }
187  }
188 
189  TFile *plots;
190 
191  plots=TFile::Open(Form("CMModelsCart_Event%d.root",ifile),"RECREATE");
192  hStripesPerBin->Write();
193 
194  for(int i = 0; i < 3; i++){
195  hCartesianForward[i]->Write();
196  hCartesianAveShift[i]->Write();
197  hCartesianCMModel[i]->Write();
198  }
199 
200  for(int i = 0; i < 2; i++){
201  hCylindricalAveShift[i]->Write();
202  hCylindricalCMModel[i]->Write();
203  }
204 
205  plots->Close();
206 
207 
208  //call to TTime after outputting TH3Fs
209  now=gSystem->Now();
210  unsigned long after = now;
211 
212  hTimePerEvent->Fill(after-before);
213 
214  /*
215  //to check histograms
216  for (int i = 0; i < 3; i++){
217  hCartesianForward[i]->SetStats(0);
218  }
219 
220  hStripesPerBin->SetStats(0);
221  canvas1->cd(1);
222  hCartesianForward[0]->Draw("colz");
223  canvas1->cd(2);
224  hCartesianForward[1]->Draw("colz");
225  canvas1->cd(3);
226  hCartesianForward[2]->Draw("colz");
227  canvas1->cd(4)->Clear();
228  canvas1->cd(5)->Clear();
229  canvas1->cd(6)->Clear();
230 
231  if(ifile == 0){
232  canvas1->Print("CMDistortionRecoCart1.pdf(","pdf");
233  } else if (ifile == nEvents - 1){
234  canvas1->Print("CMDistortionRecoCart1.pdf)","pdf");
235  } else {
236  canvas1->Print("CMDistortionRecoCart1.pdf","pdf");
237  }*/
238 
239 
240  }
241 
242  canvas->cd();
243  hTimePerEvent->Draw();
244  canvas->Print("CMDistortionRecoCart2.pdf","pdf");
245 
246  return 0;
247 }