StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
t0fitA.C
1 #include <stdlib.h>
2 #include <stdio.h>
3 #include <iostream>
4 
5 #include "TF1.h"
6 #include "TH1F.h"
7 #include "TFile.h"
8 #include "TMath.h"
9 #include "TCanvas.h"
10 #include "TTree.h"
11 #include "TSystem.h"
12 #include "tfit_CuCu.C"
13 
14 // par[]
15 // 0 - amplitude of step function
16 // 1 - edge of the step
17 // 2 - 1/sigma of step function
18 // 3 - amplitude of Gaus
19 // 4 - center of Gaus
20 // 5 - 1/sigma of Gaus
21 // 6 - constant
22 
23 Double_t tmin ( Double_t *v, Double_t *par )
24 {
25  Double_t x = v[0];
26  Double_t argS = (x - par[1])*par[2];
27  Double_t argG = (x - par[1])*par[2];
28  Double_t fitval = par[0]*TMath::Freq( argS )+
29  par[3]*TMath::Exp(-0.5*argG*argG); // +
30 // par[4];
31  return fitval;
32 }
33 const int Npar = 4;
34 
35 // par[]
36 // 0 - amplitude of step function
37 // 1 - edge of the step
38 // 2 - 1/sigma of step function
39 // 3 - constant when needed
40 
41 Double_t tmax ( Double_t *v, Double_t *par )
42 {
43  Double_t x = v[0];
44  Double_t argS = (par[1] - x )*par[2];
45  Double_t fitval = par[0]*TMath::Freq( argS );
46 // par[3];
47  return fitval;
48 }
49 const int NpaM = 3;
50 
51 class FitTree : public TObject {
52 public:
53  char * fName;
54  char * hName;
55  int barrel;
56  int ladder;
57  int wafer;
58  int hybrid;
59  double Entries;
60  int fitflag;
61  char * comment;
62  int fitID;
63  int Nparr;
64  double chi2;
65  double prr[Npar];
66  double err[Npar];
67 
68  int ftMflag;
69  char * commentM;
70  int ftMID;
71  int NparM;
72  double chi2M;
73  double prM[NpaM];
74  double erM[NpaM];
75 
76  FitTree(): Nparr(Npar),
77  NparM(NpaM){}
78 
79  ClassDef(FitTree,1)
80 };
81 
82  ClassImp(FitTree)
83 
84  struct PlotPar_t {
85  char * Name;
86  char * Title;
87  int nx;
88  int ny;
89  double xmin;
90  double xmax;
91  double ymin;
92  double ymax;
93  };
94  const PlotPar_t plotTB = // plots for time bins and anodes
95  { "timeB","time for 80 anodes", 256, 3, 0.,128., 0.,3. };
96 
97 
98 int t0fitA ( void) {
99 
100 // files *******************************************************
101 // char * fName = "afs_rhic.bnl.gov_star_users_fisyak_work_SvtSsdAlignment_Pass103_CuCuNoField_TpcOnlyPlotsTBNFP25rCut0.5cm.root";
102  char * fName = "afs_rhic.bnl.gov_star_users_fisyak_work_SvtSsdAlignment_Pass102_CuCu62PlotsTBNFP25rCut0.5cm.root";
103  TFile * f = new TFile( fName );
104 
105  char * outName = "tfit_CuCu62_L3"; // default file: fit2_par.root
106  char ftName[128] = "";
107  char foName[128] = "";
108  strcat(ftName, outName); strcat(ftName, ".root");
109  strcat(foName, outName); strcat(foName, ".txt");
110  cout << ftName << endl;
111  cout << foName << endl;
112  TFile fout(ftName,"RECREATE","fit parameters tree");
113  FILE * pFile = fopen (foName,"wt");
114  if (pFile == 0) return -1;
115 
116 // create a tree ********************************************************
117 
118  FitTree FitOut;
119  FitTree * pFitOut = &FitOut;;
120  TTree * tree = new TTree("Fit0T","t0 fit output tree");
121  tree->Branch("fit_res",&pFitOut,8000,99);
122 
123 // canvas ***************************************************************
124 
125  TCanvas *c1 = new TCanvas("fit", "t0fit", 300,10,700,650);
126 
127 //loops over barrels, ladders, wafers, hybrids
128 
129 const int NL[4] = {8, 12, 16, 20}; // ladders in the barrel
130 const int NW[4] = {4, 6, 7, 16}; // wafers on the ladder
131 
132 for (int barrel = 1; barrel <= 3; barrel++){
133 for (int ladder = 1; ladder <= NL[barrel-1]; ladder++){
134 for (int wafer = 1; wafer <= NW[barrel-1]; wafer++) {
135 for (int hybrid = 1; hybrid <= 2; hybrid++){
136 
137 // get histogram ***********************************************
138 
139  char nome[64];
140  sprintf( nome,"B%iL%02iW%iH%i",barrel,ladder,wafer,hybrid);
141 
142  TH1F *hs = new TH1F(Form("hybrid %i_%i_%i",ladder,wafer,hybrid),
143  nome,
144 // Form("B%i L%02i W%02i H%i",barrel,ladder,wafer,hybrid),
145  plotTB.nx, plotTB.xmin, plotTB.xmax );
146  hs->SetDirectory(0); // to disconnect this histogram from any files
147 
148  for (Int_t anode =1; anode <=3; anode++) {
149  TString Name("timeB");
150  Name += Form("L%02iB%iW%02iH%iA%i", ladder, barrel, wafer, hybrid, anode);
151 // TH1F *hist = (TH1F *) gDirectory->Get(Name);
152  TH1F *hist = (TH1F *) f->Get(Name);
153  if (! hist) continue;
154 // hist->SetLineColor(anode);
155  hs->Add(hist,1);
156  }
157  c1->Clear();
158  c1->cd(1);
159  hs->Draw("");
160 
161 // hs is filled up ------------------------------------------------------
162 
163  FitOut.fName = fName;
164  FitOut.hName = nome;
165  FitOut.barrel = barrel;
166  FitOut.ladder = ladder;
167  FitOut.wafer = wafer;
168  FitOut.hybrid = hybrid;
169  FitOut.Entries = hs->GetEntries(); // stats[0];
170  cout <<endl;
171  cout << " ********************************************************** " << endl;
172  cout << " Entries = " << FitOut.Entries << endl;
173  if ( FitOut.Entries < 100 ) cout << " histogram is empty: " << nome << endl;
174  if ( FitOut.Entries < 100 ) continue;
175  cout << " B " << FitOut.barrel
176  << " L " << FitOut.ladder
177  << " W " << FitOut.wafer
178  << " H " << FitOut.hybrid
179  << " nome: " << FitOut.hName << ""
180  << endl;
181 
182 // set initial values for parameters
183 
184  double par0 = FitOut.Entries/220.;
185  cout << endl;
186  int idx = barrel*10000 + ladder*100 + wafer*10 + hybrid;
187  double t0_init = 10.;
188  double tmax_init = 120.;
189  double d_length = 29928.; //in microns
190  int init = InitEdge ( idx, t0_init, tmax_init, d_length);
191  cout << " Init: " << init
192  << " t0 = " << t0_init
193  << " tmax = " << tmax_init
194  << " par0 = " << par0
195  << endl;
196  cout << endl;
197 
198 // t0 fit ****************************************
199 
200  TF1 * fit = new TF1("fit0",tmin, 0., 16., Npar);
201  fit->SetParameter(0,par0);
202  fit->SetParLimits(0,0.,2000.);
203  fit->SetParameter(1,t0_init);
204  fit->SetParameter(2,12.3);
205  fit->SetParLimits(2,0.,1000.);
206  fit->SetParameter(3,par0*2.);
207  fit->SetParLimits(3,0.,2000.);
208 // fit->SetParameter(4,1.0);
209 // fit->SetParameter(4,2.0);
210 // constant background
211 // fit->SetParameter(4,8.);
212 // fit->SetParLimits(6,0.,100.);
213 
214 
215  FitOut.fitflag = hs->Fit("fit0","LIBr");
216  cout << " fit status = " << FitOut.fitflag << endl;
217 
218  if ( FitOut.fitflag != 0 ) cout << "Fit failed: " << FitOut.fitflag << endl;
219 
220  FitOut.chi2 = fit->GetChisquare();
221  cout << " Chi**2 = " << FitOut.chi2 << endl;
222  cout << " Npar = " << FitOut.Nparr << endl;
223  for (int jpar = 0; jpar < Npar; jpar++){
224  FitOut.prr[jpar]=fit->GetParameter(jpar);
225  FitOut.err[jpar]=fit->GetParError(jpar);
226  cout << " par " << jpar
227  << " = " << FitOut.prr[jpar]
228  << " +- " << FitOut.err[jpar]
229  << endl;
230  if ( FitOut.prr[jpar] < 0.) cout << " WARNING! " << FitOut.prr[jpar] << endl;
231  }
232 
233  gPad->Update();
234  gPad->Draw();
235  c1->Update();
236  c1->Draw();
237 
238  char in;
239  cin >> in;
240  if ( in == 's') { cout << " operator intervention " << endl; return -3;}
241  if ( in == 'g') {FitOut.comment = "good "; FitOut.fitID = 1;}
242  else {FitOut.comment = "refit"; FitOut.fitID = -1;}
243  cout << " " << FitOut.hName
244  << " comment : " << FitOut.comment
245  << " fit ID : " << FitOut.fitID
246  << endl;
247  cout << endl;
248 
249 
250 // t_max fit ****************************************
251 
252  double Rmin = 116.;
253  if (init == 0) Rmin = tmax_init - 3.;
254  TF1 * ftM = new TF1("fitM",tmax, Rmin, 128., NpaM);
255  ftM->SetParameter(0, par0);
256  ftM->SetParLimits(0, 0.,2000.);
257  ftM->SetParameter(1, tmax_init);
258  ftM->SetParameter(2, 12.2);
259 // constant background
260 // ftM->SetParameter(3, 8.);
261 // ftM->SetParLimits(6,0.,100.);
262 
263 
264  FitOut.ftMflag = hs->Fit("fitM","LIBr");
265  cout << " fit status = " << FitOut.ftMflag << endl;
266 
267  if ( FitOut.ftMflag != 0 ) cout << "Fit failed: " << FitOut.ftMflag << endl;
268 
269  FitOut.chi2M = ftM->GetChisquare();
270  cout << " Chi**2 = " << FitOut.chi2M << endl;
271  cout << " Npar = " << FitOut.NparM << endl;
272  for (int jpar = 0; jpar < NpaM; jpar++){
273  FitOut.prM[jpar]=ftM->GetParameter(jpar);
274  FitOut.erM[jpar]=ftM->GetParError(jpar);
275  cout << " par " << jpar
276  << " = " << FitOut.prM[jpar]
277  << " +- " << FitOut.erM[jpar]
278  << endl;
279  if ( FitOut.prM[jpar] < 0.) cout << " WARNING! " << FitOut.prM[jpar] << endl;
280  }
281 
282  gPad->Update();
283  gPad->Draw();
284  c1->Update();
285  c1->Draw();
286 
287 // char in;
288  cin >> in;
289  if ( in == 'q') { cout << " quit " << endl; return -3;}
290  if ( in == 's') { cout << " stopped " << endl; fout.Write(); fout.Close(); fclose (pFile); return -4;}
291  if ( in == 'g') {FitOut.commentM = "good "; FitOut.ftMID = 1;}
292  else {FitOut.commentM = "refit"; FitOut.ftMID = -1;}
293  cout << " " << FitOut.hName
294  << " comment M: " << FitOut.commentM
295  << " fit ID M: " << FitOut.ftMID
296  << endl;
297 // cout << endl;
298 
299  tree->Fill();
300 
301 // fill the txt file ***********************************************************
302 
303 const double timeb = 40.; // in ns
304 const double d_veldef = 6.75; //"default" velosity constant, micron/ns
305  double d_velosity;
306 
307  int type = 0;
308 // int rdx = 0; // row index
309  int nrows = 216*2; // total no. of real rows in the table; For Db interface (where nrows = 50)
310  int npr =0;
311 // double tmin;
312 // double tmax;
313  double v[10] ={0};
314  char row_name[32];
315 
316  if ( FitOut.fitID>0 &&
317  FitOut.ftMID>0 &&
318  (FitOut.prM[1] - FitOut.prr[1])> 0. )
319  d_velosity = d_length/( FitOut.prM[1] - FitOut.prr[1])/timeb;
320  else d_velosity = d_veldef;
321  v[0] = d_velosity;
322 
323  sprintf(row_name,"B%iL%02iW%iH%i", barrel, ladder, wafer, hybrid);
324 
325  fprintf (pFile," { %i, %4i, %i, %i, %i, %2i, %i, %i, ",
326  type, idx, nrows, npr, barrel, ladder, wafer, hybrid);
327  fprintf (pFile,"%7.3f, %6.3f, %7.3f, %6.3f, %7.1f, ",
328  FitOut.prr[1], FitOut.err[1], FitOut.prM[1], FitOut.erM[1], d_length );
329  fprintf (pFile,"{%6.4f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f, %3.1f }}, ",
330  v[0],v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9] );
331  fprintf (pFile," // %2i %2i %2i %2i %s %s %s \n",
332  FitOut.fitflag, FitOut.fitID, FitOut.ftMflag, FitOut.ftMID,
333  FitOut.comment, FitOut.commentM, row_name );
334 
335 //******************************************************************************
336  delete hs;
337 
338 } //hybrid loop
339 } //wafer loop
340 } //ladder loop
341 } //barrel loop
342 
343  fout.Write();
344  fout.Close();
345 
346  fclose (pFile);
347 
348  return 0;
349 
350 } //end macro
PlotPar_t
Definition: t0fitA.C:84
FitTree
Definition: t0fitA.C:51
Generated by   doxygen 1.8.5