Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FitCathodeProjections.C
Go to the documentation of this file.
1 #include <TMath.h>
2 #include <TH2I.h>
3 #include <TF1.h>
4 #include <TCanvas.h>
5 #include <TFile.h>
6 #include <TString.h>
7 #include <fstream>
8 #include <TMath.h>
9 
10 void FitCathodeProjections(TString inputFile = "All.root"){
11 
12  bool shiftPitch = true;
13  bool shiftGap = true;
14  bool shiftOffset = false;
15 
16  TFile *file = TFile::Open(inputFile);
17  TFile *outFile = new TFile("projection_result.root","RECREATE");
18 
19  TH2I *hists2DU[24];
20  TH2I *hists2DV[24];
21  TH1D *resultHistsU[24];
22  TH1D *resultHistsV[24];
23  TF1 *fLeftU[24];
24  TF1 *fLeftV[24];
25  TF1 *fMiddleU[24];
26  TF1 *fMiddleV[24];
27  TF1 *fRightU[24];
28  TF1 *fRightV[24];
29 
30  // The old constants are stored in a histogram
31  TH1D * constantsHist = (TH1D *) file->Get("FDC_InternalAlignment/CathodeAlignmentConstants");
32 
33  // Fit Function
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);
42 
43  double deltaX=constantsHist->GetBinContent((i-1)*2+401);
44 
45  double phi_u = 75.0*TMath::DegToRad() + i_deltaPhiU;
46  double phi_v = (TMath::Pi() - 75.0*TMath::DegToRad()) + i_deltaPhiV;
47 
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);
55 
56  double deltaU=0., deltaV=0.;
57 
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);
69 
70  double deltaSP1U,deltaG1U=0.,deltaSP2U,deltaG2U=0.,deltaSP3U;
71  double deltaSP1V,deltaG1V=0.,deltaSP2V,deltaG2V=0.,deltaSP3V;
72 
73  for (unsigned int iPlane=1; iPlane<=2; iPlane++){
74 
75  char name[100];
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);
78 
79  TH2I *h = (TH2I *)file->Get(name);
80  if(iPlane == 1) hists2DU[i-1]=h;
81  else hists2DV[i-1]=h;
82 
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++){
90  TString name;
91  if(iX%5 ==0) name = Form("Proj bin (%i)",iX);
92  else name = "_z";
93  TH1D * yProj = h->ProjectionY(name,iX-2, iX);
94  if (yProj->GetEntries() < 5) {
95  //result->SetBinContent(iX, 0.0);
96  continue;
97  }
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 );
101  if ((Int_t) r == 0){
102  double mean = g->GetParameter(1);
103  double err = g->GetParError(1);
104  result->SetBinContent(iX,mean);
105  result->SetBinError(iX,err);
106  }
107  //else result->SetBinContent(iX, 0.0);
108  }
109 
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);
114  if (iPlane == 1) {
115  f_name+="u";
116  fLeft_name+="u";
117  fMiddle_name+="u";
118  fRight_name+="u";
119  }
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.);
124 
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);
131 
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);
135 
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);
139 
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);
143 
144  if (iPlane ==1) {
145  fLeftU[i-1]=fLeft;
146  fMiddleU[i-1]=fMiddle;
147  fRightU[i-1]=fRight;
148  }
149  else {
150  fLeftV[i-1]=fLeft;
151  fMiddleV[i-1]=fMiddle;
152  fRightV[i-1]=fRight;
153  }
154 
155  if (shiftGap){
156  /*
157  int nBins = h->GetNbinsX();
158  TH1D * yProj1 = h->ProjectionY("projFoil1",1, nBins/4);
159  TH1D * yProj2 = h->ProjectionY("projFoil2",nBins/4 , 3*nBins/4);
160  TH1D * yProj3 = h->ProjectionY("projFoil3",3*nBins/4, nBins);
161 
162  double mean1 =0., mean2=0., mean3=0.;
163  double max = yProj1->GetBinCenter(yProj1->GetMaximumBin());
164  g->SetParameter(1, max);
165  TFitResultPtr r = yProj1->Fit("g", "Q", "", max - 0.01, max + 0.01 );
166  if ((Int_t) r == 0) mean1 = g->GetParameter(1);
167 
168  max = yProj2->GetBinCenter(yProj2->GetMaximumBin());
169  g->SetParameter(1, max);
170  r = yProj2->Fit("g", "Q", "", max - 0.01, max + 0.01 );
171  if ((Int_t) r == 0) mean2 = g->GetParameter(1);
172 
173  max = yProj3->GetBinCenter(yProj3->GetMaximumBin());
174  g->SetParameter(1, max);
175  r = yProj3->Fit("g", "Q", "", max - 0.01, max + 0.01 );
176  if ((Int_t) r == 0) mean3 = g->GetParameter(1);
177  */
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);
182 
183  // Only adjust the gap is the sections are flat.
184  if (fabs(fLeft->GetParameter(1)-fMiddle->GetParameter(1)) < 5E-4) {
185  if (iPlane == 1){
186  deltaG1U = mean1 - mean2L;
187  }
188  else{
189  deltaG1V = mean1 - mean2L;
190  }
191  }
192  if (fabs(fMiddle->GetParameter(1)-fRight->GetParameter(1)) < 5E-4) {
193  if (iPlane == 1){
194  deltaG2U = mean2R - mean3;
195  }
196  else{
197  deltaG2V = mean2R - mean3;
198  }
199  }
200  }
201  }
202 
203  // Calculate the new offsets from the fit parameters
204  // Need to account for avg pitch from the overall slopes
205  double avgCorrection = 0.5*(slopeU+slopeV); // Should be zero when all that is left are rotations of the plane
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;
209 
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;
220 
221  double avgMag = (fabs(deltaU) + fabs(deltaV))/2.;
222  double halfAvg = avgMag/2.;
223 
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;
226  }
227 
228  // Draw the results
229  TCanvas *cResultU = new TCanvas("cResultU","cResultU",1200,900);
230  cResultU->Divide(6,4);
231  for (unsigned int i=0; i<24; i++){
232  cResultU->cd(i+1);
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");
238  }
239 
240  cResultU->SaveAs("ResultU.png");
241  cResultU->SaveAs("ResultU.pdf");
242  cResultU->Write();
243 
244  TCanvas *cResultV = new TCanvas("cResultV","cResultV",1200,900);
245  cResultV->Divide(6,4);
246  for (unsigned int i=0; i<24; i++){
247  cResultV->cd(i+1);
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");
253  }
254 
255  cResultV->SaveAs("ResultV.png");
256  cResultV->SaveAs("ResultV.pdf");
257  cResultV->Write();
258 
259  outFile->Write();
260  outFile->Close();
261  constantsFile.close();
262 
263  return;
264 }
TFile * outFile
Definition: FitGains.C:15
sprintf(text,"Post KinFit Cut")
TF1 * f
Definition: FitGains.C:21
void FitCathodeProjections(TString inputFile="All.root")
double sin(double)