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();
41 TCanvas *canvas=
new TCanvas(
"canvas",
"CMDistortionRecoCart2",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 TH2F *hStripesPerBin =
new TH2F(
"hStripesPerBin",
"CM Stripes Per Bin (z in stripes); x (cm); y (cm)",nbins,low,high,nbins,low,high);
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);
75 for (
int i=0;
i<inTree->GetEntries();
i++){
78 hStripesPerBin->Fill(position->X(),position->Y(),1);
80 deltaX = (newposition->X() - position->X())*(1e4);
81 deltaY = (newposition->Y() - position->Y())*(1e4);
82 deltaZ = (newposition->Z() - position->Z())*(1e4);
84 deltaR = (newposition->Perp() - position->Perp())*(1e4);
85 deltaPhi = newposition->DeltaPhi(*position);
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);
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);
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);
101 for (
int i = 0;
i < 3;
i ++){
102 hCartesianAveShift[
i]->Divide(hCartesianForward[
i],hStripesPerBin);
107 double x = low + ((high - low)/(1.0*nbins))*(
i+0.5);
110 double y = low + ((high - low)/(1.0*nbins))*(
j+0.5);
112 int xbin = hCartesianAveShift[0]->FindBin(x,y);
113 int ybin = hCartesianAveShift[1]->FindBin(x,y);
114 double xaveshift = (hCartesianAveShift[0]->GetBinContent(xbin))*(1
e-4);
115 double yaveshift = (hCartesianAveShift[1]->GetBinContent(ybin))*(1
e-4);
117 TVector3 shifted, original;
120 shifted.SetX(x+xaveshift);
121 shifted.SetY(y+yaveshift);
125 double raveshift = (shifted.Perp() - original.Perp())*(1e4);
126 double phiaveshift = shifted.DeltaPhi(original);
129 hCylindricalAveShift[0]->Fill(x,y,raveshift);
130 hCylindricalAveShift[1]->Fill(x,y,phiaveshift);
140 double minphi = -0.078539819;
141 double minr = 18.884615;
143 double minz = -1.3187500;
145 double maxphi = 6.3617253;
146 double maxr = 79.115387;
147 double maxz = 106.81875;
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);
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);
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);
158 double xshift, yshift, zshift, rshiftcart, phishiftcart;
161 double phi = minphi + ((maxphi - minphi)/(1.0*nphi))*(
i+0.5);
163 for(
int j = 0;
j < nr;
j++){
164 double r = minr + ((maxr - minr)/(1.0*nr))*(
j+0.5);
166 double x = r*cos(phi);
167 double y = r*sin(phi);
169 for(
int k = 0;
k <
nz;
k++){
170 double z = minz + ((maxz - minz)/(1.0*nz))*(
k+0.5);
172 xshift=(hCartesianAveShift[0]->Interpolate(x,y))*(1
e-4);
173 yshift=(hCartesianAveShift[1]->Interpolate(x,y))*(1
e-4);
174 zshift=(hCartesianAveShift[2]->Interpolate(x,y))*(1
e-4);
176 rshiftcart=(hCylindricalAveShift[0]->Interpolate(x,y))*(1
e-4);
177 phishiftcart=hCylindricalAveShift[1]->Interpolate(x,y);
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));
183 hCylindricalCMModel[0]->Fill(phi,r,z,rshiftcart*(1-z/105.5));
184 hCylindricalCMModel[1]->Fill(phi,r,z,phishiftcart*(1-z/105.5));
191 plots=TFile::Open(Form(
"CMModelsCart_Event%d.root",ifile),
"RECREATE");
192 hStripesPerBin->Write();
194 for(
int i = 0;
i < 3;
i++){
195 hCartesianForward[
i]->Write();
196 hCartesianAveShift[
i]->Write();
197 hCartesianCMModel[
i]->Write();
200 for(
int i = 0;
i < 2;
i++){
201 hCylindricalAveShift[
i]->Write();
202 hCylindricalCMModel[
i]->Write();
210 unsigned long after = now;
212 hTimePerEvent->Fill(after-before);
243 hTimePerEvent->Draw();
244 canvas->Print(
"CMDistortionRecoCart2.pdf",
"pdf");