StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
TDraw.C
1 //#define TSPECTRA
2 #if !defined(__CINT__) || defined(__MAKECINT__)
3 #include "Riostream.h"
4 #include <stdio.h>
5 #include "TSystem.h"
6 #include "TMath.h"
7 #include "TH1.h"
8 #include "TH2.h"
9 #include "TH3.h"
10 #include "TStyle.h"
11 #include "TF1.h"
12 #include "TProfile.h"
13 #include "TTree.h"
14 #include "TChain.h"
15 #include "TFile.h"
16 #include "TNtuple.h"
17 #include "TCanvas.h"
18 #include "TClassTable.h"
19 #include "TFileSet.h"
20 #include "TDataSetIter.h"
21 #include "TDataSet.h"
22 #include "TDataSetIter.h"
23 #include "TClassTable.h"
24 //#include "DeDxTree.C"
25 #include "TMinuit.h"
26 #include "TSpectrum.h"
27 #include "StBichsel/Bichsel.h"
28 #include "TROOT.h"
29 #include "TString.h"
30 #include "TObjString.h"
31 #include "TLine.h"
32 #include "TText.h"
33 #include "TROOT.h"
34 #include "TList.h"
35 #include "TPolyMarker.h"
36 #include "TLegend.h"
37 #include "TArrayI.h"
38 #include "TArrayF.h"
39 #include "TArrayD.h"
40 #include "TRVector.h"
41 #include "TRMatrix.h"
42 //#include "TQtCanvas2Html.h"
43 #include "TInterpreter.h"
44 #else
45 class TMinuit;
46 class TF1;
47 class TH1F;
48 class TH2F;
49 class TH3F;
50 class TProfile;
51 class TH2D;
52 class TCanvas;
53 class TSpectrum;
54 class TSystem;
55 class Bichsel;
56 #endif
57 Int_t npeaks = 1;
58 static const Double_t nSigMax = 10;
59 TCanvas *c1, *c2;
60 TPad *selold = 0;
61 static const Double_t WaferLength = 2.9928;
62 static const Char_t *plotName[] = {
63  "duuH", "duvP", "dutuP", "duOvertuPuP", "duOvertuPvP",
64  "dvuP", "duvH", "dvtvP", "dvOvertvPuP", "dvOvertvPvP", // 10
65  // "duuP", "duvP", "dutuP", "duOvertuPuP", "duOvertuPvP",
66  // "dvuP",/* "duvH",*/ "dvtvP", "dvOvertvPuP", "dvOvertvPvP", // 9
67  "dXvsZ","dYvsZ","dZvsZ",
68  // "dXvsX","dXvsY","dYvsX",
69  // "dYvsY","dZvsX","dZvsY",
70  // 11
71  "dX4dx","dX4dy",// "dX4dz",
72  "dX4da","dX4db","dX4dg","dY4dx","dY4dy",
73  // "dY4dz",
74  "dY4da",
75  "dY4db","dY4dg",//"dZ4dx","dZ4dy",
76  // "dZ4dz",
77  "dZ4da","dZ4db","dZ4dg"};
78 const Int_t noHist = sizeof(plotName)/sizeof(Char_t *);//34;// 26;
79 const Int_t firstHL = 0;
80 const Int_t firstHG = 10;
81 //const Int_t firstHG = 9;
82 const Int_t firstHP = firstHG + 3;
83 Int_t firstH = firstHG;
84 Int_t lastH = noHist - 1;
85 
86 struct Rot_t {
87  Double_t dX, dY, dZ;
88  Double_t alpha, beta, gamma;
89  Char_t *Comment;
90 };
91 static const Int_t NPOL = 8;
92 //________________________________________________________________________________
93 Double_t PolN(Double_t *x,Double_t *par) {
94  Double_t sum = par[NPOL-1];
95  for (Int_t i = NPOL-2; i >= 0; i--) sum = x[0]*sum + par[i];
96  return sum;
97 }
98 //________________________________________________________________________________
99 Double_t Fit(TH1 *hist, Double_t xmin=-99, Double_t xmax=99, const Char_t *opt = "qer") {
100  Double_t prob = 0;
101  if (! hist) return prob;
102  TF1 *f = (TF1 *) gROOT->GetFunction("PolN");
103  if (! f) {
104  f = new TF1("PolN",PolN,xmin,xmax,NPOL+1);
105  }
106  for (Int_t i = 0; i <= NPOL; i++) {if (i <=0) f->SetParameter(i,0.0); else f->FixParameter(i,0);}
107  hist->Fit("PolN",opt,"",xmin,xmax);
108  f->ReleaseParameter(1);
109  hist->Fit("PolN",opt,"",xmin,xmax);
110  f = hist->GetFunction("PolN");
111  if (! f) return prob;
112  f->GetRange(xmin,xmax);
113  prob = f->GetProb();
114  Int_t NfitP = f->GetNumberFitPoints();
115  // if (prob > 1e-3 || NfitP < 10) return prob;
116  if (prob > 0 || NfitP < 10) return prob;
117  for (Int_t k = 2; k <= NPOL; k++) {
118  f->ReleaseParameter(k);
119  hist->Fit("PolN",opt,"",xmin,xmax);
120  f = hist->GetFunction("PolN");
121  prob = f->GetProb();
122  if (prob > 1e-3) break;
123  }
124  return prob;
125 }
126 //________________________________________________________________________________
127 void FitPolN(TH1 *hist, Double_t Xmin=-99, Double_t Xmax=99, const Char_t *opt = "qer") {
128  Double_t xmin = Xmax;
129  Double_t xmax = Xmin;
130  const Int_t nxbin = hist->GetNbinsX();
131  Int_t np = 0;
132  Double_t x = 0;
133  for (Int_t i = 1; i <= nxbin; i++) {
134  Double_t err = hist->GetBinError(i);
135  if (err <= 1e-6) {hist->SetBinError(i,0); continue;}
136  np++;
137  x = hist->GetBinLowEdge(i);
138  if (xmin > x) xmin = x;
139  x = hist->GetBinLowEdge(i+1);
140  if (xmax < x) xmax = x;
141  }
142  if (np < 5) return;
143  Double_t prob = Fit(hist,xmin,xmax,opt);
144  TF1 *f = hist->GetFunction("PolN");
145  Int_t NfitP = f->GetNumberFitPoints();
146  if (NfitP < 5) return;
147  Int_t NMaxCut = NfitP/2;
148  Int_t NfitPCut = NMaxCut;
149  if (NfitPCut < 15) NfitPCut = 15;
150  Int_t NCut = 0;
151  Double_t dev[1000];
152  Int_t indx[1000];
153  Int_t iter = 0;
154  while (prob < 1e-3 || iter < 5) {
155  f = hist->GetFunction("PolN");
156  prob = f->GetProb();
157  iter++;
158  Double_t chisq = f->GetChisquare();
159  NfitP = f->GetNumberFitPoints();
160  if (prob >= 1e-3 || chisq < 0.0 || NfitP < NfitPCut || NCut > NMaxCut) break;
161  Double_t devCut = 3*chisq/NfitP;
162  // cout << hist->GetName() << "\t min/max " << xmin << "/" << xmax;
163  Int_t np = 0;
164  for (Int_t i = 1; i <= nxbin; i++) {
165  dev[i] = 0;
166  Double_t err = hist->GetBinError(i);
167  if (err <= 0.0) continue;
168  Double_t val = hist->GetBinContent(i);
169  dev[i] = (val - f->Eval(hist->GetBinCenter(i)))/err;
170  dev[i] *= dev[i];
171  np++;
172  }
173  TMath::Sort(nxbin, dev, indx, kTRUE); // decreasing order
174  Int_t npp = (Int_t) (0.1*np) + 1;
175  Int_t nskipped = 0;
176  for (Int_t j = 0; j < npp; j++) {
177  Int_t i = indx[j];
178  if (dev[i] < devCut && nskipped > 0) break;
179  hist->SetBinError(i,0);
180  nskipped++;
181  NCut++;
182  }
183  prob = Fit(hist,xmin,xmax,opt);
184  if (prob > 1e-3) break;
185  }
186 }
187 //________________________________________________________________________________
188 Double_t GetRMS(TH1 *h, Double_t x1, Double_t x2) {
189  if (! h) return 0;
190  TAxis *fXaxis = h->GetXaxis();
191  Int_t i1 = h->FindBin(x1); if (i1 < fXaxis->GetFirst()) i1 = fXaxis->GetFirst();
192  Int_t i2 = h->FindBin(x2); if (i2 > fXaxis->GetLast()) i2 = fXaxis->GetLast();
193  Double_t w = 0, x = 0, rms = 0;
194  Double_t sumW = 0; Double_t sumX = 0; Double_t sumX2 = 0;
195  for (Int_t binx = i1; binx <= i2; binx++) {
196  x = fXaxis->GetBinCenter(binx);
197  w = TMath::Abs(h->GetBinContent(binx));
198  sumW += w;
199  sumX += w*x;
200  sumX2 += w*x*x;
201  }
202  if (sumW == 0.0) return 0;
203  x = sumX/sumW;
204  rms = TMath::Sqrt(sumX2/sumW - x*x);
205  return rms;
206 }
207 //________________________________________________________________________________
208 Double_t g2g(Double_t *xx, Double_t *par) {
209  Double_t x = xx[0];
210  Double_t A = TMath::Exp(par[0]);
211  Double_t mu1 = par[1];
212  Double_t sig1 = par[2];
213  if (A < 1 && TMath::Abs(x -mu1) < 3*sig1) {TF1::RejectPoint(); return 0;}
214  Double_t B = TMath::Exp(par[3]);
215  Double_t mu2 = par[4];
216  Double_t sig2 = par[5];
217  Double_t gra = par[6];
218  Double_t dev1 = (x - mu1)/sig1;
219  Double_t dev2 = (x - mu2)/sig2;
220  Double_t value = A*TMath::Exp(-0.5*dev1*dev1) + B*TMath::Exp(-0.5*dev2*dev2) + gra;
221  return value;
222 }
223 //________________________________________________________________________________
224 TF1 *InitGP() {
225  struct Par_t {
226  Char_t *Name;
227  Double_t p, pmin, pmax;
228  };
229 #if 0
230  TF1 *gp = new TF1("gp","exp([0])*(exp(-0.5*((x-[1])/[2])**2)/(2.50662827463100024*[2]) + [3])");
231  const Par_t par[4] = {
232  {"LNorm", 5., 0., 25},
233  {"mu", 0., -1. , 1. },
234  {"sigma", 0.2, 0.01, 1. },
235  {"grass", 0.0, 0.00, 1. }
236  };
237  for (Int_t i = 0; i < 4; i++) {
238  gp->SetParName(i,par[i].Name);
239  gp->SetParameter(i,par[i].p);
240  gp->SetParLimits(i,par[i].pmin, par[i].pmax);
241  }
242 #else
243  TF1 *gp = new TF1("gp",g2g,-100,100,7);
244  const Par_t par[7] = {
245  {"logN1", 5., 0., 25.},
246  {"mu1", 0., -1., 1.},
247  {"sigma1",0.01, 0.001, 0.10},
248  {"logN2", 1., 0., 25.},
249  {"mu2", 0., -1., 1.},
250  {"sigma2", 0.10, 0.01, 1.},
251  {"grass", 0.0, 0.00, 1.}
252  };
253  for (Int_t i = 0; i < 7; i++) {
254  gp->SetParName(i,par[i].Name);
255  gp->SetParameter(i,par[i].p);
256  gp->SetParLimits(i,par[i].pmin, par[i].pmax);
257  }
258 #endif
259  return gp;
260 }
261 //________________________________________________________________________________
262 Double_t fbackgr(Double_t *x, Double_t *par) {
263  // return TMath::Exp(par[0])*TMath::Gaus(x[0],par[1],par[2])+par[3];
264  Double_t bw = TMath::BreitWigner(x[0],par[1],par[2]);
265  Double_t bwl = -99;
266  if (bw > 0) bwl = TMath::Log(bw);
267  return TMath::Exp(par[0] + (1. + par[4])*bwl) + par[3];
268 }
269 //________________________________________________________________________________
270 Double_t fpeaks(Double_t *x, Double_t *par) {
271  Double_t result = 0;
272  for (Int_t p=0;p<npeaks;p++) {
273  Double_t norm = TMath::Exp(par[3*p]);
274  Double_t mean = par[3*p+1];
275  Double_t sigma = par[3*p+2];
276  // result += norm*TMath::Gaus(x[0],mean,sigma);
277  result += norm*TMath::BreitWigner(x[0],mean,sigma);
278  }
279  return result + fbackgr(x,&par[3*npeaks]);
280 }
281 //________________________________________________________________________________
282 TF1 *Peaks(TH1 *h=0) {
283  if (! h) return 0;
284  Double_t allcha, sumx, sumx2, x, val, rms, mean, xmin, xmax;
285  Int_t bin;
286  const Double_t sqrtpi = 2.506628;
287  TAxis *xax = h->GetXaxis();
288  Int_t hxfirst = xax->GetFirst();
289  Int_t hxlast = xax->GetLast();
290  xmin = xax->GetXmin();
291  xmax = xax->GetXmax();
292  Double_t valmax = h->GetBinContent(hxfirst);
293  Double_t binwidx = h->GetBinWidth(hxfirst);
294  allcha = sumx = sumx2 = 0;
295  Int_t np = 0;
296  for (bin=hxfirst;bin<=hxlast;bin++) {
297  x = h->GetBinCenter(bin);
298  val = TMath::Abs(h->GetBinContent(bin));
299  if (val <= 0) continue;
300  np++;
301  if (val > valmax) valmax = val;
302  sumx += val*x;
303  sumx2 += val*x*x;
304  allcha += val;
305  }
306  if (allcha == 0 || np <= 5) return 0;
307  mean = sumx/allcha;
308  rms = sumx2/allcha - mean*mean;
309  if (rms > 0) rms = TMath::Sqrt(rms);
310  else rms = 0;
311  if (rms == 0) rms = binwidx*(hxlast-hxfirst+1)/4;
312  Double_t constant = TMath::Log(0.5*(valmax+binwidx*allcha/(sqrtpi*rms)));
313  TF1 *fback = (TF1 *) gROOT->GetFunction("back");
314  if (! fback) fback = new TF1("back",fbackgr,xmin,xmax,5);
315  fback->SetRange(mean-rms,mean+rms);
316  // fback->SetRange(xmin,xmax);
317  // fback->SetNpx(1000);
318  fback->SetParameter(0,constant); fback->SetParLimits(0,constant-10,constant+30);
319  fback->SetParameter(1,mean); fback->SetParLimits(1,-1,1);
320  fback->SetParameter(2,rms/2.);
321  fback->SetParLimits(2,20e-4,rms);
322  fback->FixParameter(3,0);
323  fback->FixParameter(4,0);
324  h->Fit(fback,"q");
325  if (fback->GetProb() > 1e-3) return fback;
326  fback->SetParameter(1,mean);
327  fback->ReleaseParameter(3);
328  h->Fit(fback,"q");
329  if (fback->GetProb() > 1e-3) return fback;
330  fback->ReleaseParameter(4);
331  fback->SetParLimits(4,0.,10.);
332  h->Fit(fback,"q");
333 #ifndef TSPECTRA
334  return fback;
335 #else
336  //________________________________________________________________________________
337  npeaks = np;
338  TSpectrum s(2*npeaks);
339  Int_t nfound = s.Search(h,2,"");
340  // cout << "Found " << nfound << " candidate peaks to fit" << endl;
341  if (! nfound || nfound > 10) return fback;
342  npeaks = 0;
343  Float_t *xpeaks = s.GetPositionX();
344  Float_t *ypeaks = s.GetPositionY();
345  TArrayI idxT(nfound); Int_t *idx = idxT.GetArray();
346  TArrayF dT(nfound,ypeaks); Float_t *d = dT.GetArray();
347  TMath::Sort(nfound,d,idx,1);
348  TArrayD Par(3*nfound+4);
349  Double_t *par = Par.GetArray();
350  Int_t p;
351  for (p=0;p<nfound;p++) {
352  Double_t xp = xpeaks[idx[p]];
353  Int_t bin = xax->FindBin(xp);
354  Double_t yp = h->GetBinContent(bin);
355  if (yp-TMath::Sqrt(yp) < fback->Eval(xp)) continue;
356  par[3*npeaks] = TMath::Log(yp);
357  par[3*npeaks+1] = xp;
358  par[3*npeaks+2] = 3*binwidx;
359  npeaks++;
360  }
361  for (p = 0; p < 4; p++) {
362  par[3*npeaks+p] = fback->GetParameter(p);
363  }
364  TF1 *fit = new TF1("fit",fpeaks,xmin,xmax,3*npeaks+1);
365  for (int p = 0; p < npeaks; p++) {
366  fit->SetParName(3*p,Form("Norm%i",p)); fit->SetParLimits(3*p,par[3*p]-10,par[3*p]+30);
367  fit->SetParName(3*p+1,Form("Mu%i",p));
368  fit->SetParName(3*p+2,Form("Sigma%i",p));
369  }
370  fit->SetParName(3*npeaks,"NormB"); fit->SetParLimits(3*npeaks,par[3*npeaks]-10,par[3*npeaks]+30);
371  fit->SetParName(3*npeaks+1,"MuB");
372  fit->SetParName(3*npeaks+2,"SigmaB");
373  fit->SetParName(3*npeaks+3,"grass");
374 
375  fit->SetNpx(1000);
376  fit->SetParameters(par);
377  // TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
378  fit->SetParameters(par);
379  fit->SetLineColor(2);
380  h->Fit("fit","q");
381  return fit;
382 #endif
383 }
384 //______________________________________________________________________________
385 void SlicesYFit(TH2* h=0, Int_t binmin=0, Int_t binmax=0, Int_t cut=10, Option_t *option="QNRI") {
386  if (! h) return;
387  Int_t nbins = h->GetXaxis()->GetNbins();
388  if (binmin < 1) binmin = 1;
389  if (binmax > nbins) binmax = nbins;
390  if (binmax < binmin) {binmin = 1; binmax = nbins;}
391  TString opt = option;
392  opt.ToLower();
393  Int_t ngroup = 1;
394  if (opt.Contains("g2")) {ngroup = 2; opt.ReplaceAll("g2","");}
395  if (opt.Contains("g3")) {ngroup = 3; opt.ReplaceAll("g3","");}
396  if (opt.Contains("g4")) {ngroup = 4; opt.ReplaceAll("g4","");}
397  if (opt.Contains("g5")) {ngroup = 5; opt.ReplaceAll("g5","");}
398  Int_t npar = 3;
399  if (npar <= 0) return;
400  Double_t *parsave = new Double_t[npar];
401 
402  //Create one histogram for each function parameter
403  Int_t ipar;
404  TH1D **hlist = new TH1D*[npar];
405  char *name = new char[2000];
406  char *title = new char[2000];
407  const TArrayD *bins = h->GetXaxis()->GetXbins();
408  for (ipar=0;ipar<npar;ipar++) {
409  sprintf(name,"%s_%d",h->GetName(),ipar);
410  sprintf(title,"Fitted value of par[%d]",ipar);
411  delete gDirectory->FindObject(name);
412  if (bins->fN == 0) {
413  hlist[ipar] = new TH1D(name,title, nbins, h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
414  } else {
415  hlist[ipar] = new TH1D(name,title, nbins,bins->fArray);
416  }
417  hlist[ipar]->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
418  }
419  sprintf(name,"%s_chi2",h->GetName());
420  delete gDirectory->FindObject(name);
421  TH1D *hchi2 = new TH1D(name,"chisquare", nbins, h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax());
422  hchi2->GetXaxis()->SetTitle(h->GetXaxis()->GetTitle());
423 
424  //Loop on all bins in X, generate a projection along Y
425  Int_t bin;
426  Int_t nentries;
427  for (bin=binmin;bin<=binmax-ngroup+1;bin++) {
428  TH1D *hpy = h->ProjectionY(Form("%s_bin_%i",h->GetName(),bin),bin,bin+ngroup-1,"e");
429  if (hpy == 0) continue;
430  nentries = Int_t(hpy->GetEntries());
431  if (nentries == 0 || nentries < cut) {delete hpy; continue;}
432  TF1 *f1 = Peaks(hpy);
433  if (f1) {
434  Int_t npfits = f1->GetNumberFitPoints();
435  // if (npfits > npar && npfits >= cut && f1->GetProb() > 1.e-3 && f1->GetParameter(1) < 0.5) {
436  Double_t mu = f1->GetParameter(1);
437  if (TMath::Abs(mu) < 0.5 && npfits > npar && npfits >= cut && (f1->GetProb() > 1.e-3 || f1->GetChisquare()/(npfits-npar) < 5.)) {
438  Int_t biny = bin + ngroup/2;
439  for (ipar=0;ipar<npar;ipar++) {
440  // hlist[ipar]->Fill(h->GetXaxis()->GetBinCenter(biny),f1->GetParameter(ipar));
441  hlist[ipar]->SetBinContent(biny,f1->GetParameter(ipar));
442  hlist[ipar]->SetBinError(biny,f1->GetParError(ipar));
443  }
444  hchi2->Fill(h->GetXaxis()->GetBinCenter(biny),f1->GetChisquare()/(npfits-npar));
445  }
446 #if 0
447  hpy->Draw();
448  c1->Update();
449  cout << "type something" << endl;
450  Int_t i;
451  cin >> i;
452 #endif
453  delete hpy;
454  }
455  }
456  delete [] parsave;
457  delete [] name;
458  delete [] title;
459  delete [] hlist;
460 }
461 //________________________________________________________________________________
462 Double_t STcheb(Int_t N, Double_t *par, Double_t x) {// N polynome degree, dimension is par[N+1]
463  if (N < 0 || N > 12) return 0;
464  Double_t T0 = 1;
465  Double_t T1 = 2*x - 1;
466  Double_t T2;
467  Double_t Sum = par[0]*T0;
468  if (N >= 1) {
469  T1 = 2*x - 1;
470  Sum += par[1]*T1;
471  for (int n = 2; n <= N; n++) {
472  T2 = 2*(2*x - 1)*T1 - T0;
473  Sum += par[n]*T2;
474  T0 = T1;
475  T1 = T2;
476  }
477  }
478  return Sum;
479 }
480 //________________________________________________________________________________
481 Double_t STchebN(Double_t *x, Double_t *par) {
482  Int_t N = (Int_t) par[0];
483  return STcheb(N,par+1,-x[0]-0.1);
484 }
485 //________________________________________________________________________________
486 Double_t STchebP(Double_t *x, Double_t *par) {
487  Int_t N = (Int_t) par[0];
488  return STcheb(N,par+1,x[0]-0.1);
489 }
490 //________________________________________________________________________________
491 void FitG(TFile *f, Int_t i, Int_t j, Int_t nx, Int_t s, Int_t &s1, Int_t &s2, TF1 *gp, ofstream &out,
492  Double_t FitR[6], Double_t dFitR[6], Double_t LSFit[6], Double_t dLSFit[6],TString &name) {
493  TString line("");
494  TString comment("");
495  TH2F *h = 0;
496  name = plotName[j];
497  if (i == 6) {name += "AllSvt";}
498  if (i == 7) {name += "AllSsd";}
499  if (i == 8) {name += "AllSvtSsd";}
500  if (c2) c2->cd();
501  if (i < 6) h = (TH2F *) f->Get(Form("%s%i",plotName[j],i));
502  else {
503  TH2F *h2 = (TH2F *) f->Get(Form("%s%i",plotName[j],s1));
504  if (! h2) { cout << "Histogram " << Form("%s%i",plotName[s1]) << " is not found" << endl; return;}
505  h = new TH2F(*h2);
506  h->SetName(name);
507  TString Title(h->GetTitle());
508  Int_t index = Title.Index("Sector");
509  if (index < 0) index = Title.Index("Clam");
510  if (index > 0) {
511  TString t(Title,index);
512  t += "All";
513  h->SetTitle(t);
514  }
515  for (s = s1+1; s <= s2; s++) {
516  h2 = (TH2F *) f->Get(Form("%s%i",plotName[j],s));
517  h->Add(h2);
518  }
519  }
520  if (! h) {cout << "Histogram for i/j = " << i << "/" << j << endl; return;}
521  h->SetXTitle(f->GetName());
522  TProfile *prof = h->ProfileX();
523  prof->SetMarkerStyle(24);
524  prof->SetMarkerColor(6);
525  TH1 *sp = h->ProjectionY("_py",-1,-1,"e");
526  sp->Fit("gaus","q");
527  Double_t Mu = 0;
528  Double_t dMu = 0;
529  TString HistName(h->GetName());
530  TF1 *gaus = sp->GetFunction("gaus");
531  gp->SetParameters(TMath::Log(gaus->GetParameter(0)),gaus->GetParameter(1),TMath::Abs(gaus->GetParameter(2)),0.);
532  sp->Fit(gp,"q");
533  Mu = sp->GetFunction("gp")->GetParameter(1);
534  if (j >= firstHG && j < firstHP) dMu = sp->GetFunction("gp")->GetParError(1);
535  Double_t *params = gp->GetParameters();
536  params[0] -= TMath::Log(100.);
537  gp->SetParameters(params);
538  // SlicesYFit(h,0,0,10,"qnig3"); //g3
539  SlicesYFit(h,0,0,10,"qni"); //g3
540  TH1 *fit = (TH1 *) gDirectory->Get(Form("%s_1",h->GetName()));
541  // TH1 *sig = (TH1 *) gDirectory->Get(Form("%s_2",h->GetName()));
542  // TH1 *gra = (TH1 *) gDirectory->Get(Form("%s_3",h->GetName()));
543  Double_t slope = 0;
544  Double_t dslope = 0;
545  TLegend *leg = new TLegend(0.2,0.2,0.8,0.3,"");
546  leg->SetTextSize(0.033);
547  if (fit) {
548  fit->SetTitle(h->GetTitle());
549  fit->SetMarkerStyle(20);
550  fit->SetMarkerColor(1);
551  fit->SetMaximum(0.2);
552  fit->SetMinimum(-.2);
553  fit->SetStats(1);
554  Double_t zmax = 99;
555  FitPolN(fit,-zmax,zmax);
556  // fit->Fit("PolN","eqr","",-zmax,zmax);
557  TF1 *pol1 = fit->GetFunction("PolN");
558  Double_t prob = 0;
559  if (pol1) {
560  prob = pol1->GetProb();
561  // if (prob <= 1e-3) {
562  // fit = prof;
563  // FitPolN(fit,-zmax,zmax);
564  // pol1 = fit->GetFunction("PolN");
565  // prob = pol1->GetProb();
566  // }
567  // if (prob > 1.e-7) {
568  slope = pol1->GetParameter(1);
569  dslope = pol1->GetParError(1);
570  // }
571  }
572  static const Char_t *dXYZ[3] = {"dX", "dY", "dZ"};
573  static const Char_t *abc[6] = {"=> dx", "=> dy", "=> dz", "=> alpha","=> beta","=> gamma"};
574  TString Name(h->GetName());
575  TString Title(h->GetTitle());
576  line = "";
577  if (j >= firstHG && j < firstHP) {
578  if (dMu > 0) {
579  for (Int_t m = 0; m < 3; m++) {
580  if (Name.BeginsWith(dXYZ[m])) {
581  Double_t mu = -1e4*Mu;
582  Double_t dmu = 1e4*dMu;
583  line += Form("|%7.2f+-%5.2f",mu,dmu);
584  Double_t dev = mu - LSFit[m];
585  Double_t sdev = TMath::Sqrt(dmu*dmu+dLSFit[m]*dLSFit[m]);
586  if (dLSFit[m] == 0 || sdev > 0 && TMath::Abs(dev/sdev) < nSigMax) {
587  Double_t dMu2 = dMu*dMu;
588  FitR[m] += -Mu/dMu2;
589  dFitR[m] += 1./dMu2;
590  line += "A";
591  } else line +="R";
592  if (m == 2) comment = Form("| slope = %7.2f+-%5.2f",1e3*slope, 1e3*dslope);
593  }
594  else line += Form("| ");
595  // cout << line << endl;
596  }
597  for (Int_t m = 3; m < 6; m++) {
598  if (dslope > 0 && Title.Contains(abc[m])) {
599  Double_t mu = 1e3*slope;
600  Double_t dmu = 1e3*dslope;
601  line += Form("|%7.2f+-%5.2f",mu,dmu);
602  Double_t dev = mu - LSFit[m];
603  Double_t sdev = TMath::Sqrt(dmu*dmu+dLSFit[m]*dLSFit[m]);
604  if (dLSFit[m] == 0 || sdev > 0 && TMath::Abs(dev/sdev) < nSigMax) {
605  Double_t dslope2 = dslope*dslope;
606  FitR[m] += slope/dslope2;
607  dFitR[m] += 1./dslope2;
608  line += "A";
609  } else line +="R";
610  }
611  else line += Form("| ");
612  // cout << line << endl;
613  }
614  }
615  if (pol1)
616  leg->AddEntry(pol1,Form("Mu = %7.2f +- %5.2f (mkm) Slope = %7.2f +- %5.2f (mrad)",
617  1e4*Mu, 1e4*dMu, 1e3*slope, 1e3*dslope));
618  } else {
619  if (j >= firstHP) {
620  for (Int_t m = 0; m < 6; m++) {
621  if (dslope > 0 && Title.Contains(abc[m])) {
622  Double_t scale = 1e4;
623  if (m >= 3) scale = 1.e3;
624  Double_t mu = scale*slope;
625  Double_t dmu = scale*dslope;
626  line += Form("|%7.2f+-%5.2f",mu,dmu);
627  Double_t dev = mu - LSFit[m];
628  Double_t sdev = TMath::Sqrt(dmu*dmu+dLSFit[m]*dLSFit[m]);
629  if (dLSFit[m] == 0 || sdev > 0 && TMath::Abs(dev/sdev) < nSigMax) {
630  Double_t dslope2 = dslope*dslope;
631  FitR[m] += slope/dslope2;
632  dFitR[m] += 1./dslope2;
633  line += "A";
634  } else line +="R";
635  if (pol1) {
636  TString legT(abc[m]);
637  legT.ReplaceAll("=> d","#Delta ");
638  legT.ReplaceAll("=> ","#");
639  legT += Form(" = %8.2f +- %8.2f",mu,dmu);
640  if (scale < 5.e3) legT += "(mrad)";
641  else legT += "(mkm)";
642  legT += Form(" prob = %4.3f",prob);
643  leg->AddEntry(pol1,legT);
644  }
645  }
646  else line += Form("| ");
647  }
648  }
649  }
650  line += "|"; line += fit->GetName(); line += "/"; line += h->GetTitle();// line += "\t"; line += f->GetName();
651  line += comment;
652  cout << line << endl;
653  out << line << endl;
654 
655  Int_t ij = i + nx*(j-firstH) + 1;
656  c1->cd(ij)->SetLogz(1);
657  h->SetMinimum(1);
658 #if 0
659  if (h) h->DrawCopy("colz");
660  if (prof) prof->DrawCopy("same");
661 #else
662  if (h) h->Draw("colz");
663  if (prof) prof->Draw("same");
664 #endif
665  if (fit) {
666 #if 0
667  fit->DrawCopy("same");
668 #else
669  fit->Draw("same");
670 #endif
671  TF1 *pol1 = fit->GetFunction("PolN");
672  if (pol1) {pol1->SetLineColor(2); pol1->Draw("same");}
673  // TPaveStats *st = (TPaveStats*) fit->FindObject("stats");
674  // if (st) {
675  // st->SetX1NDC(0.1);
676  // st->SetX2NDC(0.5);
677  // st->Draw();
678  // }
679  }
680  if (leg->GetEntry()) leg->Draw();
681  // cout << f->GetName() << "\t" << h->GetName() << "/" << h->GetTitle()
682  // << "\tMu = " << Mu << " +/- " << dMu
683  // << "\tSlope = " << slope << " +/- " << dslope << endl;
684  // static const Char_t *blank10 = Form("| ");
685  }
686 }
687 //________________________________________________________________________________
688 void TDrawG(Int_t sectors=8) {
689  Int_t s, s1, s2;
690  s = s1 = s2 = sectors;
691  Int_t nx = 1;
692  if (s > 5) nx = s; // 8 => sum for SVT and SSD
693  Int_t ny = lastH - firstH + 1;
694  Int_t scaleX = 60; //600/nx;
695  Int_t scaleY = 40; //800/ny;
696  // Int_t scale = TMath::Min(scaleX,scaleY);
697  c2 = new TCanvas();
698  c1 = new TCanvas("Global","SVT ClamShell and SSD sectors",10,10,10+scaleX*nx,10+scaleY*ny);
699  cout << "nx/ny = " << nx << "/" << ny << endl;
700  c1->Divide(nx,ny);
701  ofstream out;
702  TString Out("Results.");
703  Out += gSystem->BaseName(gDirectory->GetName());
704  Out.ReplaceAll(".root","");
705  Out.ReplaceAll(" ","");
706  if (gSystem->AccessPathName(Out)) out.open(Out, ios::out); //"Results.list",ios::out | ios::app);
707  else out.open(Out, ios::app);
708  TF1 *gp = InitGP();
709  TCollection *files = gROOT->GetListOfFiles();
710  TIter next(files);
711  TFile *f = (TFile *) next();
712 // TFile *f = (TFile *) gDirectory;
713  if (! f) return;
714  out << "____________________________________________________________________________________________________" << endl;
715  out << "|dX mkm |dY mkm |dZ mkm |alpha mrad |beta mrad |gamma mrad |Comment" << endl;
716  cout << "____________________________________________________________________________________________________" << endl;
717  cout << "|dX mkm |dY mkm |dZ mkm |alpha mrad |beta mrad |gamma mrad |Comment" << endl;
718  for (Int_t i = 0; i < nx; i++) {
719  f->cd();
720  out << "______________________________________________________________________________________________ " << f->GetName() << endl;
721  cout << "______________________________________________________________________________________________ " << f->GetName() << endl;
722  Double_t FitR[6], dFitR[6], LSFit[6], dLSFit[6];
723  memset (FitR, 0, 6*sizeof(Double_t));
724  memset (dFitR, 0, 6*sizeof(Double_t));
725  memset (LSFit, 0, 6*sizeof(Double_t));
726  memset (dLSFit, 0, 6*sizeof(Double_t));
727  TString line("");
728  TString lTitle("");
729  TH1D *LSF = (TH1D *) f->Get("LSF");
730  TString name;
731  s = s1 = s2 = i;
732  if (i == 6) {s1 = 0; s2 = 1;}
733  if (i == 7) {s1 = 2; s2 = 5;}
734  if (i == 8) {s1 = 0; s2 = 5;}
735  if (LSF) {
736  Double_t *array = LSF->GetArray();
737  TRVector AmX(6);
738  TRSymMatrix S(6);
739  for (s = s1; s <= s2; s++) {
740  Int_t im = 1 + 28*s;
741  Int_t is = im + 6;
742  AmX += TRVector(6,array+im);
743  S += TRSymMatrix(6,array+is);// cout << "S " << S << endl;
744  }
745  TRSymMatrix SInv(S,TRArray::kInverted);// cout << "SInv " << SInv << endl;
746  TRVector X(SInv,TRArray::kSxA,AmX); //cout << "X " << X << endl;
747  TString line("");
748  for (Int_t l = 0; l < 6; l++) {
749  Double_t scale = 1e4;
750  if (l > 2) scale = 1e3;
751  LSFit[l] = scale*X(l);
752  dLSFit[l] = scale*TMath::Sqrt(SInv(l,l));
753  line += Form("|%7.2f+-%5.2f ",LSFit[l],dLSFit[l]);
754  if (SInv(l,l) > 0) {
755  FitR[l] += X(l)/SInv(l,l);
756  dFitR[l] += 1./SInv(l,l);
757  }
758  }
759  line += "|"; line += LSF->GetName(); line += "/";
760  if (i < 6) line += LSF->GetTitle();
761  else {
762  if (i == 6) line += "Sum Over Svt Shells";
763  if (i == 7) line += "Sum Over Ssd Sectors";
764  if (i == 8) line += "Sum Over Svt+Ssd";
765  }
766  cout << line << endl;
767  out << line << endl;
768  } // LSF
769  for (Int_t j = firstH; j <= lastH; j++) {//cout << "Plot:" << plotName[j] << endl;
770  FitG(f,i,j,nx,s,s1,s2,gp, out, FitR, dFitR, LSFit, dLSFit, name);
771  }
772  for (Int_t m = 0; m < 6; m++) {
773  if (dFitR[m] > 0) {
774  Double_t scale = 1e4;
775  if (m > 2) scale = 1e3;
776  FitR[m] = scale*FitR[m]/dFitR[m];
777  dFitR[m] = scale/TMath::Sqrt(dFitR[m]);
778  line += Form("|%7.2f+-%5.2f ", FitR[m],dFitR[m]);
779  } else {
780  line += Form("| ");
781  }
782  }
783  line += "| Average for ";
784  if (i == 8) line += "All Svt + Ssd";
785  else if (i == 7) line += "All Ssd";
786  else if (i == 6) line += "All Svt";
787  else if (i < 2) line += Form("SVT ClamShell %i",i);
788  else if (i < 6) line += Form("SSD Sector %i",i-1);
789  cout << line << endl;
790  out << line << endl;
791  }
792  out.close();
793 }
794 //________________________________________________________________________________
795 void TDrawL(Int_t iHist=-1, Int_t barrel = 4, Int_t ladder = 0, Int_t wafer = 0) {
796  static const Int_t NlPerBarrel[4] = {8, 12, 16, 20};
797  if (barrel < 1 || barrel > 4) {cout << "Wrong barrel no. " << barrel << endl; return;}
798  Int_t NL = NlPerBarrel[barrel-1];
799  Int_t l1 = 1;
800  Int_t l2 = NL;
801  if (ladder) {l1 = l2 = ladder;}
802  Int_t nx = l2 - l1 + 1;
803  firstH = firstHL; lastH = firstHG - 1;
804  if (iHist >= 0) {firstH = lastH = iHist;}
805  Int_t ny = lastH - firstH + 1;
806  Int_t scaleX = 60; //600/nx;
807  Int_t scaleY = 40; //800/ny;
808  // Int_t scale = TMath::Min(scaleX,scaleY);
809  c1 = new TCanvas(Form("Ladder_Barrel_%i",barrel),Form("Barrel %i, Ladder %i, Wafer %i",barrel,ladder,wafer) ,10,10,10+scaleX*nx,10+scaleY*ny);
810  cout << "nx/ny = " << nx << "/" << ny << endl;
811  c1->Divide(nx,ny);
812  ofstream out;
813  ofstream outC;
814  TString Out("Results.Barrel_");
815  Out += Form("%i_",barrel);
816  Out += gSystem->BaseName(gDirectory->GetName());
817  Out.ReplaceAll(".root","");
818  TDatime t;
819  // Out += t.AsString();
820  Out.ReplaceAll(" ","");
821  if (gSystem->AccessPathName(Out)) out.open(Out, ios::out); //"Results.list",ios::out | ios::app);
822  else out.open(Out, ios::app);
823  Out += ".h";
824  if (gSystem->AccessPathName(Out)) outC.open(Out, ios::out); //"Results.list",ios::out | ios::app);
825  else outC.open(Out, ios::app);
826  TF1 *gp = InitGP();
827  TH2 *h = 0;
828  Int_t head = 0;
829  for (Int_t i = 0; i < nx; i++) {
830  if (! head) {
831  out << "________________________________________________________________________________________________________" << endl;
832  out << "|du mkm |dv mkm |dw mkm |alpha mrad |beta mrad |gamma mrad |Comment" << endl;
833  cout << "________________________________________________________________________________________________________" << endl;
834  cout << "|du mkm |dv mkm |dw mkm |alpha mrad |beta mrad |gamma mrad |Comment" << endl;
835  outC << "struct data_t {" << endl;
836  outC << "\tInt_t barrel, layer, ladder, wafer, type;" << endl;
837  outC << "\tDouble_t u, Du, v, Dv, w, Dw, alpha, Dalpha, beta, Dbeta, gamma, Dgamma;" << endl;
838  outC << "\tChar_t *Comment;" << endl;
839  outC << "};" << endl;
840  outC << "data_t Data[] = {" << endl;
841  }
842  head++;
843  out << "__________________________________________________________________________________________________ " << endl;
844  cout << "__________________________________________________________________________________________________ " << endl;
845  ladder = l1 + i;
846  Int_t layer = 2*barrel - 1;
847  if (ladder%2 && barrel < 4) layer++;
848 
849  Double_t FitR[6], dFitR[6], LSFit[6], dLSFit[6];
850  memset (FitR, 0, 6*sizeof(Double_t));
851  memset (dFitR, 0, 6*sizeof(Double_t));
852  memset (LSFit, 0, 6*sizeof(Double_t));
853  memset (dLSFit, 0, 6*sizeof(Double_t));
854  TString line("");
855  TString lTitle("");
856  TString lineC("");
857  TH1D *LSFB = (TH1D *) gDirectory->Get(Form("LSFB%i",barrel));
858  if (LSFB) {
859  Double_t *array = LSFB->GetArray();
860  Int_t im = 1 + 28*(ladder-1);
861  Int_t is = im + 6;
862  TRVector AmX(6,array+im);// cout << "AmX " << AmX << endl;
863  TRSymMatrix S(6,array+is);// cout << "S " << S << endl;
864  TRSymMatrix SInv(S,TRArray::kInverted);// cout << "SInv " << SInv << endl;
865  TRVector X(SInv,TRArray::kSxA,AmX); //cout << "X " << X << endl;
866  TString line("");
867  for (Int_t l = 0; l < 6; l++) {
868  // if (l == 4) X(l) = -X(l); // switch sign
869  Double_t scale = 1e4;
870  if (l > 2) scale = 1e3;
871  LSFit[l] = scale*X(l);
872  dLSFit[l] = scale*TMath::Sqrt(SInv(l,l));
873  line += Form("|%7.2f+-%5.2f ",LSFit[l],dLSFit[l]);
874  if (SInv(l,l) > 0) {
875  FitR[l] += X(l)/SInv(l,l);
876  dFitR[l] += 1./SInv(l,l);
877  }
878  }
879  line += "|"; line += LSFB->GetName(); line += "/"; line += LSFB->GetTitle();// line += "\t"; line += f->GetName();
880  cout << line << endl;
881  out << line << endl;
882  }
883  for (Int_t j = firstH; j <= lastH; j++) {
884  Int_t Id = ladder + 100*(wafer + 10*layer);
885  h = (TH2 *) gDirectory->Get(Form("%s%04i",plotName[j],Id));
886  if (! h) continue;
887  h->SetMinimum(1);
888  h->SetXTitle(gDirectory->GetName());
889  Int_t ij = i + nx*(j-firstH) + 1;
890  c1->cd(ij)->SetLogz(1);
891  TH1 *fit = 0;
892  TProfile *prof = h->ProfileX();
893  prof->SetMarkerStyle(24);
894  prof->SetMarkerColor(6);
895  TH1 *sp = h->ProjectionY("_py",-1,-1,"e");
896  sp->Fit("gaus","q");
897  Double_t Mu = 0;
898  Double_t dMu = 0;
899  TF1 *gaus = sp->GetFunction("gaus");
900  gp->SetParameters(TMath::Log(gaus->GetParameter(0)),gaus->GetParameter(1),TMath::Abs(gaus->GetParameter(2)),0.);
901  sp->Fit(gp,"q");
902  Mu = sp->GetFunction("gp")->GetParameter(1);
903  dMu = sp->GetFunction("gp")->GetParError(1);
904 
905  Double_t *params = gp->GetParameters();
906  params[0] -= TMath::Log(100.);
907  gp->SetParameters(params);
908  // SlicesYFit(h,0,0,10,"qnig3"); // g3
909  SlicesYFit(h,0,0,10,"qni"); // g3
910  fit = (TH1 *) gDirectory->Get(Form("%s_1",h->GetName()));
911  // TH1 *sig = (TH1 *) gDirectory->Get(Form("%s_2",h->GetName()));
912  // TH1 *gra = (TH1 *) gDirectory->Get(Form("%s_3",h->GetName()));
913  Double_t slope = 0;
914  Double_t dslope = 0;
915  TLegend *leg = new TLegend(0.1,0.2,0.9,0.3,"");
916  leg->SetTextSize(0.025);
917  if (fit) {
918  fit->SetTitle(h->GetTitle());
919  fit->SetMarkerStyle(20);
920  fit->SetMarkerColor(1);
921  fit->SetMaximum(0.2);
922  fit->SetMinimum(-.2);
923  fit->SetStats(1);
924  Double_t zmax = 99;
925  FitPolN(fit,-zmax,zmax);
926  // fit->Fit("PolN","eqr","",-zmax,zmax);
927  TF1 *pol1 = fit->GetFunction("PolN");
928  if (! pol1 ) goto endhLoop;
929  Double_t prob = pol1->GetProb();
930 // if (prob <= 1e-3) {
931 // fit = prof;
932 // FitPolN(fit,-zmax,zmax);
933 // pol1 = fit->GetFunction("PolN");
934 // prob = pol1->GetProb();
935 // }
936  if (prob >= 0) {
937  Mu = pol1->GetParameter(0);
938  dMu = pol1->GetParError(0);
939  if (dMu > 99.99e-4) dMu= 99.99e-4;
940  slope = pol1->GetParameter(1);
941  dslope = pol1->GetParError(1);
942  if (dslope > 99.99e-3) dslope = 99.99e-3;
943  } else {
944  Mu = slope = 0;
945  dMu = dslope = 0;
946  }
947  static const Char_t *duv[2] = {"du", "dv"};
948  TString Name(h->GetName());
949  TString Title(h->GetTitle());
950  line = "";
951  lTitle = "";
952  lineC = Form("\t{%1i,%1i,%2i,%2i,%2i",barrel, layer, ladder, wafer, j);
953  for (Int_t m = 0; m < 2; m++) {
954  if (Name.BeginsWith(duv[m]) && ! Name.Contains("Over") && dMu > 0) {
955  Double_t mu = -1e4*Mu;
956  Double_t dmu = 1e4*dMu;
957  line += Form("|%7.2f+-%5.2f",mu,dmu);
958  lineC += Form(",%7.2f,%5.2f",mu,dmu);
959  lTitle += Form("%s = %7.2f +- %5.2f (mkm)", duv[m], -mu,dmu);
960  Double_t dev = mu - LSFit[m];
961  Double_t sdev = TMath::Sqrt(dmu*dmu+dLSFit[m]*dLSFit[m]);
962  if (dLSFit[m] == 0 || sdev > 0 && TMath::Abs(dev/sdev) < nSigMax) {
963  Double_t dMu2 = dMu*dMu;
964  FitR[m] += -Mu/dMu2;
965  dFitR[m] += 1./dMu2;
966  line += "A";
967  } else line +="R";
968  }
969  else {
970  line += Form("| ");
971  lineC += ", 0,-9.99";
972  }
973  }
974  Int_t index = Title.Index("=>");
975  TString tag("");
976  if (index >= 0) {
977  index = index+2;
978  static TString separator("[^ ;,]+");
979  TString t(Title.Data()+index);
980  TObjArray *array = t.Tokenize(separator);
981  tag = ((TObjString *) array->At(0))->GetString();
982  delete array;
983  }
984  for (Int_t m = 2; m < 6; m++) {
985  Double_t dslope2, scale;
986  Double_t mu, dmu;
987  if (dslope <= 0) goto Empty;
988  scale = 1e3;
989  switch (m) {
990  case 2:
991  if (! tag.Contains("dw")) goto Empty;
992  scale = 1e4;
993  mu = scale*slope;
994  dmu = scale*dslope;
995  line += Form("|%7.2f+-%5.2f",mu,dmu);
996  lineC += Form(",%7.2f,%5.2f",mu,dmu);
997  lTitle += Form(" dw = %7.2f +- %5.2f (mkm)", mu, dmu);
998  break;
999  case 3:
1000  if (! tag.Contains("alpha")) goto Empty;
1001  // if (tag.Contains("-alpha")) slope = - slope;
1002  mu = scale*slope;
1003  dmu = scale*dslope;
1004  line += Form("|%7.2f+-%5.2f",mu,dmu);
1005  lineC += Form(",%7.2f,%5.2f",mu,dmu);
1006  lTitle += Form(" alpha = %7.2f +- %5.2f (mrad)", mu, dmu);
1007  break;
1008  case 4:
1009  if (! tag.Contains("beta")) goto Empty;
1010  // if (! tag.Contains("-beta")) slope = - slope;
1011  mu = scale*slope;
1012  dmu = scale*dslope;
1013  line += Form("|%7.2f+-%5.2f", mu,dmu);
1014  lineC += Form(",%7.2f,%5.2f", mu,dmu);
1015  lTitle += Form(" beta = %7.2f +- %5.2f (mrad)", mu, dmu);
1016  break;
1017  case 5:
1018  if (! tag.Contains("gamma")) goto Empty;
1019  // if (tag.Contains("-gamma")) slope = - slope;
1020  mu = scale*slope;
1021  dmu = scale*dslope;
1022  line += Form("|%7.2f+-%5.2f", mu,dmu);
1023  lineC += Form(",%7.2f,%5.2f", mu,dmu);
1024  lTitle += Form(" gamma = %7.2f +- %5.2f (mrad)", mu, dmu);
1025  break;
1026  default:
1027  goto Empty;
1028  };
1029  if (! Name.Contains("duOvertuP") && dslope > 0) {
1030  Double_t dev = mu - LSFit[m];
1031  Double_t sdev = TMath::Sqrt(dmu*dmu+dLSFit[m]*dLSFit[m]);
1032  if (dLSFit[m] == 0 || sdev > 0 && TMath::Abs(dev/sdev) < nSigMax) {
1033  dslope2 = dslope*dslope;
1034  FitR[m] += slope/dslope2;
1035  dFitR[m] += 1./dslope2;
1036  line += "A";
1037  } else line +="R";
1038  } else line +="R";
1039  continue;
1040  Empty:
1041  line += Form("| ");
1042  lineC += ", 0,-9.99";
1043  }
1044  lTitle += Form(" prob = %5.3f",prob);
1045  leg->AddEntry(pol1,lTitle);
1046  line += "|"; line += fit->GetName(); line += "/"; line += h->GetTitle();
1047  lineC += ",\""; lineC += fit->GetName(); lineC += "\"},";
1048  cout << line << endl;
1049  out << line << endl;
1050  outC << lineC << endl;
1051  }
1052  endhLoop:
1053 #if 0
1054  if (h) h->DrawCopy("colz");
1055  if (prof) prof->DrawCopy("same");
1056 #else
1057  if (h) h->Draw("colz");
1058  if (prof) prof->Draw("same");
1059 #endif
1060  if (fit) {
1061 #if 0
1062  fit->DrawCopy("same");
1063 #else
1064  fit->Draw("same");
1065 #endif
1066  TF1 *pol1 = fit->GetFunction("PolN");
1067  if (pol1) {pol1->SetLineColor(2); pol1->Draw("same");}
1068  }
1069  leg->Draw();
1070  }
1071  line = ""; lineC = "";
1072  lineC = Form("\t{%1i,%1i,%2i,%2i,%2i",barrel, layer, ladder, wafer, -1);
1073  for (Int_t m = 0; m < 6; m++) {
1074  if (dFitR[m] > 0) {
1075  Double_t scale = 1e4;
1076  if (m > 2) scale = 1e3;
1077  FitR[m] = scale*FitR[m]/dFitR[m];
1078  dFitR[m] = scale/TMath::Sqrt(dFitR[m]);
1079  line += Form("|%7.2f+-%5.2f ", FitR[m],dFitR[m]);
1080  lineC += Form(",%7.2f,%5.2f", FitR[m],dFitR[m]);
1081  } else {
1082  line += Form("| ");
1083  lineC += ", 0,-9.99";
1084  }
1085  }
1086  line += "| Average for "; line += Form("barrel %1i,layer %1i,ladder %2i,wafer %2i",barrel, layer, ladder, wafer);
1087  lineC += ",\"Average\"},";
1088  cout << line << endl;
1089  out << line << endl;
1090  outC << lineC << endl;
1091  }
1092  out.close();
1093  outC.close();
1094 }
1095 //________________________________________________________________________________
1096 //void TDrawD(const Char_t *tag="duuH", Int_t barrel = 1, Int_t ladder = 0, Int_t wafer = 0) {// fit drift velocity
1097 void TDrawD(const Char_t *tag="duuH", Int_t barrel = 1, Int_t ladder = 0, Int_t wafer = 0) {// fit drift velocity
1098  static const Int_t NlPerBarrel[3] = {8, 12, 16};
1099  static const Int_t NwPerBarrel[3] = {4, 6, 7};
1100  if (barrel < 1 || barrel > 3) {cout << "Wrong barrel no. " << barrel << endl; return;}
1101  Int_t NL = NlPerBarrel[barrel-1];
1102  Int_t l1 = 1;
1103  Int_t l2 = NL;
1104  if (ladder) {l1 = l2 = ladder;}
1105  Int_t nx = l2 - l1 + 1;
1106  Int_t NW = NwPerBarrel[barrel-1];
1107  Int_t w1 = 1;
1108  Int_t w2 = NW;
1109  if (wafer) {w1 = w2 = wafer;}
1110  Int_t ny = w2 - w1 + 1;
1111  Int_t scaleX = 60; //600/nx;
1112  Int_t scaleY = 40; //800/ny;
1113  // Int_t scale = TMath::Min(scaleX,scaleY);
1114 // gStyle->SetPadBottomMargin( 0.01);
1115 // gStyle->SetPadLeftMargin( 0.01);
1116 // gStyle->SetPadRightMargin( 0.01);
1117 // gStyle->SetPadBottomMargin( 0.01);
1118  TH1D *fitPN[2];
1119  c1 = new TCanvas(Form("Drift_Barrel_%i",barrel),Form("Barrel %i, Ladder %i, Wafer %i",barrel,ladder,wafer) ,10,10,10+scaleX*nx,10+scaleY*ny);
1120  cout << "nx/ny = " << nx << "/" << ny << endl;
1121  if (nx > 1 || ny > 1) c1->Divide(nx,ny);
1122  ofstream out;
1123  ofstream outC;
1124  TString Out("Results.");
1125  Out += Form("DriftBarrel_%i",barrel);
1126  Out += gSystem->BaseName(gDirectory->GetName());
1127  Out.ReplaceAll(".root","");
1128  TDatime t;
1129  // Out += t.AsString();
1130  Out.ReplaceAll(" ","");
1131  Out += ".h";
1132  if (gSystem->AccessPathName(Out)) outC.open(Out, ios::out); //"Results.list",ios::out | ios::app);
1133  else outC.open(Out, ios::app);
1134  TH2 *h = 0;
1135  Int_t head = 0;
1136  Int_t NPol1 = NPOL;// first parameter reserved for degree of polynomial
1137  if (TString(tag) == "duvH") {NPol1 = 12 - NPOL;} // only 3 power for anodes
1138  const Int_t NPMax = NPol1;
1139  Double_t params[NPol1];
1140  Double_t dparams[NPol1];
1141  TLegend *leg = 0;
1142  Double_t Ymnx[2], Y;
1143  Int_t p10;
1144  Double_t xFmin[2] = {-1.08, 0.12};
1145  Double_t xFmax[2] = {-0.12, 1.08};
1146  for (Int_t i = 0; i < nx; i++) {
1147  if (! head) {
1148  outC << "struct data_t {" << endl;
1149  outC << "\tInt_t type, idx, nrows, barrel, layer, ladder, wafer, hybrid, Npar;" << endl;
1150  outC << "\tDouble_t ";
1151  for (Int_t i = 0; i < NPol1-1; i++) {
1152  outC << "v" << i;
1153  if (i < NPol1 - 2) {outC << ", ";}
1154  else {outC << ";";}
1155  }
1156  outC << endl;
1157  // outC << "\tDouble_t param[" << NPol1 - 1 << "];"<< endl;
1158  // outC << "\tDouble_t dparam[" << NPol1 - 1 << "];"<< endl;
1159  outC << "\tChar_t Comment[10];" << endl;
1160  outC << "};" << endl;
1161  outC << "data_t Data[] = {" << endl;
1162  }
1163  head++;
1164  ladder = l1 + i;
1165  Int_t layer = 2*barrel - 1;
1166  if (ladder%2) layer++;
1167  TString lTitle("");
1168  TString lineC("");
1169  TF1 *f[2];
1170  static const Char_t *funName[2] = {"fSTchebN","fSTchebP"};
1171  for (Int_t k = 0; k < 2; k++) {
1172  f[k] = (TF1 *) gROOT->GetFunction(funName[k]);
1173  if (! f[k]) {
1174  if (k == 0) f[k] = new TF1(funName[k],STchebN, xFmin[k], xFmax[k],NPol1);
1175  else f[k] = new TF1(funName[k],STchebP, xFmin[k], xFmax[k],NPol1);
1176  }
1177  f[k]->FixParameter(0,1); // N == 1
1178  f[k]->SetParameter(1,0);
1179  f[k]->SetParameter(2,0);
1180  for (Int_t i = 3; i <= NPol1; i++) f[k]->FixParameter(i,0);
1181  f[k]->SetLineColor(2*k+1);
1182  // f[k]->Print();
1183  }
1184  TH1D *py, *fit;
1185  TProfile *prof;
1186  for (Int_t j = 0; j < ny; j++) {
1187  Int_t ij = i + nx*j + 1;
1188  c1->cd(ij)->SetLogz(1);
1189  wafer = w1 + j;
1190  Int_t Id = ladder + 100*(wafer + 10*layer);
1191  h = (TH2 *) gDirectory->Get(Form("%s%04i",tag,Id));
1192  if (!h) continue;
1193  h->SetMinimum(1);
1194  prof = 0;
1195  fit = 0;
1196  leg = 0;
1197  if (! h) continue;
1198  TAxis *yax = h->GetYaxis();
1199  Double_t ymin = yax->GetXmin();
1200  Double_t ymax = yax->GetXmax();
1201  Int_t nybins = yax->GetNbins();
1202  if (h->GetEntries() < 100) continue;;
1203  h->SetXTitle(gDirectory->GetName());
1204  if ( h->GetRMS(2) > 0.9*(ymax-ymin)/TMath::Sqrt(12.)) goto ENDL;
1205  prof = h->ProfileX();
1206  prof->SetMarkerStyle(24);
1207  prof->SetMarkerColor(6);
1208  // SlicesYFit(h,0,0,10,"qnig3");
1209  SlicesYFit(h,0,0,10,"qni");
1210  fit = (TH1D *) gDirectory->Get(Form("%s_1",h->GetName()));
1211  if (! fit) goto ENDL;
1212  Ymnx[0] = fit->GetMinimum();
1213  Ymnx[1] = fit->GetMaximum();
1214  for (Int_t l = 0; l < 2; l++) {
1215  p10 = 0;
1216  if (TMath::Abs(Ymnx[l]) > 0) p10 = (Int_t) TMath::Log10(TMath::Abs(Ymnx[l]));
1217  Y = TMath::Power(10.,p10);
1218  p10 = (Int_t) (TMath::Abs(Ymnx[l]/Y)) + 1;
1219  Ymnx[l] = TMath::Sign(p10*Y,Ymnx[l]);
1220  }
1221  // yax->SetRange(yax->FindBin(Ymnx[0]),yax->FindBin(Ymnx[1]));
1222  yax->SetRange(yax->FindBin(-.25),yax->FindBin(0.25));
1223  fit->SetMarkerStyle(20);
1224  fit->SetMarkerColor(1);
1225  fitPN[0] = new TH1D(*fit);
1226  fitPN[1] = new TH1D(*fit);
1227  for (Int_t k = 0; k < 2 ; k++) {
1228  if (k == 0) {py = h->ProjectionY("_py1",1,nybins/2-1,"e");}
1229  else {py = h->ProjectionY("_py2",nybins/2+2,nybins,"e");}
1230  if (py->GetEntries() > 100 || py->GetRMS() < 0.9*(ymax-ymin)/TMath::Sqrt(12.)) {
1231 #ifdef __PROB_SELECTION__
1232  Double_t oldProb = 0;
1233 #else
1234  Double_t chi2Old = 1e20;
1235 #endif
1236  Int_t na = 0;
1237  Int_t np = 0;
1238  for (Int_t n = 0; n < NPMax-1; n++) {
1239  f[k]->FixParameter(0,n);
1240  f[k]->ReleaseParameter(n+1);
1241  fitPN[k]->Fit(f[k],"erq","",xFmin[k],xFmax[k]);
1242  if (f[k]->GetNumberFitPoints() < 10) break;
1243 #ifdef __PROB_SELECTION__
1244  Double_t prob = f[k]->GetProb();
1245  if (prob > 1e-3) {
1246  if (oldProb < 1e-3) {oldProb = prob;}
1247  else {
1248  if (prob < 2*oldProb) {
1249  f[k]->FixParameter(0,n-1);
1250  f[k]->FixParameter(n+1,0);
1251  break;
1252  } else oldProb = prob;
1253  }
1254  }
1255 #else
1256  Double_t chi2 = f[k]->GetChisquare();
1257  Double_t dChi2 = chi2Old - chi2;
1258  Double_t par = f[k]->GetParameter(n+1);
1259  Double_t dpar = f[k]->GetParError(n+1);
1260  if (n == 0 || TMath::Abs(par) > 2*dpar && dChi2 > chi2Old/(NPMax - 1 - na)) {chi2Old = chi2; na++; np = n;} // accepted
1261  else f[k]->FixParameter(n+1,0); // rejected
1262 #endif
1263  }
1264  f[k]->FixParameter(0,np);
1265  fitPN[k]->Fit(f[k],"erq");
1266  c1->Update();
1267 
1268  if (f[k]->GetNumberFitPoints() >= 10) {
1269  Int_t NPar = (Int_t) f[k]->GetParameter(0);
1270  memset (params, 0, NPol1*sizeof(Double_t));
1271  memset (dparams, 0, NPol1*sizeof(Double_t));
1272  for (Int_t l = 1; l <= NPar+1; l++) {
1273  params[l] = f[k]->GetParameter(l);
1274  dparams[l] = f[k]->GetParError(l);
1275  }
1276 
1277  lTitle = Form("Shift = %8.2f +- %4.2f (mkm) Slope = %8.2f +- %4.2f (mrad) ",1e4*params[1],1e4*dparams[1],1e3*params[2],1e3*dparams[3]);
1278  lTitle += Form(" N = %i prob = %4.3f", NPar, f[k]->GetProb());
1279 
1280  if (! leg) {
1281  leg = new TLegend(0.1,0.2,0.9,0.3,"");
1282  leg->SetTextSize(0.020);
1283  }
1284  leg->AddEntry(f[k],lTitle);
1285  static const Int_t type = 2, idx = 0, nrows = 0;
1286  // lineC = Form("\t{%1i,%1i,%1i,%1i,%1i,%2i,%2i,%2i,%2i,\n{",type,idx,nrows,barrel, layer, ladder, wafer, k+1, NPar);
1287  lineC = Form("\t{%1i,%1i,%1i,%1i,%1i,%2i,%2i,%2i,%2i,",type,idx,nrows,barrel, layer, ladder, wafer, k+1, NPar);
1288  // for (Int_t l = 1; l < NPol1; l++) {lineC += Form("%8.5f",params[l]); if (l != NPol1 - 1) lineC += ",";} lineC += "},\n{";
1289  for (Int_t l = 1; l < NPol1; l++) {lineC += Form("%8.5f",params[l]); if (l != NPol1 - 1) lineC += ",";} lineC += "";
1290  // for (Int_t l = 1; l < NPol1; l++) {lineC += Form("%8.5f",dparams[l]); if (l != NPol1 - 1) lineC += ",";} lineC += "},\n";
1291  lineC += ",\""; lineC += h->GetName(); lineC += "\"},";
1292  outC << lineC << endl;
1293  }
1294  }
1295  }
1296  ENDL:
1297 #if 0
1298  h->DrawCopy("colz");
1299  if (prof) prof->DrawCopy("same");
1300  if (fit) {
1301 
1302  for (Int_t k = 0; k < 2; k++) {fitPN[k]->DrawCopy("same"); f[k]->DrawCopy("same");}
1303  }
1304 #else
1305  h->Draw("colz");
1306  if (prof) prof->Draw("same");
1307  if (fit) {
1308  for (Int_t k = 0; k < 2; k++) {fitPN[k]->Draw("same");}
1309  // f[k]->Draw("same");}
1310  }
1311 #endif
1312  if (leg) leg->Draw();
1313  }
1314  }
1315  outC.close();
1316 }
1317 //________________________________________________________________________________
1318 void TDraw(Int_t k = 0) {
1319  if (k == 0) {TDrawG(); return;}
1320  if (k == -3) {
1321  for (Int_t i = 3; i >= 1; i--) TDrawL(-1,i);
1322  return;
1323  }
1324  if (k < 5) {TDrawL(-1,k); return;}
1325  if (k == 5) {
1326  for (Int_t i = 4; i >= 1; i--) TDrawL(-1,i);
1327  return;
1328  }
1329  if (k == 10) {
1330  for (Int_t i = 1; i <= 3; i++) TDrawD("duuH",i);
1331  return;
1332  }
1333  if (k == 20) {
1334  for (Int_t i = 1; i <= 3; i++) TDrawD("duvH",i);
1335  return;
1336  }
1337 #if 0
1338  TInterpreter::EErrorCode error = TInterpreter::kNoError;
1339  TDrawL(-1,4); gInterpreter->ProcessLine(Form("TQtCanvas2Html d1((TPad *) %p,900,650);",c1), &error);
1340  for (Int_t i = 1; i <= 3; i++) {
1341  TDrawD("duuP",i); gInterpreter->ProcessLine(Form("TQtCanvas2Html d((TPad *) %p,900,650);",c1), &error);
1342  }
1343 #endif
1344 }
Rot_t
Definition: TDraw.C:86
Generated by   doxygen 1.8.5