12 bool shiftPitch =
true;
14 bool shiftOffset =
false;
16 TFile *file = TFile::Open(inputFile);
17 TFile *
outFile =
new TFile(
"projection_result.root",
"RECREATE");
21 TH1D *resultHistsU[24];
22 TH1D *resultHistsV[24];
31 TH1D * constantsHist = (TH1D *) file->Get(
"FDC_InternalAlignment/CathodeAlignmentConstants");
34 ofstream constantsFile, pitchFile;
35 constantsFile.open(
"CathodeAlignment.txt");
36 pitchFile.open(
"StripPitchesV2.txt");
37 for(
int i=1; i<=24; i++){
38 double i_deltaPhiU = constantsHist->GetBinContent((i-1)*4+1);
39 double i_deltaU = constantsHist->GetBinContent((i-1)*4+2);
40 double i_deltaPhiV = constantsHist->GetBinContent((i-1)*4+3);
41 double i_deltaV = constantsHist->GetBinContent((i-1)*4+4);
43 double deltaX=constantsHist->GetBinContent((i-1)*2+401);
45 double phi_u = 75.0*TMath::DegToRad() + i_deltaPhiU;
46 double phi_v = (TMath::Pi() - 75.0*TMath::DegToRad()) + i_deltaPhiV;
48 double sinPhiU =
sin(phi_u);
49 double cosPhiU = cos(phi_u);
50 double sinPhiV =
sin(phi_v);
51 double cosPhiV = cos(phi_v);
52 double sinPhiUmPhiV =
sin(phi_u-phi_v);
53 double sinPhiUmPhiV2 = sinPhiUmPhiV*sinPhiUmPhiV;
54 double cosPhiUmPhiV = cos(phi_u-phi_v);
56 double deltaU=0., deltaV=0.;
58 double slopeU, slopeV;
59 double i_SP1_U = constantsHist->GetBinContent((i-1)*10+101);
60 double i_G1_U = constantsHist->GetBinContent((i-1)*10+102);
61 double i_SP2_U = constantsHist->GetBinContent((i-1)*10+103);
62 double i_G2_U = constantsHist->GetBinContent((i-1)*10+104);
63 double i_SP3_U = constantsHist->GetBinContent((i-1)*10+105);
64 double i_SP1_V = constantsHist->GetBinContent((i-1)*10+106);
65 double i_G1_V = constantsHist->GetBinContent((i-1)*10+107);
66 double i_SP2_V = constantsHist->GetBinContent((i-1)*10+108);
67 double i_G2_V = constantsHist->GetBinContent((i-1)*10+109);
68 double i_SP3_V = constantsHist->GetBinContent((i-1)*10+110);
70 double deltaSP1U,deltaG1U=0.,deltaSP2U,deltaG2U=0.,deltaSP3U;
71 double deltaSP1V,deltaG1V=0.,deltaSP2V,deltaG2V=0.,deltaSP3V;
73 for (
unsigned int iPlane=1; iPlane<=2; iPlane++){
76 if(iPlane == 1 )
sprintf(name,
"FDC_InternalAlignment/Cathode Projections/Plane %.2i u_res vs u",i);
77 else sprintf(name,
"FDC_InternalAlignment/Cathode Projections/Plane %.2i v_res vs v",i);
79 TH2I *
h = (TH2I *)file->Get(name);
80 if(iPlane == 1) hists2DU[i-1]=
h;
83 TF1 *g =
new TF1(
"g",
"gaus", -0.5, 0.5);
84 if(iPlane == 1)
sprintf(name,
"Plane %.2i u residual", i);
85 else sprintf(name,
"Plane %.2i v residual", i);
86 TH1D *result =
new TH1D(name,
"meas-pred", h->GetNbinsX(), -50., 50.);
87 if(iPlane == 1) resultHistsU[i-1]=result;
88 else resultHistsV[i-1]=result;
89 for (
int iX=3; iX <= h->GetNbinsX(); iX++){
91 if(iX%5 ==0) name = Form(
"Proj bin (%i)",iX);
93 TH1D * yProj = h->ProjectionY(name,iX-2, iX);
94 if (yProj->GetEntries() < 5) {
98 double max = yProj->GetBinCenter(yProj->GetMaximumBin());
99 g->SetParameter(1, max);
100 TFitResultPtr r = yProj->Fit(
"g",
"Q",
"", max - 0.01, max + 0.01 );
102 double mean = g->GetParameter(1);
103 double err = g->GetParError(1);
104 result->SetBinContent(iX,mean);
105 result->SetBinError(iX,err);
110 TString f_name = Form(
"f%.2i",i);
111 TString fLeft_name = Form(
"fLeft%.2i",i);
112 TString fMiddle_name = Form(
"fMiddle%.2i",i);
113 TString fRight_name = Form(
"fRight%.2i",i);
120 TF1 *
f =
new TF1(f_name,
"[0]+[1]*x", -50.,50.);
121 TF1 *fLeft =
new TF1(fLeft_name,
"[0]+[1]*(x+23.75)", -46.,-25.);
122 TF1 *fMiddle =
new TF1(fMiddle_name,
"[0]+[1]*x", -21.,21.);
123 TF1 *fRight =
new TF1(fRight_name,
"[0]+[1]*(x-23.75)", 25.,46.);
125 result->Fit(f_name,
"Q+rob=0.80");
126 if (iPlane ==1) deltaU = f->GetParameter(0);
127 else deltaV = f->GetParameter(0);
128 double slope = f->GetParameter(1);
129 if (iPlane == 1) slopeU = f->GetParameter(1);
130 else slopeV = f->GetParameter(1);
132 result->Fit(fLeft_name,
"QMR+rob=0.8");
133 if (iPlane == 1) deltaSP1U=0.5*(fLeft->GetParameter(1)-slope);
134 else deltaSP1V=0.5*(fLeft->GetParameter(1)-slope);
136 result->Fit(fMiddle_name,
"QMR+rob=0.85");
137 if (iPlane == 1) deltaSP2U=0.5*(fMiddle->GetParameter(1)-slope);
138 else deltaSP2V=0.5*(fMiddle->GetParameter(1)-slope);
140 result->Fit(fRight_name,
"QMR+rob=0.8");
141 if (iPlane == 1) deltaSP3U=0.5*(fRight->GetParameter(1)-slope);
142 else deltaSP3V=0.5*(fRight->GetParameter(1)-slope);
146 fMiddleU[i-1]=fMiddle;
151 fMiddleV[i-1]=fMiddle;
178 double mean1 =fLeft->GetParameter(0);
179 double mean2L = fMiddle->GetParameter(0)-23.75*fMiddle->GetParameter(1);
180 double mean2R = fMiddle->GetParameter(0)+23.75*fMiddle->GetParameter(1);
181 double mean3 =fRight->GetParameter(0);
184 if (fabs(fLeft->GetParameter(1)-fMiddle->GetParameter(1)) < 5E-4) {
186 deltaG1U = mean1 - mean2L;
189 deltaG1V = mean1 - mean2L;
192 if (fabs(fMiddle->GetParameter(1)-fRight->GetParameter(1)) < 5E-4) {
194 deltaG2U = mean2R - mean3;
197 deltaG2V = mean2R - mean3;
205 double avgCorrection = 0.5*(slopeU+slopeV);
206 cout <<
" shift plane " << i <<
" slopeU " << slopeU <<
" slopeV " << slopeV <<
" avgCorrection " << avgCorrection << endl;
207 cout <<
"dSP1U " << -deltaSP1U <<
" dG1U " << deltaG1U <<
" dSP2U " << -deltaSP2U <<
" dG2U " << deltaG2U <<
" dSP3U " << -deltaSP3U << endl;
208 cout <<
"dSP1V " << -deltaSP1V <<
" dG1V " << deltaG1V <<
" dSP2V " << -deltaSP2V <<
" dG2V " << deltaG2V <<
" dSP3V " << -deltaSP3V << endl;
210 pitchFile << (shiftPitch ? (i_SP1_U - deltaSP1U - avgCorrection):i_SP1_U) <<
" ";
211 pitchFile << (shiftGap ? (i_G1_U + deltaG1U) :i_G1_U) <<
" ";
212 pitchFile << (shiftPitch ? (i_SP2_U - deltaSP2U - avgCorrection):i_SP2_U) <<
" ";
213 pitchFile << (shiftGap ? (i_G2_U + deltaG2U) :i_G2_U) <<
" ";
214 pitchFile << (shiftPitch ? (i_SP3_U - deltaSP3U - avgCorrection):i_SP3_U) <<
" ";
215 pitchFile << (shiftPitch ? (i_SP1_V - deltaSP1V - avgCorrection):i_SP1_V) <<
" ";
216 pitchFile << (shiftGap ? (i_G1_V + deltaG1V) :i_G1_V) <<
" ";
217 pitchFile << (shiftPitch ? (i_SP2_V - deltaSP2V - avgCorrection):i_SP2_V) <<
" ";
218 pitchFile << (shiftGap ? (i_G2_V + deltaG2V) :i_G2_V) <<
" ";
219 pitchFile << (shiftPitch ? (i_SP3_V - deltaSP3V - avgCorrection):i_SP3_V) << endl;
221 double avgMag = (fabs(deltaU) + fabs(deltaV))/2.;
222 double halfAvg = avgMag/2.;
224 if(shiftOffset)constantsFile << i_deltaPhiU <<
" " << i_deltaU - deltaU/fabs(deltaU)*halfAvg <<
" " << i_deltaPhiV <<
" " << i_deltaV - deltaV/fabs(deltaV)*halfAvg << endl;
225 else constantsFile << i_deltaPhiU <<
" " << i_deltaU <<
" " << i_deltaPhiV <<
" " << i_deltaV << endl;
229 TCanvas *cResultU =
new TCanvas(
"cResultU",
"cResultU",1200,900);
230 cResultU->Divide(6,4);
231 for (
unsigned int i=0; i<24; i++){
233 hists2DU[i]->Draw(
"colz");
234 resultHistsU[i]->Draw(
"same");
235 fLeftU[i]->Draw(
"same");
236 fMiddleU[i]->Draw(
"same");
237 fRightU[i]->Draw(
"same");
240 cResultU->SaveAs(
"ResultU.png");
241 cResultU->SaveAs(
"ResultU.pdf");
244 TCanvas *cResultV =
new TCanvas(
"cResultV",
"cResultV",1200,900);
245 cResultV->Divide(6,4);
246 for (
unsigned int i=0; i<24; i++){
248 hists2DV[i]->Draw(
"colz");
249 resultHistsV[i]->Draw(
"same");
250 fLeftV[i]->Draw(
"same");
251 fMiddleV[i]->Draw(
"same");
252 fRightV[i]->Draw(
"same");
255 cResultV->SaveAs(
"ResultV.png");
256 cResultV->SaveAs(
"ResultV.pdf");
261 constantsFile.close();
sprintf(text,"Post KinFit Cut")
void FitCathodeProjections(TString inputFile="All.root")