Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CMDistortionRecoPhiR.C
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CMDistortionRecoPhiR.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 CMDistortionRecoPhiR(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  TCanvas *canvas1=new TCanvas("canvas1","CMDistortionReco1",1200,800);
37  canvas1->Divide(3,2);
38 
39  //canvas for time plot
40  TCanvas *canvas=new TCanvas("canvas","CMDistortionReco2",400,400);
41 
42  TVector3 *position, *newposition;
43  position = new TVector3(1.,1.,1.);
44  newposition = new TVector3(1.,1.,1.);
45 
46  //histogram to compare times
47  TH1F *hTimePerEvent = new TH1F("hTimePerEvent","Time Per Event; time (ms)",20,0,10000);
48 
49  for (int ifile=0;ifile < nEvents;ifile++){
50  //call to TTime before opening ttree
51  TTime now;
52  now=gSystem->Now();
53  unsigned long before = now;
54 
55  //get data from ttree
56  sourcefilename=((TFileInfo*)(filelist->GetList()->At(ifile)))->GetCurrentUrl()->GetFile();
57 
58  char const *treename="cmDistHitsTree";
59  TFile *input=TFile::Open(sourcefilename, "READ");
60  TTree *inTree=(TTree*)input->Get("tree");
61 
62  inTree->SetBranchAddress("position",&position);
63  inTree->SetBranchAddress("newposition",&newposition);
64 
65  //hardcoded numbers from average distortion file's hIntDistortionPosX
66  int nbinsphi = 30; //when using 35, blank spots at around r = 22 cm, phi just above n below pi
67  double lowphi = -0.078539819;
68  double highphi = 6.3617253;
69  int nbinsr = 30; // when using 35, blank stripe around r = 58 cm
70  double lowr = 0.0;
71  double highr = 90.0;
72 
73  //for forward only
74 
75  TH2F *hStripesPerBinPhiR = new TH2F("hStripesPerBinPhiR","CM Stripes Per Bin (z in stripes); phi (rad); r (cm)",nbinsphi,lowphi,highphi,nbinsr,lowr,highr);
76 
77  TH2F *hCartesianForwardPhiR[3];
78  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);
79  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);
80  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);
81 
82  TH2F *hCylindricalForwardPhiR[2];
83  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);
84  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);
85 
86  for (int i=0;i<inTree->GetEntries();i++){
87  inTree->GetEntry(i);
88 
89  double r = position->Perp();
90  double phi = position->Phi();
91 
92  if(position->Phi() < 0.0){
93  phi = position->Phi() + 2.0*TMath::Pi();
94  }
95 
96  hStripesPerBinPhiR->Fill(phi,r,1);
97 
98  deltaX = (newposition->X() - position->X())*(1e4); //convert from cm to micron
99  deltaY = (newposition->Y() - position->Y())*(1e4);
100  deltaZ = (newposition->Z() - position->Z())*(1e4);
101 
102  deltaR = (newposition->Perp() - position->Perp())*(1e4);
103  deltaPhi = newposition->DeltaPhi(*position);
104 
105  hCartesianForwardPhiR[0]->Fill(phi,r,deltaX);
106  hCartesianForwardPhiR[1]->Fill(phi,r,deltaY);
107  hCartesianForwardPhiR[2]->Fill(phi,r,deltaZ);
108 
109  hCylindricalForwardPhiR[0]->Fill(phi,r,deltaR);
110  hCylindricalForwardPhiR[1]->Fill(phi,r,deltaPhi);
111  }
112 
113  TH2F *hCartesianAveShiftPhiR[3];
114  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);
115  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);
116  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);
117 
118  TH2F *hCylindricalAveShiftPhiR[2];
119  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);
120  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);
121 
122  for (int i = 0; i < 3; i ++){
123  hCartesianAveShiftPhiR[i]->Divide(hCartesianForwardPhiR[i],hStripesPerBinPhiR);
124  }
125 
126  hCylindricalAveShiftPhiR[0]->Divide(hCylindricalForwardPhiR[0],hStripesPerBinPhiR);
127  hCylindricalAveShiftPhiR[1]->Divide(hCylindricalForwardPhiR[1],hStripesPerBinPhiR);
128 
129  //same range and bins for each coordinate, binned in cm
130  //hardcoded numbers from average distortion file's hIntDistortionPosX
131  int nphi = 82;
132  int nr = 54;
133  int nz = 82;
134 
135  double minphi = -0.078539819;
136  double minr = 18.884615;
137  double minz = -1.3187500;
138 
139  double maxphi = 6.3617253;
140  double maxr = 79.115387;
141  double maxz = 106.81875;
142 
143  TH3F *hCartesianCMModelPhiR[3];
144  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
145  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);
146  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);
147 
148  TH3F *hCylindricalCMModelPhiR[2];
149  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);
150  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);
151 
152  double xshiftPhiR, yshiftPhiR, zshiftPhiR, rshiftcartPhiR, phishiftcartPhiR;
153 
154  for(int i = 0; i < nphi; i++){
155  double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(i+0.5); //center of bin
156 
157  for(int j = 0; j < nr; j++){
158  double r = minr + ((maxr - minr)/(1.0*nr))*(j+0.5); //center of bin
159 
160  for(int k = 0; k < nz; k++){
161  double z = minz + ((maxz - minz)/(1.0*nz))*(k+0.5); //center of bin
162 
163  xshiftPhiR=(hCartesianAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);//coordinate of your stripe
164  yshiftPhiR=(hCartesianAveShiftPhiR[1]->Interpolate(phi,r))*(1e-4);//convert micron to cm
165  zshiftPhiR=(hCartesianAveShiftPhiR[2]->Interpolate(phi,r))*(1e-4);
166 
167  rshiftcartPhiR=(hCylindricalAveShiftPhiR[0]->Interpolate(phi,r))*(1e-4);
168  phishiftcartPhiR=hCylindricalAveShiftPhiR[1]->Interpolate(phi,r);
169 
170  hCartesianCMModelPhiR[0]->Fill(phi,r,z,xshiftPhiR*(1-z/105.5));
171  hCartesianCMModelPhiR[1]->Fill(phi,r,z,yshiftPhiR*(1-z/105.5));
172  hCartesianCMModelPhiR[2]->Fill(phi,r,z,zshiftPhiR*(1-z/105.5));
173 
174  hCylindricalCMModelPhiR[0]->Fill(phi,r,z,rshiftcartPhiR*(1-z/105.5));
175  hCylindricalCMModelPhiR[1]->Fill(phi,r,z,phishiftcartPhiR*(1-z/105.5));
176  }
177  }
178  }
179 
180  TFile *plots;
181 
182  plots=TFile::Open(Form("CMModelsPhiR_Event%d.root",ifile),"RECREATE");
183 
184  for(int i = 0; i < 3; i++){
185  hCartesianForwardPhiR[i]->Write();
186  hCartesianAveShiftPhiR[i]->Write();
187  hCartesianCMModelPhiR[i]->Write();
188  }
189 
190  for(int i = 0; i < 2; i++){
191  hCylindricalForwardPhiR[i]->Write();
192  hCylindricalAveShiftPhiR[i]->Write();
193  hCylindricalCMModelPhiR[i]->Write();
194  }
195 
196  plots->Close();
197 
198 
199  //call to TTime after outputting TH3Fs
200  now=gSystem->Now();
201  unsigned long after = now;
202 
203  hTimePerEvent->Fill(after-before);
204 
205  //to check histograms
206  for (int i = 0; i < 3; i++){
207  hCartesianForwardPhiR[i]->SetStats(0);
208  }
209 
210  canvas1->cd(1);
211  hStripesPerBinPhiR->Draw("colz");
212  canvas1->cd(2)->Clear();
213  canvas1->cd(3)->Clear();
214  canvas1->cd(4);
215  hCartesianForwardPhiR[0]->Draw("colz");
216  canvas1->cd(5);
217  hCartesianForwardPhiR[1]->Draw("colz");
218  canvas1->cd(6);
219  hCartesianForwardPhiR[2]->Draw("colz");
220 
221  if(ifile == 0){
222  canvas1->Print("CMDistortionReco1.pdf(","pdf");
223  } else if (ifile == nEvents - 1){
224  canvas1->Print("CMDistortionReco1.pdf)","pdf");
225  } else {
226  canvas1->Print("CMDistortionReco1.pdf","pdf");
227  }
228 
229  }
230 
231  canvas->cd();
232  hTimePerEvent->Draw();
233  canvas->Print("CMDistortionReco2.pdf","pdf");
234 
235  return 0;
236 }