Hall-D Software  alpha
FitTimeToDistance.C
Go to the documentation of this file.
1 double long_drift_func[3][3];
2 double short_drift_func[3][3];
3 double magnet_correction[2][2];
4
5 vector<double> cdc_drift_table;
6 double Bz_avg;
7
8 // Set values for the region cut
9 float deltaMin = -0.175, deltaMax = 0.175, tMin = 300, tMax = 1200;
10
11 unsigned int Locate(vector<double>&xx, double x){
12  int n=xx.size();
13  if (x==xx[0]) return 0;
14  else if (x==xx[n-1]) return n-2;
15
16  int jl=-1;
17  int ju=n;
18  int ascnd=(xx[n-1]>=xx[0]);
19  while(ju-jl>1){
20  int jm=(ju+jl)>>1;
21  if ( (x>=xx[jm])==ascnd)
22  jl=jm;
23  else
24  ju=jm;
25  }
26  return jl;
27 }
28
29 Double_t TimeToDistance( Double_t *x, Double_t *par){
30  Double_t d=0.0;
31  double delta = x[1]; // yAxis
32  double t = x[0]; // xAxis
33
34  // Cut out region in fit.
35  if (delta > deltaMax || delta < deltaMin) return 0.0;
36  if (delta < (((deltaMax - deltaMin) / (tMax - tMin))*(t - tMin) + deltaMin)) return 0.0;
37  // Variables to store values for time-to-distance functions for delta=0
38  // and delta!=0
39  double f_0=0.;
40  double f_delta=0.;
41  if (t > 0){
42  if (delta>0){
43  double a1=par[0];
44  double a2=par[1];
45  double a3=par[2];
46  double b1=par[3];
47  double b2=par[4];
48  double b3=par[5];
49  double c1=par[6];
50  double c2=par[7];
51  double c3=par[8];
52  // use "long side" functional form
53  double my_t=0.001*t;
54  double sqrt_t=sqrt(my_t);
55  double t3=my_t*my_t*my_t;
56  double delta_mag=fabs(delta);
57  double a=a1+a2*delta_mag;
58  double b=b1+b2*delta_mag;
59  double c=c1+c2*delta_mag+c3*delta*delta;
60  f_delta=a*sqrt_t+b*my_t+c*t3;
61  f_0=a1*sqrt_t+b1*my_t+c1*t3;
62
63  }
64  else{
65  double my_t=0.001*t;
66  double sqrt_t=sqrt(my_t);
67  double delta_mag=fabs(delta);
68
69  // use "short side" functional form
70  double a1=par[9];
71  double a2=par[10];
72  double a3=par[11];
73  double b1=par[12];
74  double b2=par[13];
75  double b3=par[14];
76  double c1=par[15];
77  double c2=par[16];
78  double c3=par[17];
79  double delta_sq=delta*delta;
80  double a=a1+a2*delta_mag+a3*delta_sq;
81  double b=b1+b2*delta_mag+b3*delta_sq;
82  f_delta=a*sqrt_t+b*my_t;
83  f_0=a1*sqrt_t+b1*my_t;
84
85  }
86
87  unsigned int max_index=cdc_drift_table.size()-1;
88  if (t>cdc_drift_table[max_index]){
89  d=f_delta;
90  return d;
91  }
92
93  // Drift time is within range of table -- interpolate...
94  unsigned int index=0;
95  index=Locate(cdc_drift_table,t);
96  double dt=cdc_drift_table[index+1]-cdc_drift_table[index];
97  double frac=(t-cdc_drift_table[index])/dt;
98  double d_0=0.01*(double(index)+frac);
99
100  double P=0.;
101  double tcut=250.0; // ns
102  if (t<tcut) {
103  P=(tcut-t)/tcut;
104  }
105  d=f_delta*(d_0/f_0*P+1.-P);
106  }
107  return d;
108 }
109
110
111 void FitTimeToDistance(TString inputROOTFile = "hd_root.root")
112 {
113  // Script for fitting the time to distance relation from data
114  TFile *thisFile = TFile::Open(inputROOTFile);
115  TH1I *Bz_hist = (TH1I *) thisFile->Get("/CDC_TimeToDistance/Bz");
116  TF2 *f1,*f2; const Int_t npar = 18;
117  bool isFieldOff = false;
118  if (Bz_hist == 0) isFieldOff = true;
119  f1 = new TF2("f1",TimeToDistance, 10, 1500, -0.3, 0.3, npar);
120  f2 = new TF2("f2",TimeToDistance, 10, 1500, -0.3, 0.3, npar);
121  //f1 = new TF2("f1",TimeToDistanceFieldOff, 0, 200, -0.18, 0.18, npar);
122  //f2 = new TF2("f2",TimeToDistanceFieldOff, 0, 200, -0.18, 0.18, npar);
123
124  TProfile *constants = (TProfile*) thisFile->Get("/CDC_TimeToDistance/CDC_TD_Constants");
125  long_drift_func[0][0] = constants->GetBinContent(101);
126  long_drift_func[0][1] = constants->GetBinContent(102);
127  long_drift_func[0][2] = constants->GetBinContent(103);
128  long_drift_func[1][0] = constants->GetBinContent(104);
129  long_drift_func[1][1] = constants->GetBinContent(105);
130  long_drift_func[1][2] = constants->GetBinContent(106);
131  long_drift_func[2][0] = constants->GetBinContent(107);
132  long_drift_func[2][1] = constants->GetBinContent(108);
133  long_drift_func[2][2] = constants->GetBinContent(109);
134  magnet_correction[0][0] = constants->GetBinContent(110);
135  magnet_correction[0][1] = constants->GetBinContent(111);
136  short_drift_func[0][0] = constants->GetBinContent(112);
137  short_drift_func[0][1] = constants->GetBinContent(113);
138  short_drift_func[0][2] = constants->GetBinContent(114);
139  short_drift_func[1][0] = constants->GetBinContent(115);
140  short_drift_func[1][1] = constants->GetBinContent(116);
141  short_drift_func[1][2] = constants->GetBinContent(117);
142  short_drift_func[2][0] = constants->GetBinContent(118);
143  short_drift_func[2][1] = constants->GetBinContent(119);
144  short_drift_func[2][2] = constants->GetBinContent(120);
145  magnet_correction[1][0] = constants->GetBinContent(121);
146  magnet_correction[1][1] = constants->GetBinContent(122);
147
148  for (unsigned int i=1; i<=78; i++){
149  cdc_drift_table.push_back(constants->GetBinContent(i));
150  }
151
152  int run = (int) constants->GetBinContent(125);
153
154  // So now you have the input to your function
155
156  Double_t parameters[npar] =
157  {long_drift_func[0][0], long_drift_func[0][1], long_drift_func[0][2],
158  long_drift_func[1][0], long_drift_func[1][1], long_drift_func[1][2],
159  long_drift_func[2][0], long_drift_func[2][1], long_drift_func[2][2],
160  short_drift_func[0][0], short_drift_func[0][1], short_drift_func[0][2],
161  short_drift_func[1][0], short_drift_func[1][1], short_drift_func[1][2],
162  short_drift_func[2][0], short_drift_func[2][1], short_drift_func[2][2]};
163
164
165  f1->SetParameters(parameters);
166  f2->SetParameters(parameters);
167
168  /*
169  double fractionalRange=0.15;
170  for(unsigned int i=0; i < npar; i++){
171  f1->SetParLimits(i,parameters[i]-fractionalRange*parameters[i],parameters[i]+fractionalRange*parameters[i]);
172  f2->SetParLimits(i,parameters[i]-fractionalRange*parameters[i],parameters[i]+fractionalRange*parameters[i]);
173  }
174  */
175
176  TProfile2D *profile = (TProfile2D *) thisFile->Get("/CDC_TimeToDistance/Predicted Drift Distance Vs Delta Vs t_drift");
177
178  //profile->Fit("f2");
179  gStyle->SetOptStat(0);
180
181  TCanvas *c1 = new TCanvas ("c1", "c1", 800, 600);
182  Double_t contours[21] =
183  { 0.00, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45,
184  0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0};
185
186  profile->SetContour(21, contours);
187  profile->Draw("colz");
188  f1->SetContour(21, contours);
189  profile->Draw("cont2 list same");
190  f1->Draw("cont2 list same");
191  c1->Update();
192  c1->SaveAs(Form("Before/Before_Run%i.png", run));
193
194  TProfile2D *profileRebin = profile->Rebin2D(4,4, "Rebin");
195  TFitResultPtr fr = profileRebin->Fit("f2", "S");
196  TCanvas *c2 = new TCanvas ("c2", "c2", 800, 600);
197  c2->cd();
198  profile->Draw("colz");
199  f2->SetContour(21, contours);
200  profile->Draw("cont2 list same");
201  f2->Draw("cont2 list same");
202  c2->Update();
203  c2->SaveAs(Form("After/After_Run%i.png", run));
204
205  TCanvas *c3 = new TCanvas ("c3", "c3", 800, 600);
206  f1->Draw("cont2 list");
207  f1->SetLineColor(3);
208  f2->Draw("cont2 list same");
209  c3->Update();
210  c3->SaveAs(Form("Combined/Combined_Run%i.png", run));
211
212  if ((Int_t) fr == 0){ // Fit converged with no errors
213  ofstream outputTextFile;
214  outputTextFile.open(Form("ccdb/ccdb_Format_%i.txt",run));
215  outputTextFile << fr->Parameter(0) << " " << fr->Parameter(1) << " " << fr->Parameter(2) << " " ;
216  outputTextFile << fr->Parameter(3) << " " << fr->Parameter(4) << " " << fr->Parameter(5) << " " ;
217  outputTextFile << fr->Parameter(6) << " " << fr->Parameter(7) << " " << fr->Parameter(8) << " " ;
218  outputTextFile << magnet_correction[0][0] << " " << magnet_correction[0][1] << endl;
219  outputTextFile << fr->Parameter(9) << " " << fr->Parameter(10) << " " << fr->Parameter(11) << " " ;
220  outputTextFile << fr->Parameter(12) << " " << fr->Parameter(13) << " " << fr->Parameter(14) << " " ;
221  outputTextFile << fr->Parameter(15) << " " << fr->Parameter(16) << " " << fr->Parameter(17) << " " ;
222  outputTextFile << magnet_correction[1][0] << " " << magnet_correction[1][1] << endl;
223  outputTextFile.close();
224  }
225
226  // Get Residual vs Drift Time
227  TH2I *resVsT = (TH2I*)thisFile->Get("/CDC_TimeToDistance/Residual Vs. Drift Time");
228  TCanvas *c4 = new TCanvas("c4", "c4", 800, 600);
229  resVsT->Draw("colz");
230  c4->SaveAs(Form("ResVsT/ResVsT_Run%i.png", run));
231
232  TCanvas *c5 = new TCanvas("c5", "c5", 800, 600);
233  gStyle->SetOptFit();
234  p = (TH1D*)resVsT->ProjectionY();
235  //p = (TH1D*)resVsT->ProjectionY("_py", resVsT->GetXaxis()->FindBin(100.), -1);
236  p->Draw();
237  p->Fit("gaus", "sqWR", "", -0.01, 0.01);
238  c5->SaveAs(Form("Proj/Proj_ResVsT_Run%i.png", run));
239
240  TCanvas *c6 = new TCanvas("c6", "c6", 1600, 900);
241  c6->Divide(3, 2, 0.001, 0.001);
242  c6->cd(1);
243  profile->SetContour(21, contours);
244  profile->Draw("colz");
245  f1->SetContour(21, contours);
246  profile->Draw("cont2 list same");
247  f1->Draw("cont2 list same");
248  c6->Update();
249
250  c6->cd(2);
251  profile->Draw("colz");
252  f2->SetContour(21, contours);
253  profile->Draw("cont2 list same");
254  f2->Draw("cont2 list same");
255  c6->Update();
256
257  c6->cd(3);
258  f1->Draw("cont2 list");
259  f1->SetLineColor(3);
260  f2->Draw("cont2 list same");
261  c6->Update();
262
263  c6->cd(4);
264  resVsT->Draw("colz");
265  c6->Update();
266
267  c6->cd(5);
268  p->Draw();
269  c6->Update();
270
271  c6->SaveAs(Form("Monitoring/Monitoring_Run%i.png", run));
272
273  gApplication->Terminate();
274  //return;
275 }
276
unsigned int Locate(vector< double > &xx, double x)
double short_drift_func[3][3]
double Bz_avg
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
#define c
float tMin
TFile * thisFile
static char index(char c)
Definition: base64.cpp:115
Double_t TimeToDistance(Double_t *x, Double_t *par)
Double_t c1[2][NMODULES]
Definition: tw_corr.C:68
float deltaMax
TFile f1(fnam)
Double_t c2[2][NMODULES]
Definition: tw_corr.C:69
double long_drift_func[3][3]
vector< double > cdc_drift_table
double magnet_correction[2][2]
float deltaMin
double sqrt(double)
void FitTimeToDistance(TString inputROOTFile="hd_root.root")
float tMax
TF1 * f2
Definition: FitGains.C:28
TCanvas * c4
TCanvas * c3