StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
t0fitM.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_CuCu200NF.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 t0fitM ( int hybrid,
99  int wafer,
100  int ladder,
101  int barrel,
102  int edge,
103  double Rmin,
104  double Rmax )
105  {
106  cout << endl;
107  cout << " H " << hybrid
108  << " W " << wafer
109  << " L " << ladder
110  << " B " << barrel
111  << " E " << edge
112  << endl;
113 
114  char o_name[64];
115  sprintf( o_name,"B%iL%02iW%iH%iE%i_leaf.root",barrel,ladder,wafer,hybrid,edge);
116  cout << " o_name: " << o_name << endl;
117 
118 // files *******************************************************
119 // char * fName = "afs_rhic.bnl.gov_star_users_fisyak_work_SvtSsdAlignment_Pass103_CuCuNoField_TpcOnlyPlotsTBNFP25rCut0.5cm.root";
120  char * fName = "afs_rhic.bnl.gov_star_users_fisyak_work_SvtSsdAlignment_Pass102_CuCu62PlotsTBNFP25rCut0.5cm.root";
121  TFile * f = new TFile( fName );
122  TFile fout(o_name,"RECREATE","fit parameters tree");
123 
124 // create a tree ********************************************************
125 
126  FitTree FitOut;
127  FitTree * pFitOut = &FitOut;;
128  TTree * tree = new TTree("Fit0T","t0 fit output tree");
129  tree->Branch("fit_res",&pFitOut,8000,99);
130 
131 // get histogram ***********************************************
132 
133  char nome[64];
134  sprintf( nome,"B%iL%02iW%iH%i",barrel,ladder,wafer,hybrid);
135 
136  TH1F *hs = new TH1F(Form("hybrid %i_%i_%i",ladder,wafer,hybrid),
137  nome,
138  plotTB.nx, plotTB.xmin, plotTB.xmax );
139  hs->SetDirectory(0); // to disconnect this histogram from any files
140 
141  for (Int_t anode =1; anode <=3; anode++) {
142  TString Name("timeB");
143  Name += Form("L%02iB%iW%02iH%iA%i", ladder, barrel, wafer, hybrid, anode);
144 // TH1F *hist = (TH1F *) gDirectory->Get(Name);
145  TH1F *hist = (TH1F *) f->Get(Name);
146  if (! hist) continue;
147 // hist->SetLineColor(anode);
148  hs->Add(hist,1);
149  }
150 
151 // canvas ***************************************************************
152 
153  TCanvas *c1 = new TCanvas("fit", "t0fit", 300,10,700,650);
154  c1->Clear();
155  c1->cd(1);
156  hs->Draw("");
157 
158 // hs is filled up ------------------------------------------------------
159 
160  FitOut.fName = fName;
161  FitOut.hName = nome;
162  FitOut.barrel = barrel;
163  FitOut.ladder = ladder;
164  FitOut.wafer = wafer;
165  FitOut.hybrid = hybrid;
166  FitOut.Entries = hs->GetEntries(); // stats[0];
167  cout << " Entries = " << FitOut.Entries << endl;
168  if ( FitOut.Entries < 100 ) cout << " histogram is empty: " << nome << endl;
169  cout << " B " << FitOut.barrel
170  << " L " << FitOut.ladder
171  << " W " << FitOut.wafer
172  << " H " << FitOut.hybrid
173  << " nome: " << FitOut.hName << ""
174  << endl;
175 
176 // t0 fit ****************************************
177 
178  double par0 = FitOut.Entries/220.;
179  cout << " par[0] set to " << par0 << endl;
180  cout << endl;
181  int idx = barrel*10000 + ladder*100 + wafer*10 + hybrid;
182  double t0_init = 10.;
183  double tmax_init = 117.;
184  int init = InitEdge ( idx, t0_init, tmax_init);
185  cout << " Init: " << init
186  << " t0 = " << t0_init
187  << " tmax = " << tmax_init
188  << endl;
189 
190  TF1 * fit = new TF1("fit0",tmin, Rmin, Rmax, Npar);
191  fit->SetParameter(0,par0);
192  fit->SetParLimits(0,0.,2000.);
193  fit->SetParameter(1,t0_init);
194  fit->SetParameter(2,3.3);
195  fit->SetParLimits(2,0.,1000.);
196  fit->SetParameter(3,par0*2.);
197  fit->SetParLimits(3,0.,2000.);
198  fit->SetParameter(4,10.0);
199  fit->SetParameter(5,3.3);
200 // constant background
201 // fit->SetParameter(4,8.);
202 // fit->SetParLimits(6,0.,100.);
203 
204 
205  FitOut.fitflag = 0;
206  if (edge ==0 ) FitOut.fitflag = hs->Fit("fit0","LIBr");
207  cout << " fit status = " << FitOut.fitflag << endl;
208 
209  if ( FitOut.fitflag != 0 ) cout << "Fit failed: " << FitOut.fitflag << endl;
210 
211  FitOut.chi2 = fit->GetChisquare();
212  cout << " Chi**2 = " << FitOut.chi2 << endl;
213  cout << " Npar = " << FitOut.Nparr << endl;
214  for (int jpar = 0; jpar < Npar; jpar++){
215  FitOut.prr[jpar]=fit->GetParameter(jpar);
216  FitOut.err[jpar]=fit->GetParError(jpar);
217  cout << " par " << jpar
218  << " = " << FitOut.prr[jpar]
219  << " +- " << FitOut.err[jpar]
220  << endl;
221  if ( FitOut.prr[jpar] < 0.) cout << " WARNING! " << FitOut.prr[jpar] << endl;
222  }
223 
224  gPad->Update();
225  gPad->Draw();
226  c1->Update();
227  c1->Draw();
228 
229  if ( edge == 0 && FitOut.fitflag == 0)
230  {FitOut.comment = "good"; FitOut.fitID = 2;}
231  else {FitOut.comment = "nofit"; FitOut.fitID = -2;}
232  cout << " " << FitOut.hName
233  << " comment: " << FitOut.comment
234  << "; fit ID: " << FitOut.fitID
235  << endl;
236  cout << endl;
237 
238 
239 // t_max fit ****************************************
240 
241  if (init == 0) Rmin = tmax_init - 3.;
242  TF1 * ftM = new TF1("fitM",tmax, Rmin, Rmax, NpaM);
243  ftM->SetParameter(0, par0);
244  ftM->SetParLimits(0, 0.,2000.);
245  ftM->SetParameter(1, tmax_init);
246  ftM->SetParameter(2, 2.2);
247 // constant background
248 // ftM->SetParameter(3, 8.);
249 // ftM->SetParLimits(3,0.,100.);
250 
251  FitOut.ftMflag = 0;
252  if (edge ==1 ) FitOut.ftMflag = hs->Fit("fitM","LIBr");
253  cout << " fit status = " << FitOut.ftMflag << endl;
254 
255  if ( FitOut.ftMflag != 0 ) cout << "Fit failed: " << FitOut.ftMflag << endl;
256 
257  FitOut.chi2M = ftM->GetChisquare();
258  cout << " Chi**2 = " << FitOut.chi2M << endl;
259  cout << " Npar = " << FitOut.NparM << endl;
260  for (int jpar = 0; jpar < NpaM; jpar++){
261  FitOut.prM[jpar]=ftM->GetParameter(jpar);
262  FitOut.erM[jpar]=ftM->GetParError(jpar);
263  cout << " par " << jpar
264  << " = " << FitOut.prM[jpar]
265  << " +- " << FitOut.erM[jpar]
266  << endl;
267  if ( FitOut.prM[jpar] < 0.) cout << " WARNING! " << FitOut.prM[jpar] << endl;
268  }
269 
270  gPad->Update();
271  gPad->Draw();
272  c1->Update();
273  c1->Draw();
274 
275  if ( edge == 1 && FitOut.ftMflag == 0 )
276  {FitOut.commentM = "good"; FitOut.ftMID = 2;}
277  else {FitOut.commentM = "nofit"; FitOut.ftMID = -2;}
278  cout << " " << FitOut.hName
279  << " comment M: " << FitOut.commentM
280  << "; fit ID M: " << FitOut.ftMID
281  << endl;
282  cout << endl;
283  cout << endl;
284 
285  tree->Fill();
286 
287 // delete hs;
288 
289  fout.Write();
290  fout.Close();
291 
292  return 0;
293 
294 } //end macro
295 
PlotPar_t
Definition: t0fitA.C:84
FitTree
Definition: t0fitA.C:51
Generated by   doxygen 1.8.5