StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
dEdxFit.C
1 /*
2  root.exe -q -b lBichsel.C pionMIP.root 'dEdxFit.C+("SecRow3C","GF")'
3 */
4 #if !defined(__CINT__)
5 // code that should be seen ONLY by the compiler
6 #else
7 #if !defined(__CINT__) || defined(__MAKECINT__)
8 // code that should be seen by the compiler AND rootcint
9 #else
10 // code that should always be seen
11 #endif
12 #endif
13 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,34,18)
14 //#define __USE_ROOFIT__
15 #endif
16 //________________________________________________________________________________
17 #if !defined(__CINT__) || defined(__MAKECINT__)
18 #include "Riostream.h"
19 #include <stdio.h>
20 #include "TROOT.h"
21 #include "TSystem.h"
22 #include "TMath.h"
23 #include "TH1.h"
24 #include "TH2.h"
25 #include "TH3.h"
26 #include "THnSparse.h"
27 #include "TStyle.h"
28 #include "TF1.h"
29 #include "TProfile.h"
30 #include "TTree.h"
31 #include "TChain.h"
32 #include "TFile.h"
33 #include "TNtuple.h"
34 #include "TFitResult.h"
35 #include "TCanvas.h"
36 #include "TFileSet.h"
37 #include "TDataSetIter.h"
38 #include "TDataSet.h"
39 #include "TClassTable.h"
40 //#include "DeDxTree.C"
41 #include "TMinuit.h"
42 #include "TSpectrum.h"
43 #include "StBichsel/Bichsel.h"
44 #include "StBichsel/StdEdxModel.h"
45 #include "TString.h"
46 #include "TLine.h"
47 #include "TText.h"
48 #include "TList.h"
49 #include "TPolyMarker.h"
50 #include "TKey.h"
51 #include "TLegend.h"
52 #ifdef __USE_ROOFIT__
53 #include "RooRealVar.h"
54 #include "RooDataSet.h"
55 #include "RooGaussian.h"
56 #include "RooFFTConvPdf.h"
57 #include "RooPlot.h"
58 #include "RooCFunction1Binding.h"
59 #include "RooCFunction3Binding.h"
60 #include "RooTFnBinding.h"
61 #include "RooDataHist.h"
62 #include "RooAbsPdf.h"
63 #include "RooRealProxy.h"
64 #include "RooFit.h"
65 #include "RooRandom.h"
66 #include "RooFitResult.h"
67 #include "RooWorkspace.h"
68 using namespace RooFit ;
69 #endif /* __USE_ROOFIT__ */
70 #include "TObjectTable.h"
71 #else
72 class TMinuit;
73 class TF1;
74 class TH1F;
75 class TH2F;
76 class TH3F;
77 class TProfile;
78 class TH2D;
79 class TCanvas;
80 class TSpectrum;
81 class TSystem;
82 class Bichsel;
83 // Refer to a class implemented in libRooFit to force its loading
84 // via the autoloader.
85 #ifdef __USE_ROOFIT__
86 class Roo2DKeysPdf;
87 #endif /* __USE_ROOFIT__ */
88 #endif
89 #include "Names.h"
90 using namespace std;
91 const Int_t NH = NHYPS/2;
92 Int_t N = 0;
93 const Int_t noPoints = 1000;
94 Double_t X[noPoints];
95 Double_t dX[noPoints];
96 Double_t Nu[noPoints];
97 Double_t Mu[noPoints];
98 Double_t dMu[noPoints];
99 Double_t Sigma[noPoints];
100 Double_t dSigma[noPoints];
101 TCanvas *canvas = 0;
102 Double_t Xlog10bg, Ylog2dx, Z;
103 static TFile *newf = 0;
104 static TFile *fOut = 0;
105 static TNtuple *FitP = 0;
106 static TH1 *projNs[5];
107 // peak postion at p = 0.475 GeV/c wrt pion
108 // Z pion
110 struct peak_t {Double_t peak, sigma, mass; const Char_t *Name;};
111 // mean RMS
112 static const peak_t Peaks[6] = {
113 // {0. , 0., 0.13956995, "pion"}, // pion
114 // {1.425822, 0.101693, 0.93827231, "proton"}, // proton - pion
115 // {0.565455, 0.061626, 0.493677, "kaon" }, // Kaon - pi
116 // {0.424916, 0.004081, 0.51099907e-3,"e"}, // e - pi
117 // {2.655586, 0.123754, 1.875613, "d"}, // d - pi
118 // {0.004178, 0.002484, 0.105658, "mu"}};// mu - pi
119  // 06/25/10
120  { 0. , 0., 0.13956995, "pion"}, // pion
121  { 1.42574, 0.101741, 0.938272, "proton"},
122  { 0.565411, 0.0616611, 0.493677, "kaon"},
123  { 0.424919, 0.00408318, 0.000510999, "e"},
124  { 2.65548, 0.123809, 1.87561, "deuteron"},
125  { 0.000717144, 0.00490783, 0.105658, "mu"}};
126 
127 Bichsel *gBichsel = 0;
128 //#define PRINT 1
129 TF1 *func = 0;
130 TF1 *LandauF = 0;
131 class St_TpcSecRowCor;
132 #ifdef __USE_ROOFIT__
133 Double_t landauZ(Double_t *x, Double_t *par) {
134  Double_t xd = x[0];
135  Double_t meand = par[0];
136  Double_t sigmad = par[1];
137  Double_t mpshift = -0.22278298; // 1.0844535734 + 0.61417; // LandauZ maximum location
138  // MP shift correction
139  Double_t mpc = TMath::Exp(meand) - mpshift * sigmad;
140  Double_t xx = TMath::Exp(xd);
141  return xx*TMath::Landau(xx,mpc,sigmad);
142 }
143 
144 class RooRealVar;
145 
146 class RooLandauZ : public RooAbsPdf {
147 public:
148  RooLandauZ() {} ;
149  RooLandauZ(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma);
150  RooLandauZ(const RooLandauZ& other, const char* name=0);
151  virtual TObject* clone(const char* newname) const { return new RooLandauZ(*this,newname); }
152  inline virtual ~RooLandauZ() { }
153 
154  Int_t getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t staticInitOK=kTRUE) const;
155  void generateEvent(Int_t code);
156 
157 protected:
158 
159  RooRealProxy x ;
160  RooRealProxy mean ;
161  RooRealProxy sigma ;
162 
163  Double_t evaluate() const ;
164 
165 private:
166 
167  ClassDef(RooLandauZ,1) // LandauZ Distribution PDF
168 };
169 
170 ClassImp(RooLandauZ)
171 
172 
173 //_____________________________________________________________________________
174 RooLandauZ::RooLandauZ(const char *name, const char *title, RooAbsReal& _x, RooAbsReal& _mean, RooAbsReal& _sigma) :
175  RooAbsPdf(name,title),
176  x("x","Dependent",this,_x),
177  mean("mean","Mean",this,_mean),
178  sigma("sigma","Width",this,_sigma)
179 {
180 }
181 
182 
183 //_____________________________________________________________________________
184 RooLandauZ::RooLandauZ(const RooLandauZ& other, const char* name) :
185  RooAbsPdf(other,name),
186  x("x",this,other.x),
187  mean("mean",this,other.mean),
188  sigma("sigma",this,other.sigma)
189 {
190 }
191 
192 
193 //_____________________________________________________________________________
194 Double_t RooLandauZ::evaluate() const
195 {
196  // return TMath::Landau(x, mean, sigma);
197  const Double_t par[2] = {mean, sigma};
198  Double_t xd = x;
199  return landauZ(&xd,(Double_t *) par);
200  // return xx*TMath::Landau(xx,mpc,sigmad);
201  // return xx*TMath::Landau(xx,mpc,sigmad) / sigmad;
202 
203 }
204 
205 
206 //_____________________________________________________________________________
207 Int_t RooLandauZ::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
208 {
209  if (matchArgs(directVars,generateVars,x)) return 1 ;
210  return 0 ;
211 }
212 
213 
214 //_____________________________________________________________________________
215 void RooLandauZ::generateEvent(Int_t code)
216 {
217  assert(code==1) ;
218  Double_t xgen ;
219  Double_t mpshift = -0.22278298; // 1.0844535734 + 0.61417; // Landau maximum location
220  // MP shift correction
221  Double_t meand = mean;
222  Double_t sigmad = sigma;
223  Double_t mpc = TMath::Exp(meand) - mpshift * sigmad;
224  while(1) {
225  xgen = RooRandom::randomGenerator()->Landau(mpc,sigma);
226  if (xgen <= 0) continue;
227  Double_t xgl = TMath::Log(xgen) ;
228  if (xgl<x.max() && xgl>x.min()) {
229  x = xgl;
230  break;
231  }
232  }
233  return;
234 }
235 //_______________________________________________________________________________
237 //--------------------------B E G I N N I N G ------------------------------//
238 //------------------------L A N D A U Z 5 function ------------------//
239 //---------------------------------------------------------------//
240 Double_t landauZ5(Double_t *x, Double_t *par) {
241  // Double_t Norm = par[0];
242  Double_t xd = x[0];
243  Double_t meand_pi = par[1];
244  Double_t meand_pr = par[1] + Peaks[1].peak;
245  Double_t sigmad_pi = par[8];
246  Double_t sigmad_pr = sigmad_pi;
247  Double_t meand_k = par[1] + Peaks[2].peak;
248  Double_t sigmad_k = sigmad_pi;
249  Double_t meand_el = par[1] + Peaks[3].peak;
250  Double_t sigmad_el = sigmad_pi;
251  Double_t meand_d = par[1] + Peaks[4].peak;
252  Double_t sigmad_d = sigmad_pi;
253  Double_t mpshift = -0.22278298; // 1.0844535734 + 0.61417; // LandauZ maximum location
254  Double_t frac[5];
255  frac[1] = TMath::Power(TMath::Sin(par[3]),2);
256  frac[2] = TMath::Power(TMath::Sin(par[4]),2);
257  // frac[3] = 0.;
258  // frac[4] = 0.;
259  frac[3] = TMath::Power(TMath::Sin(par[5]),2);
260  frac[4] = TMath::Power(TMath::Sin(par[6]),2);
261 
262  if( ((par[4]==0.)&&(par[5]==0.)&&(par[6]==0.)&&(par[3]!=0.)) ||
263  ((par[3]==0.)&&(par[5]==0.)&&(par[6]==0.)&&(par[4]!=0.)) ||
264  ((par[3]==0.)&&(par[4]==0.)&&(par[6]==0.)&&(par[5]!=0.)) ||
265  ((par[3]==0.)&&(par[4]==0.)&&(par[5]==0.)&&(par[6]!=0.)) )
266  {frac[0]=0.;
267  }
268 
269  else {frac[0] = 1. - frac[1] - frac[2] - frac[3] - frac[4];}
270  // MP shift correction
271  Double_t mpc_pi = TMath::Exp(meand_pi) - mpshift * sigmad_pi;
272  Double_t mpc_pr = TMath::Exp(meand_pr) - mpshift * sigmad_pr;
273  Double_t mpc_k = TMath::Exp(meand_k) - mpshift * sigmad_k;
274  Double_t mpc_el = TMath::Exp(meand_el) - mpshift * sigmad_el;
275  Double_t mpc_d = TMath::Exp(meand_d) - mpshift * sigmad_d;
276 
277  Double_t xx = TMath::Exp(xd);
278 
279  return (xx*TMath::Landau(xx,mpc_pi,sigmad_pi)*frac[0] +
280  xx*TMath::Landau(xx,mpc_pr,sigmad_pr)*frac[1] +
281  xx*TMath::Landau(xx,mpc_k, sigmad_k) *frac[2] +
282  xx*TMath::Landau(xx,mpc_el,sigmad_el)*frac[3] +
283  xx*TMath::Landau(xx,mpc_d, sigmad_d) *frac[4])*par[7];
284 }
285 
286 
287 //________________________________________________________________________________
288 
289 //-------------------------L A N D A U 5 C L A S S --------------------------------------------------------------------------------//
290 //-------------------------------------------------------------------------------------------------------------------------//
291 
292 class RooLandauZ5:public RooAbsPdf
293 {
294 
295 public:
296  RooLandauZ5(){}
297  RooLandauZ5(const char *name, const char *title, RooAbsReal& _x,
298  RooAbsReal& _norm, RooAbsReal& _mean_pi,
299  RooAbsReal& _sigma_pi, RooAbsReal& _sigma_pr,
300  RooAbsReal& _sigma_k, RooAbsReal& _sigma_el,
301  RooAbsReal& _sigma_d, RooAbsReal& _total, RooAbsReal& _width);
302  RooLandauZ5(const RooLandauZ5& other, const char* name=0);
303  virtual TObject *clone(const char *newname) const
304  { return new RooLandauZ5(*this,newname); }
305  inline virtual ~RooLandauZ5() {}
306 
307 protected:
308  RooRealProxy x;
309  RooRealProxy norm;
310  RooRealProxy mu;
311  RooRealProxy sigma_pi;
312  RooRealProxy sigma_pr;
313  RooRealProxy sigma_k;
314  RooRealProxy sigma_el;
315  RooRealProxy sigma_d;
316  RooRealProxy total;
317  RooRealProxy width;
318 
319  Double_t evaluate() const;
320 
321 private:
322  ClassDef(RooLandauZ5,1)
323 };
324 
325 ClassImp(RooLandauZ5)
326 
327 RooLandauZ5::RooLandauZ5(const char *name, const char *title, RooAbsReal& _x,
328  RooAbsReal& _norm, RooAbsReal& _mu,
329  RooAbsReal& _sigma_pi, RooAbsReal& _sigma_pr,
330  RooAbsReal& _sigma_k, RooAbsReal& _sigma_el,
331  RooAbsReal& _sigma_d, RooAbsReal& _total,
332  RooAbsReal& _width):
333  RooAbsPdf(name,title),
334  x("x","Dependent",this,_x),
335  norm("norm","Norm",this,_norm),
336  mu("mu","peak position",this,_mu),
337  sigma_pi("sigma_pi","Width pion",this,_sigma_pi),
338  sigma_pr("sigma_pr","Width proton",this,_sigma_pr),
339  sigma_k("sigma_k","Width kaon",this,_sigma_k),
340  sigma_el("sigma_el","Width electron",this,_sigma_el),
341  sigma_d("sigma_d","Width deutron",this,_sigma_d),
342  total("total","Total",this,_total),
343  width("width","Width",this,_width){}
344 
345 RooLandauZ5::RooLandauZ5(const RooLandauZ5& other, const char* name):
346  RooAbsPdf(other,name),
347  x("x",this,other.x),
348  norm("norm",this,other.norm),
349  mu("mu",this,other.mu),
350  sigma_pi("sigma_pi",this,other.sigma_pi),
351  sigma_pr("sigma_pr",this,other.sigma_pr),
352  sigma_k("sigma_k",this,other.sigma_k),
353  sigma_el("sigma_el",this,other.sigma_el),
354  sigma_d("sigma_d",this,other.sigma_d),
355  total("total",this,other.total),
356  width("width",this,other.width){}
357 
358 
359 Double_t RooLandauZ5::evaluate() const{
360  const Double_t par[10] = {norm,mu,sigma_pi,sigma_pr,sigma_k,sigma_el,sigma_d,total,width};
361  Double_t xd =x;
362  return landauZ5(&xd, (Double_t *)par);
363 }
364 
365 //---------------------------------------------------------------------//
366 //------Global variables needed to create TF1 for particle fraction ----//
367 TF1 *l5xg_func = 0, *func_pi=0, *func_pr=0, *func_k=0, *func_el=0, *func_de=0;
368 Double_t frac_pr=0.,frac_k=0.,frac_el=-0.,frac_de=-0., frac_pi=0.;
369 
370 //-------------Function for scaling TF1 ---------------//
371 Double_t func_lz5xg_mult(Double_t *x, Double_t *par) {
372  if (! l5xg_func) return 0;
373  Double_t xd = x[0];
374  Double_t mult = par[7];
375  Int_t Case = (int)par[9];
376  switch(Case){
377  case 0:
378  if (! func_pi) return 0;
379  return ((func_pi->Eval(xd))*mult*frac_pi);
380  break;
381  case 1:
382  if (! func_pr) return 0;
383  return ((func_pr->Eval(xd))*mult*frac_pr);
384  break;
385  case 2:
386  if (! func_k) return 0;
387  return ((func_k->Eval(xd))*mult*frac_k);
388  break;
389  case 3:
390  if (! func_el) return 0;
391  return ((func_el->Eval(xd))*mult*frac_el);
392  break;
393  case 4:
394  if (! func_de) return 0;
395  return ((func_de->Eval(xd))*mult*frac_de);
396  break;
397  case 5:
398  return ((l5xg_func->Eval(xd))*mult);
399  break;
400  default:
401  return ((l5xg_func->Eval(xd))*mult);
402  break;
403  }
404 }
405 
406 //--------------------Function for creating TF1 for particle fractions -----//
407 void Sep_func(RooFFTConvPdf *l5xg, RooRealVar *t, RooRealVar *norm, RooRealVar *mu, RooRealVar *sg, RooRealVar *fProton, RooRealVar *fKaon, RooRealVar *fElectron, RooRealVar *fDeuteron, RooRealVar *total, RooRealVar *width, Int_t i )
408 {
409  switch(i){
410  case 0://Pion//
411  if (func_pi) delete func_pi;
412  func_pi = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
413  func_pi->FixParameter(3,0.);
414  func_pi->FixParameter(4,0.);
415  func_pi->FixParameter(5,0.);
416  func_pi->FixParameter(6,0.);
417  break;
418  case 1://Proton//
419  if (func_pr) delete func_pr;
420  func_pr = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
421  frac_pr = TMath::Power(TMath::Sin(func_pr->GetParameter(i+2)),2);
422  func_pr->FixParameter(4,0.);
423  func_pr->FixParameter(5,0.);
424  func_pr->FixParameter(6,0.);
425  break;
426  case 2: //Kaon//
427  if (func_k) delete func_k;
428  func_k = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
429  frac_k = TMath::Power(TMath::Sin(func_k->GetParameter(i+2)),2);
430  func_k->FixParameter(3,0.);
431  func_k->FixParameter(5,0.);
432  func_k->FixParameter(6,0.);
433  break;
434  case 3: //Elektron//
435  if (func_el) delete func_el;
436  func_el = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
437  frac_el = TMath::Power(TMath::Sin(func_el->GetParameter(i+2)),2);
438  func_el->FixParameter(3,0.);
439  func_el->FixParameter(4,0.);
440  func_el->FixParameter(6,0.);
441  break;
442  case 4: //Deuteron//
443  if (func_de) delete func_de;
444  func_de = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
445  frac_de = TMath::Power(TMath::Sin(func_de->GetParameter(i+2)),2);
446  func_de->FixParameter(3,0.);
447  func_de->FixParameter(4,0.);
448  func_de->FixParameter(5,0.);
449  break;
450  default: break;
451  }
452 }
453 
454 //----------------------------------------------------------------------//
455 
456 TF1 *FitRL5(TH1 *hist, Bool_t outer = kFALSE)
457 {
458  if(!hist) return 0;
459 // cout<<"TEST (FUNKCJA): OUTER sector = "<<outer<<", "<<hist->GetName()<<endl;
460  l5xg_func = 0; func_pi=0; func_pr=0; func_k=0; func_el=0;func_de=0; frac_pr=0.;frac_k=0.;frac_el=-0.;frac_de=-0.; frac_pi=0.;
461 
462  //-------------------------chose Inner or Outer sector----------------//
463  Double_t sec_w=0.;
464  if(outer) sec_w = 0.0699;//cout<<"!!!!!!!!!!!!WIDTH = "<<sec_w;}
465  else sec_w = 0.0762;
466 
467  //-------------------------------Variables---------------------------//
468  static RooRealVar *t = 0, *norm = 0, *mu = 0, *sigma = 0, *fProton = 0, *fKaon=0, *fElectron = 0, *fDeuteron = 0, *total = 0, *width = 0;
469  static RooRealVar *mg = 0, *sg = 0;
470  static RooLandauZ5 *landauZ5 = 0;
471  static RooGaussian *gauss = 0;
472  static RooFFTConvPdf *l5xg = 0;
473  static TF1 *l5xg_mult = 0;
474  // cout<<"kasowanie"<<endl;
475  if (t) delete t;
476  t = new RooRealVar("t" ,"t" ,-2.,5.) ;
477  t->setBins(1000,"cache");
478  if (norm) delete norm;
479  norm = new RooRealVar("norm" ,"norm" ,0.) ;
480  if (mu) delete mu;
481  mu = new RooRealVar("mu" ,"mu" ,0.,-1.,1.) ;
482  if (sigma) delete sigma;
483  sigma = new RooRealVar("sigma" ,"sigma" ,0.) ;
484  if (fProton) delete fProton;
485  fProton = new RooRealVar("fProton" ,"fProton" ,0.,0.,1.57) ;
486  if (fKaon) delete fKaon;
487  fKaon = new RooRealVar("fKaon" ,"fKaon" ,0.,0.,1.57) ;
488  if (fElectron) delete fElectron;
489  fElectron = new RooRealVar("fElectron","fElectron",0.) ;
490  if (fDeuteron) delete fDeuteron;
491  fDeuteron = new RooRealVar("fDeuteron","fDeuteron",0.) ;
492  if (total) delete total;
493  total = new RooRealVar("total" ,"total" ,0.,0.,10.);
494  if (width) delete width;
495  width = new RooRealVar("width" ,"width" ,sec_w);//0.1,0.01,2.5) ;
496  *width = sec_w;
497 
498  if (mg) delete mg;
499  mg = new RooRealVar("mg","mg",0.) ;
500  if (sg) delete sg;
501  sg = new RooRealVar("sg","sg",0.25,0.01,10.) ;
502 
503  if (landauZ5) delete landauZ5;
504  landauZ5= new RooLandauZ5("landauZ5","landauZ5",*t,*norm,*mu,*sigma,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width) ;
505  if (gauss) delete gauss;
506  gauss = new RooGaussian ("gauss","gauss",*t,*mg,*sg) ;
507  if (l5xg) delete l5xg;
508  l5xg = new RooFFTConvPdf("l5xg","landauZ5 (X) gauss",*t,*landauZ5,*gauss) ;
509  //Import data
510  //--------------------------------------------------------------
511  RooDataHist data("data","data",*t,Import(*hist));
512 
513  //------- Fit to data -----------------------//
514  l5xg->fitTo(data,Save());
515  //----- Create TF1 -------------//
516  if (l5xg_func) delete l5xg_func;
517  l5xg_func = (TF1*)l5xg->asTF(RooArgList(*t), RooArgList(*norm,*mu,*sg,*fProton,*fKaon,*fElectron,*fDeuteron,*total,*width),*t );
518  l5xg_func->SetParError(0,norm->getError());
519  l5xg_func->SetParError(1,mu->getError());
520  l5xg_func->SetParError(2,sg->getError());
521  l5xg_func->SetParError(3,fProton->getError());
522  l5xg_func->SetParError(4,fKaon->getError());
523  l5xg_func->SetParError(5,fElectron->getError());
524  l5xg_func->SetParError(6,fDeuteron->getError());
525  l5xg_func->SetParError(7,total->getError());
526  l5xg_func->SetParError(8,width->getError());
527 
528  //----- Scale TF1 fit function ------//
529  Double_t mult_ev = hist->Integral()*hist->GetBinWidth(5);
530  if (l5xg_mult) delete l5xg_mult;
531  l5xg_mult = new TF1(Form("mult_%s",hist->GetName()),func_lz5xg_mult,-2.,5.,11);
532  l5xg_mult->SetParent(hist);
533  l5xg_mult->SetParName(0,"norm"); //l5xg_mult->SetParLimits(0,-80,80);
534  l5xg_mult->SetParName(1,"mu"); //l5xg_mult->SetParLimits(1,-1.5,1.5);
535  l5xg_mult->SetParName(2,"Sigma"); //l5xg_mult->SetParLimits(2,0.2,0.8);
536  l5xg_mult->SetParName(3,"P");
537  l5xg_mult->SetParName(4,"K"); //l5xg_mult->SetParLimits(4,0.0,0.5);
538  l5xg_mult->SetParName(5,"e"); //l5xg_mult->SetParLimits(5,0.0,0.5);
539  l5xg_mult->SetParName(6,"d"); //l5xg_mult->SetParLimits(6,0.0,0.5);
540  l5xg_mult->SetParName(7,"Total");//l5xg_mult->SetParLimits(7,0.01,2.0);
541  l5xg_mult->SetParName(8,"WidthL");
542  l5xg_mult->SetParName(9,"case");
543  l5xg_mult->SetParName(10,"Mu"); //l5xg_mult->SetParLimits(1,-1.5,1.5);
544  Double_t pars[11], errs[11];
545  memset(pars, 0, sizeof(pars));
546  memset(errs, 0, sizeof(errs));
547  l5xg_func->GetParameters(pars);
548 
549  pars[7] = mult_ev; cout<<"MULT_EV = "<<mult_ev<<endl;
550  pars[9] = 5;
551  pars[10] = pars[1];
552  l5xg_mult->SetParameters(pars);
553 
554  l5xg_mult->SetLineColor(1);
555  // l5xg_mult->SetLineStyle(2);
556  memcpy(errs, l5xg_func->GetParErrors(), 8*sizeof(Double_t));
557  errs[10] = errs[1];
558  l5xg_mult->SetParErrors(errs);
559  hist->GetListOfFunctions()->Add(l5xg_mult);
560  l5xg_mult->SetParent(hist);
561  //-- Fits for particles: 0.Pion, 1.Proton, 2.Kaon, 3.Elektron, 4.Deuteron --//
562  for (Int_t i = 1; i <= 5; i++) {//<=5
563  Int_t j = 0;
564  if(i==5)j = 0;//i==5
565  else j = i;
566  //cout<<"(!!:TEST:!!) i = "<<i<<", j = "<<j<<endl;
567  Sep_func(l5xg,t,norm,mu,sg,fProton,fKaon,fElectron,fDeuteron,total,width,j);
568  TF1 *g0 = new TF1(Form("%s_%s",Peaks[j].Name,hist->GetName()),func_lz5xg_mult,-2.,5.,11);
569  Double_t pars0[11], errs0[11];
570  memset(pars0, 0, sizeof(pars0));
571  memset(errs0, 0, sizeof(errs0));
572  g0->GetParameters(pars0);
573  g0->SetParameter(7,mult_ev);
574  g0->SetParameter(9,j);
575  if(i==1)frac_pi+=frac_pr; //cout<<"(!!:TEST:!!) PROTON : frac_pr = "<<frac_pr<<endl;}
576  else if(i==2)frac_pi+=frac_k;// cout<<"(!!:TEST:!!) KAON : frac_k = "<<frac_k<<endl;}
577  else if(i==3)frac_pi+=frac_el;// cout<<"(!!:TEST:!!) ELEKTRON : frac_el = "<<frac_el<<endl;}
578  else if(i==4)frac_pi+=frac_de;// cout<<"(!!:TEST:!!) DEUTERON : frac_de = "<<frac_de<<endl;}
579  else if(i==5){//5
580  frac_pi=1.-frac_pi;
581  //cout<<"(!!:TEST:!!) PION : frac_pi = "<<frac_pi<<endl;
582  //cout<<"!!!!!!!!!!!!!!!!!!PEAK POSS = "<<func_pi->GetParameter(1)<<", max"<<func_pi->GetMaximumX()<<endl;
583  pars[1] = func_pi->GetMaximumX();
584  l5xg_mult->SetParameter(1,pars[1]);
585  }
586  g0->SetLineColor(j+2);
587  g0->SetParent(hist);
588  hist->GetListOfFunctions()->Add(g0);
589  }
590  Double_t X = l5xg_mult->GetParameter(1);
591  Double_t Y = 0;
592  static TPolyMarker *pm = 0;
593  if (pm) delete pm;
594  pm = new TPolyMarker(1, &X, &Y);
595  hist->GetListOfFunctions()->Add(pm);
596  pm->SetMarkerStyle(23);
597  pm->SetMarkerColor(kRed);
598  pm->SetMarkerSize(1.3);
599  hist->Draw();
600 
601  return l5xg_mult;
602 }
603 //________________________________________________________________________________
604 TF1 *FitRL5(const Char_t *hName = "f1_1", Bool_t outer = kFALSE) {
605  TH1 *hist = (TH1 *) gDirectory->Get(hName);
606  if (! hist) return 0;
607  return FitRL5(hist,outer);
608 }
609 //---------------------------------------------------------------------------//
610 //----------------------------LANDAU Z 5 END --------------------------------//
611 //---------------------------------------------------------------------------//
612 
613 
614 //----- Fit the line to plot width(row) from MC data -----------------------//
615 //----- Fit is used in FitRL5 as w width for Inner and Outer sector -------//
616 void pol_fit(TNtuple *FitP)
617 {
618  //fit for MonteCarlo MIP pion function to find landau width
619  TString command("a4:y>>hist(");
620  command += 47;
621  command += ",0,46)";
622  FitP->Draw(command, "", "");
623  TH1* hist = (TH1*) gDirectory->Get("hist");
624  hist->SetDirectory(0);
625  hist->GetXaxis()->SetTitle("t");
626  hist->GetYaxis()->SetTitle("width");
627 
628  hist->Draw();
629 
630  TF1 *f1 = new TF1("f1","[0]",0.,13.);
631  f1->SetParameter(0,0.076); //Inner: 0.0762
632  f1->SetLineColor(2);
633 
634  TF1 *f2 = new TF1("f2","[0]",13.,50.);
635  f2->SetParameter(0,0.07);//Outer: 0.0699
636  f2->SetLineColor(2);
637 
638  hist->Fit("f1", "R");
639  hist->Fit("f2", "R");
640  f1->Draw("same");
641  f2->Draw("same");
642 }
643 
644 
645 //________________________________________________________________________________
646 TF1 *FitRL1lxg() {
647  static TF1 *flxg = 0;
648  if (flxg) return flxg;
649  static RooFFTConvPdf *lxg = 0;
650  static RooRealVar *t = 0, *ml = 0, *sl = 0, *mg = 0, *sg = 0;
651  static RooLandauZ *landauz = 0;
652  static RooGaussian *gauss = 0;
653  if (! lxg) {
654  // S e t u p c o m p o n e n t p d f s
655  // ---------------------------------------
656  // Construct observable
657  t = new RooRealVar("t","t",-2,6) ;
658  // Construct landauz(t,ml,sl) ;
659  ml = new RooRealVar("ml","log mean landauz",0.0,-20,20) ;
660  sl = new RooRealVar("sl","sigma landauz",0.10,0.01,10) ;
661  landauz = new RooLandauZ("lx","lx",*t,*ml,*sl) ;
662  // Construct gauss(t,mg,sg)
663  mg = new RooRealVar("mg","mg",0) ;
664  sg = new RooRealVar("sg","sg",0.25,0.01,10) ;
665  gauss = new RooGaussian ("gauss","gauss",*t,*mg,*sg) ;
666  // C o n s t r u c t c o n v o l u t i o n p d f
667  // ---------------------------------------
668  // Set #bins to be used for FFT sampling to 10000
669  t->setBins(10000,"cache") ;
670  // Construct landauz (x) gauss
671  lxg = new RooFFTConvPdf("lxg","landauZ (X) gauss",*t,*landauz,*gauss) ;
672  }
673  flxg = (TF1*)lxg->asTF(RooArgList(*t), RooArgList(*ml,*sl,*sg),*t );
674  return flxg;
675 }
676 //________________________________________________________________________________
677 Double_t gfR5Func(Double_t *x, Double_t *par) {
678  // par[0] - norm
679  // par[1] - pion position wrt Z_pion (Bichsel prediction)
680  // par[2] - sigma
681  // par[3] - proton signal
682  // par[4] - Kaon -"-
683  // par[5] - electorn -"-
684  // par[6] - deuteron -"-
685  // par[7] - Total
686  // par[8] - width of Landau
687  Double_t sigma = par[2];
688  Double_t frac[5];
689  Int_t i;
690  frac[0] = 1;
691  for (i = 1; i < 5; i++) {
692  frac[i] = TMath::Sin(par[2+i]);
693  frac[i] *= frac[i];
694  frac[0] -= frac[i];
695  }
696  if (frac[0] < 0.4) return 0;
697  Double_t Value = 0;
698  Int_t i1 = 0;
699  Int_t i2 = 4;
700  Int_t icase = (Int_t ) par[9];
701  if (icase >= 0) {i1 = i2 = icase;}
702  // MP, WidthL, SigmaG
703  Double_t parLI[3] = { par[1], par[8], 0.24};
704  TF1 *lxg = FitRL1lxg();
705  for (i = i1; i <= i2; i++) {
706  Double_t Sigma = TMath::Sqrt(sigma*sigma + Peaks[i].sigma*Peaks[i].sigma);
707  parLI[0] = par[1] + Peaks[i].peak;
708  parLI[1] = par[8];
709  parLI[2] = Sigma;
710  Value += frac[i]*lxg->EvalPar(x,parLI);
711  // cout << "i\t" << i << "\tx = " << x[0] << " frac " << frac[i] << "\t" << Value << endl;
712  }
713  return par[7]*TMath::Exp(par[0])*Value;
714 }
715 //________________________________________________________________________________
716 TF1 *FitR5(TH1 *proj, Option_t *opt="", Int_t nhyps = 5) { // fit by 5 landau convoluted with gauss via RooFit
717  // fit in momentum range p = 0.45 - 0.50 GeV/c
718  if (! proj) return 0;
719  TString Opt(opt);
720  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
721  TF1 *g2 = (TF1*) gROOT->GetFunction("R5");
722  if (! g2) {
723  g2 = new TF1("R5",gfR5Func, -5, 5, 10);
724  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
725  g2->SetParName(1,"mu"); //g2->SetParLimits(1,-1.5,1.5);
726  g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.2,0.8);
727  g2->SetParName(3,"P");
728  g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.5);
729  g2->SetParName(5,"e"); g2->SetParLimits(5,0.0,0.5);
730  g2->SetParName(6,"d"); g2->SetParLimits(6,0.0,0.5);
731  g2->SetParName(7,"Total");
732  g2->SetParName(8,"WidthL"); g2->SetParLimits(8,0.01,2.0);
733  g2->SetParName(9,"case");
734  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
735  }
736 
737  Double_t total = proj->Integral()*proj->GetBinWidth(5);
738  g2->SetParameters(0, proj->GetMean(), proj->GetRMS(), 0.0, 0.0, 0.0, 0.0,0.0,0.5,-1);
739  g2->FixParameter(3,0);
740  g2->FixParameter(4,0);
741  g2->FixParameter(5,0);
742  g2->FixParameter(6,0);
743  g2->FixParameter(7,total);
744  g2->FixParameter(9,-1);
745  Int_t iok = proj->Fit(g2,Opt.Data());
746  if (nhyps == 5) {
747  g2->ReleaseParameter(3); g2->SetParLimits(3,0.0,TMath::Pi()/2);
748  g2->ReleaseParameter(4); g2->SetParLimits(4,0.0,TMath::Pi()/2);
749  g2->ReleaseParameter(5); g2->SetParLimits(5,0.0,TMath::Pi()/2);
750  g2->ReleaseParameter(6); g2->SetParLimits(6,0.0,TMath::Pi()/2);
751  iok = proj->Fit(g2,Opt.Data());
752  }
753  if ( iok < 0) {
754  cout << g2->GetName() << " fit has failed with " << iok << " for "
755  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
756  proj->Fit(g2,Opt.Data());
757  }
758  Opt += "m";
759  iok = proj->Fit(g2,Opt.Data());
760  if (! Opt.Contains("q",TString::kIgnoreCase)) {
761  Double_t params[10];
762  g2->GetParameters(params);
763  Double_t X = params[1];
764  Double_t Y = TMath::Exp(params[0]);
765  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
766  proj->GetListOfFunctions()->Add(pm);
767  pm->SetMarkerStyle(23);
768  pm->SetMarkerColor(kRed);
769  pm->SetMarkerSize(1.3);
770  if (nhyps == 5) {
771  for (int i = 0; i < nhyps; i++) {
772  TF1 *f = new TF1(*g2);
773  f->SetName(Form("L5%s",Peaks[i].Name));
774  f->FixParameter(9,i);
775  f->SetLineColor(i+2);
776  proj->GetListOfFunctions()->Add(f);
777  }
778  proj->Draw();
779  }
780  }
781  return g2;
782 }
783 
784 //________________________________________________________________________________
785 //RooFitResult *FitRL1(TH1 *hist) {
786 TF1 *FitRL1(TH1 *hist) {
787  if (! hist) return 0;
788  // S e t u p c o m p o n e n t p d f s
789  // ---------------------------------------
790  // Construct observable
791  RooRealVar t("t","t",-2,6) ;
792  // Construct landauz(t,ml,sl) ;
793  RooRealVar ml("ml","log mean landauz",0.0,-20,20) ;
794  RooRealVar sl("sl","sigma landauz",0.10,0.01,10) ;
795  RooLandauZ landauz("lx","lx",t,ml,sl) ;
796  // Construct gauss(t,mg,sg)
797  RooRealVar mg("mg","mg",0) ;
798  RooRealVar sg("sg","sg",0.25,0.01,10) ;
799  RooGaussian gauss("gauss","gauss",t,mg,sg) ;
800  // C o n s t r u c t c o n v o l u t i o n p d f
801  // ---------------------------------------
802  // Set #bins to be used for FFT sampling to 10000
803  t.setBins(10000,"cache") ;
804  // Construct landauz (x) gauss
805  RooFFTConvPdf lxg("lxg","landauZ (X) gauss",t,landauz,gauss) ;
806  // S a m p l e , f i t a n d p l o t c o n v o l u t e d p d f
807  RooDataHist* data = new RooDataHist("data","data",t,Import(*hist));
808  // Fit gxlx to data
809  //RooFitResult* r =
810  lxg.fitTo(*data,Save()) ;
811  // Plot data, landauz pdf, landauz (X) gauss pdf
812  RooPlot* frame = t.frame(Title("landauz (x) gauss convolution")) ;
813  data->plotOn(frame) ;
814  lxg.plotOn(frame) ;
815  // landauz.plotOn(frame,LineStyle(kDashed)) ;
816  // Draw frame on canvas
817  TCanvas *ca = (TCanvas *) gROOT->GetListOfCanvases()->FindObject("FitRL1");
818  if (! ca) ca = new TCanvas("FitRL1","FitRL1",600,600) ;
819  else ca->Clear();
820  // gPad->SetLeftMargin(0.15) ; frame->GetYaxis()->SetTitleOffset(1.4) ;
821  frame->Draw() ;
822  delete data;
823  return 0;//r;
824 }
825 //________________________________________________________________________________
826 //RooFitResult *FitRL1(const Char_t *hName = "f1_1") {
827 TF1 *FitRL1(const Char_t *hName = "f1_1") {
828  TH1 *hist = (TH1 *) gDirectory->Get(hName);
829  if (! hist) return 0;
830  return FitRL1(hist);
831 }
832 #endif /* __USE_ROOFIT__ */
833 //________________________________________________________________________________
834 TF1* Landau(){
835  if (!LandauF)
836  LandauF = // P10
837  new TF1("LandauF","exp([0]-0.5*((x-[1])/[2])**2+exp([3]-0.5*((x-[4])/[5])**2+exp([6]-0.5*((x-[7])/[8])**2)))",-5,10);
838  // dEdxP->Draw("phi*sigma_z:(z-zm)/sigma_z>>Lan(100,-3,7)","bg>1&&x>1&&x<3","prof");
839  // Lan->Fit("LandauF","i")
840  Double_t params[9] = {
841  -7.74975e+00,
842  6.53414e+00,
843  1.21524e+00,
844  3.31409e+00,
845  -2.58291e+00,
846  3.51463e+00,
847  -3.47755e+00,
848  3.77698e-02,
849  6.67913e-01};
850  LandauF->SetParameters(params);
851  return LandauF;
852 }
853 //________________________________________________________________________________
854 Double_t fithfcn(Double_t *x,Double_t *par) {
855  Double_t z = x[0];
856  Double_t sigma = par[0];
857  Double_t norm = 1./TMath::Sqrt(2*TMath::Pi())/sigma;
858  Double_t value = 0;
859  // Int_t hyp = (Int_t) par[1+2*NH];
860  // Double_t ref = par[hyp+1];
861  for (int k = 0; k < NH; k++) {
862  Double_t mu = par[k+1];
863  // if (k != hyp) mu -= ref;
864  Double_t dev = (z - mu)/sigma;
865  value += norm*TMath::Exp(par[k+1+NH]-dev*dev/2.);
866  }
867  return value;
868 }
869 //________________________________________________________________________________
870 void FitH(const Char_t *set="z", Int_t Hyp = -1, Int_t Bin=-1) {
871  if (!gBichsel) {
872  gSystem->Load("StBichsel");
873  gBichsel = Bichsel::Instance();
874  }
875  TString Set(set);
876  const Double_t window = 0.4;
877  const Double_t range[2] = {-2., 4.};
878  const Double_t LFrMin = -10;
879  if (! canvas) {
880  canvas = new TCanvas("FitHC","FitH Canvas");
881  canvas->SetGrid();
882  }
883  const Int_t nHYPS = NHYPS;
884  TH2 *hists[nHYPS];
885  TProfile *histp[nHYPS];
886  TFile *fRootFile = (TFile *) gDirectory->GetFile();
887  if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
888  else printf("%s found\n",fRootFile->GetName());
889  TString newfile("FitH");
890  newfile += Set;
891  newfile += gSystem->BaseName(fRootFile->GetName());
892  TString FileN("FitPars");
893  FileN += Set;
894  FileN += ".h";
895  TFile * f = 0;
896  if (Bin < 0) f = new TFile(newfile.Data(),"recreate");
897  TH1 *proj = 0;
898  TF1 *g = new TF1("g",fithfcn,range[0],range[1],2+2*NH);
899  g->SetParName(0,"Sigma"); g->SetParLimits(0,0.06,0.12);
900  g->SetParName(1+2*NH,"Hyp");
901  Int_t Nfit = 0;
902  Int_t nh1 = 0;
903  Int_t nh2 = NHYPS-1; //Int_t nh1 = 3, nh2 = 3;
904  if (Hyp >=0 ) {nh1 = nh2 = Hyp;}
905  Int_t parstatus[NH];
906  for (Int_t hyp = nh1; hyp<=nh2; hyp++) {
907  Int_t kh = hyp%NH;
908  Char_t *HistN = (Char_t *) HistNames[hyp];
909  if (Set.Contains("70")) HistN = (Char_t *) HistNames70[hyp];
910  hists[hyp] = (TH2 *) fRootFile->Get(HistN);
911  if (!hists[hyp]) {printf("Cannot histogram %s\n",HistNames[hyp]); continue;}
912  histp[hyp] = (TProfile *) fRootFile->Get(HistNameP[hyp]);
913  if (!histp[hyp]) {printf("Cannot histogram %s\n",HistNameP[hyp]); continue;;}
914  Int_t k = 0;
915  if (hyp > NH) k = NH;
916  for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
917  Int_t lh = hypl%NH;
918  TString name1("LF.");
919  name1 += Names[hypl];
920  TString name2("");
921  name2 += Names[hypl];
922  g->SetParName(lh+1+NH,name1.Data());
923  g->SetParName(lh+1,name2.Data());
924  }
925  Int_t nx = hists[hyp]->GetNbinsX();
926  Int_t jbin1 = 1;
927  Int_t jbin2 = nx;
928  if (Bin > 0) {jbin1 = jbin2 = Bin;}
929  for (int jbin=jbin1; jbin<=jbin2; jbin++) {
930  TString name(Form("%s_%i_%i",hists[hyp]->GetName(),hyp,jbin));
931  proj = hists[hyp]->ProjectionY(name.Data(),jbin,jbin);
932  Int_t ix1=proj->GetXaxis()->FindBin(-window);
933  Int_t ix2=proj->GetXaxis()->FindBin( window);
934  Double_t dinT = 5*proj->Integral();
935  Double_t dint = proj->Integral(ix1,ix2);
936  if (dint < 100.) {
937  printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hists[hyp]->GetName(),jbin,hyp,dint);
938  delete proj;
939  continue;
940  }
941  for (Int_t lh = 0; lh < NH; lh++) {
942  g->ReleaseParameter(lh+1);
943  g->ReleaseParameter(lh+1+NH);
944  parstatus[lh] = 0;
945  }
946  printf("hist:%s bin %i for hyp %i has %10.0f entries\n",hists[hyp]->GetName(),jbin,hyp,dint);
947  Double_t bg = TMath::Power(10., hists[hyp]->GetBinCenter(jbin));
948  Double_t pmom = Masses[hyp]*bg;
949  if (pmom > 0.1) {parstatus[4] = 1; printf("fix muon\n");}
950  if (kh == 4 && pmom > 0.1) continue;
951  Double_t devZ[NH];
952  Double_t zref = 0;
953  if (TString(set) == "70") zref = 1.e-6*gBichsel->GetI70(TMath::Log10(bg),1.0);
954  else zref = 1.e-6*TMath::Exp(gBichsel->GetMostProbableZ(TMath::Log10(bg),1.0));
955  // printf("zref = %f\n",zref);
956  Int_t iok = 1;
957  for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
958  Int_t lh = hypl%NH;
959  Double_t z = 0;
960  bg = pmom/Masses[hypl];
961  if (Set.Contains("70")) z = 1.e-6*gBichsel->GetI70(TMath::Log10(bg),1.0);
962  else z = 1.e-6*TMath::Exp(gBichsel->GetMostProbableZ(TMath::Log10(bg),1.0));
963  devZ[lh] = TMath::Log(z/zref);
964  }
965  for (int l = 0; l < NH; l++) {
966  if (l == kh) continue;
967  if (devZ[l] < -2.0 || devZ[l] > 4.0) {
968  printf("Fix %s dZ = %f\n", Names[NH*(hyp/NH)+l],devZ[l]);
969  parstatus[l] = 1;
970  continue;
971  }
972  if (TMath::Abs(devZ[l]) < 0.01 ||
973  (hyp%NH == 3 && l == 4 && TMath::Abs(devZ[l]) < 0.04)) {
974  printf("Fix %s dZ = %f\n", Names[NH*(hyp/NH)+l],devZ[l]);
975  parstatus[l] = 1;
976  continue;
977  }
978  }
979  for (Int_t hypl = (hyp/NH)*NH; hypl < (hyp/NH+1)*NH; hypl++) {
980  AGAIN:
981  Int_t lh = hypl%NH;
982  if (parstatus[lh]) continue;
983  Double_t windowN = devZ[lh]-window;
984  Double_t windowP = devZ[lh]+window;
985  for (int m = 0; m < NH; m++) {
986  if (m != lh && ! parstatus[m]) {
987  Double_t dev = 0.5*(devZ[m] - devZ[lh]);
988  if (dev < 0 && windowN < devZ[lh] + dev) windowN = devZ[lh] + dev;
989  if (dev > 0 && windowP > devZ[lh] + dev) windowP = devZ[lh] + dev;
990  }
991  }
992  g->SetParameter(lh+1,devZ[lh]);
993  g->SetParLimits(lh+1,windowN, windowP);
994  g->SetParLimits(lh+1+NH,LFrMin, TMath::Log(dinT));
995  if (hypl == hyp)
996  g->SetParameter(lh+1+NH,TMath::Log(dint));
997  else {
998  ix1=proj->GetXaxis()->FindBin(windowN);
999  ix2=proj->GetXaxis()->FindBin(windowP);
1000  Double_t din= proj->Integral(ix1,ix2);
1001  // printf("%s dZ = %f in = %f\n", Names[hypl], devZ[lh], din);
1002  if (din > 0) g->SetParameter(lh+1+NH,TMath::Log(din));
1003  else {
1004  printf("Fix %s din = %f\n", Names[hypl],din);
1005  parstatus[lh] = 1;
1006  goto AGAIN;
1007  }
1008  }
1009  }
1010  for (int l = 0; l < NH; l++)
1011  if (parstatus[l]) {
1012  printf("Fix %s\n",Names[NH*(hyp/NH)+l]);
1013  g->FixParameter(l+1+NH, LFrMin);
1014  g->FixParameter(l+1, devZ[l]);
1015  }
1016  g->FixParameter(1+2*NH, hyp);
1017  if (parstatus[kh]) iok = 0;
1018  if (! iok) {printf("Too close\n"); continue;}
1019  // proj->Fit(g->GetName(),"RIM");
1020  proj->Fit(g->GetName(),"R");
1021  canvas->Update();
1022  Int_t lh = (hyp%NH)+1;
1023  Double_t Nul = g->GetParameter(lh);
1024  printf("hyp = %i lh = %i Nul = %f zref =%f %f\n",hyp,lh,Nul,zref,g->GetParameter(lh));
1025  Double_t Mul = 1.e6*zref*TMath::Exp(Nul);
1026  Double_t dMul = g->GetParError(lh);
1027  Int_t NFitPoints = g->GetNumberFitPoints();
1028  Int_t NDF = g->GetNDF();
1029  Double_t prob = g->GetProb();//TMath::Prob(chisq, NDF);
1030  Double_t chisq = g->GetChisquare();
1031  Double_t X = hists[hyp]->GetXaxis()->GetBinCenter(jbin);
1032  printf("Nul = %f Mul = %f dMul = %f\n",Nul,Mul,dMul);
1033  printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
1034  "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
1035  Names[hyp],hyp,jbin,Nfit,X,pmom,Nul,Mul,dMul,chisq,NFitPoints,NDF,prob);
1036  printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
1037  Names[hyp],hyp,jbin,Nfit,X,pmom,Nul,Mul,dMul,chisq,NFitPoints,NDF,prob,g->GetName());
1038  if (f) {
1039  FILE *fp = fopen(FileN.Data(),"a");
1040  if (fp) {
1041  if (Nfit == 0) {
1042  TDatime time;
1043  fprintf (fp,"dEdxPoint_t dEdxZ[] = {\n");
1044  fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
1045  fprintf(fp,
1046  "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
1047  }
1048  if (Nfit == 0) fprintf(fp," {");
1049  else fprintf(fp,",{");
1050  fprintf(fp,
1051  "\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
1052  Names[hyp],hyp,jbin,Nfit,X,pmom,Nul,Mul,dMul,chisq,NFitPoints,NDF,prob,g->GetName());
1053  fclose(fp);
1054  Nfit++;
1055  }
1056  proj->Write();
1057  delete proj;
1058  }
1059  }
1060  }
1061  if (f) {
1062  FILE *fp = fopen(FileN.Data(),"a");
1063  if (fp) fprintf(fp,"};\n");
1064  fclose(fp);
1065  delete f;
1066  }
1067 }
1068 
1069 
1070 //________________________________________________________________________________
1071 TH2F *Project(TH3F *hist) {
1072  if (!hist) return 0;
1073  Int_t nx = hist->GetNbinsX();
1074  Int_t ny = hist->GetNbinsY();
1075  TString name(hist->GetName());
1076  name += "_xy";
1077  TH2F *h = new TH2F(name.Data(),hist->GetTitle(),
1078  nx,hist->GetXaxis()->GetXmin(),hist->GetXaxis()->GetXmax(),
1079  ny,hist->GetYaxis()->GetXmin(),hist->GetYaxis()->GetXmax());
1080  Double_t params[5];
1081  Double_t error;
1082  TF1 *g = new TF1("g","gaus(0)+pol1(3)");
1083  for (int i=1;i<=nx;i++){
1084  for (int j=1;j<=ny;j++){
1085  TH1D *proj = hist->ProjectionZ("f_11",i,i,j,j);
1086  Double_t mean = proj->GetMean();
1087  Double_t sum = proj->Integral();
1088  Double_t rms = proj->GetRMS();
1089  params[0] = sum;
1090  params[1] = mean;
1091  params[2] = rms;
1092  params[3] = 0;
1093  params[4] = 0;
1094  g->SetParameters(params);
1095  // g->SetRange(mean-2*rms,mean+2*rms);
1096  g->SetRange(-1.,1.);
1097  error = 0;
1098  if (sum > 100) {
1099  proj->Fit("g","RQ");
1100  g->GetParameters(params);
1101  error = g->GetParError(1);
1102  }
1103  else error = 0;
1104  h->SetCellContent(i,j,params[1]);
1105  h->SetCellError(i,j,error);
1106  printf("i:%i j:%i sum:%f mean:%f rms:%f mu:%f sigma:%f\n",
1107  i,j,sum,mean,rms,params[1],params[2]);
1108  delete proj;
1109  }
1110  }
1111  return h;
1112 }
1113 //________________________________________________________________________________
1114 TH2D *GetRelYZError(TH3 *hist, const Char_t *Name="_yz") {
1115  if (!hist) return 0;
1116  Int_t nx = hist->GetNbinsX();
1117  Int_t ny = hist->GetNbinsY();
1118  Int_t nz = hist->GetNbinsZ();
1119  TString name(hist->GetName());
1120  name += Name;
1121  // TAxis *fXaxis = hist->GetXaxis();
1122  TAxis *fYaxis = hist->GetYaxis();
1123  TAxis *fZaxis = hist->GetZaxis();
1124  TH2D *h = new TH2D(name.Data(),"Relative error in Space charge (%)",
1125  ny,fYaxis->GetXmin(),fYaxis->GetXmax(),
1126  nz,fZaxis->GetXmin(),fZaxis->GetXmax());
1127  h->SetXTitle("Row number");
1128  h->SetYTitle("Z (cm)");
1129  Int_t bin;
1130  Double_t Scale = TMath::Sqrt(24.*3726./1153.);
1131  for (int iybin=0;iybin<ny;iybin++){
1132  for (int izbin=0;izbin<nz;izbin++){
1133  Double_t cont = 0, err = 0;
1134  for (int ixbin=1;ixbin<=nx;ixbin++){
1135  bin = hist->GetBin(ixbin,iybin,izbin);
1136  cont += hist->GetBinContent(bin);
1137  err += hist->GetBinError(bin)*hist->GetBinError(bin);
1138  }
1139  if (cont > 0) {
1140  bin = h->GetBin(iybin,izbin);
1141  Double_t val = 100*TMath::Sqrt(err)/cont*Scale;
1142  h->SetBinContent(bin,val);
1143  }
1144  }
1145  }
1146  return h;
1147 }
1148 //________________________________________________________________________________
1149 TH2F *ProjectX(TH3F *hist, const Char_t *Name="_yz",const Int_t binx1=0,const Int_t binx2=10000) {
1150  if (!hist) return 0;
1151  Int_t nx = hist->GetNbinsX();
1152  if (nx > binx2) nx = binx2;
1153  Int_t ny = hist->GetNbinsY();
1154  Int_t nz = hist->GetNbinsZ();
1155  TString name(hist->GetName());
1156  name += Name;
1157  // TAxis *fXaxis = hist->GetXaxis();
1158  TAxis *fYaxis = hist->GetYaxis();
1159  TAxis *fZaxis = hist->GetZaxis();
1160  TH2F *h = new TH2F(name.Data(),hist->GetTitle(),
1161  ny,fYaxis->GetXmin(),fYaxis->GetXmax(),
1162  nz,fZaxis->GetXmin(),fZaxis->GetXmax());
1163  // Double_t params[3];
1164  // Double_t error;
1165  // TF1 *g = new TF1("g","gaus");
1166  for (int iybin=0;iybin<ny;iybin++){
1167  for (int izbin=0;izbin<nz;izbin++){
1168  for (int ixbin=binx1;ixbin<nx;ixbin++){
1169  Int_t bin = hist->GetBin(ixbin,iybin,izbin);
1170  Double_t cont = hist->GetBinContent(bin);
1171  if (cont) h->Fill(fYaxis->GetBinCenter(iybin),fZaxis->GetBinCenter(izbin), cont);
1172  }
1173  }
1174  }
1175  return h;
1176 }
1177 //________________________________________________________________________________
1178 TF1 *FitGP(TH1 *proj, Option_t *opt="RQ", Double_t nSigma=3, Int_t pow=3) {
1179  if (! proj) return 0;
1180  TString Opt(opt);
1181  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
1182  TF1 *g = 0, *g0 = 0;
1183  TF1 *gaus = (TF1*) gROOT->GetFunction("gaus");
1184  if (pow >= 0) g0 = new TF1("g0",Form("gaus(0)+pol%i(3)",pow),-0.2,0.2);
1185  else g0 = new TF1("g0","gaus",-0.2,0.2);
1186  g0->SetParName(0,"Constant");
1187  g0->SetParName(1,"Mean");
1188  g0->SetParName(2,"Sigma");
1189  for (int i=0; i<=pow;i++) g0->SetParName(3+i,Form("a%i",i));
1190  TF1 *g1 = new TF1("g1",Form("gaus(0)+pol%i(3)",pow+1),-0.2,0.2);
1191  g1->SetParName(0,"Constant");
1192  g1->SetParName(1,"Mean");
1193  g1->SetParName(2,"Sigma");
1194  for (int i=0; i<=pow+1;i++) g1->SetParName(3+i,Form("a%i",i));
1195  TF1 *g2 = new TF1("g2",Form("gaus(0)+pol%i(3)",pow+2),-0.2,0.2);
1196  g2->SetParName(0,"Constant");
1197  g2->SetParName(1,"Mean");
1198  g2->SetParName(2,"Sigma");
1199  for (int i=0; i<=pow+2;i++) g2->SetParName(3+i,Form("a%i",i));
1200  Double_t params[9];
1201  Int_t peak = proj->GetMaximumBin();
1202  Double_t peakX = proj->GetBinCenter(peak);
1203  params[0] = proj->GetBinContent(peak);
1204  if (peakX > 0.5) {
1205  params[1] = 0;
1206  params[2] = 0.2;
1207  }
1208  else {
1209  params[1] = peakX;
1210  params[2] = proj->GetRMS();
1211  if (params[2] > 0.25) params[2] = 0.25;
1212  }
1213  params[3] = 0;
1214  params[4] = 0;
1215  params[5] = 0;
1216  params[6] = 0;
1217  params[7] = 0;
1218  params[8] = 0;
1219  if (gaus) {
1220  g = gaus;
1221  g->SetParameters(params);
1222  g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
1223  proj->Fit(g,opt);
1224  g->GetParameters(params);
1225  if (g->GetProb() > 0.01) return g;
1226  params[2] = TMath::Abs(params[2]);
1227  }
1228  g = g0;
1229  g->SetParameters(params);
1230  g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
1231  proj->Fit(g,opt);
1232  if (g->GetProb() > 0.01) return g;
1233  g->GetParameters(params);
1234  g = g1;
1235  params[2] = TMath::Abs(params[2]);
1236  g->SetParameters(params);
1237  g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
1238  proj->Fit(g,opt);
1239  if (g->GetProb() > 0.01) return g;
1240  g->GetParameters(params);
1241  g = g2;
1242  params[2] = TMath::Abs(params[2]);
1243  g->SetParameters(params);
1244  g->SetRange(params[1]-nSigma*params[2],params[1]+nSigma*params[2]);
1245  proj->Fit(g,opt);
1246  if (! Opt.Contains("q",TString::kIgnoreCase)) {
1247  g->GetParameters(params);
1248  Double_t X = params[1];
1249  Double_t Y = params[0]/TMath::Sqrt(2*TMath::Pi()*params[2]);
1250  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
1251  proj->GetListOfFunctions()->Add(pm);
1252  pm->SetMarkerStyle(23);
1253  pm->SetMarkerColor(kRed);
1254  pm->SetMarkerSize(1.3);
1255  proj->Draw();
1256  }
1257  return g;
1258 }
1259 //________________________________________________________________________________
1260 Double_t gfFunc(Double_t *x, Double_t *par) {
1261  // par[0] - norm
1262  // par[1] - pion position wrt Z_pion (Bichsel prediction)
1263  // par[2] - sigma
1264  // par[3] - proton signal
1265  // par[4] - Kaon -"-
1266  // par[5] - electorn -"-
1267  // par[6] - deuteron -"-
1268  // par[7] - Total
1269  Double_t sigma = par[2];
1270  Double_t frac[5];
1271  Int_t i;
1272  frac[0] = 1;
1273  for (i = 1; i < 5; i++) {
1274  frac[i] = TMath::Sin(par[2+i]);
1275  frac[i] *= frac[i];
1276  frac[0] -= frac[i];
1277  }
1278  if (frac[0] < 0.4) return 0;
1279  Double_t Value = 0;
1280  Int_t icase = (Int_t) par[8];
1281  Int_t i1 = 0;
1282  Int_t i2 = 4;
1283  if (icase >= 0) {i1 = i2 = icase;}
1284  for (i = i1; i <= i2; i++) {
1285  Double_t Sigma = TMath::Sqrt(sigma*sigma + Peaks[i].sigma*Peaks[i].sigma);
1286  Value += frac[i]*TMath::Gaus(x[0],par[1]+Peaks[i].peak,Sigma,1);
1287  // cout << "i\t" << i << "\tx = " << x[0] << " frac " << frac[i] << "\t" << Value << endl;
1288  }
1289  return par[7]*TMath::Exp(par[0])*Value;
1290 }
1291 //________________________________________________________________________________
1292 TF1 *FitGF(TH1 *proj, Option_t *opt="") {
1293  static TSpectrum *fSpectrum = 0;
1294  if (! fSpectrum) {
1295  fSpectrum = new TSpectrum(6);
1296  }
1297  // fit in momentum range p = 0.45 - 0.50 GeV/c
1298  if (! proj) return 0;
1299  TString Opt(opt);
1300  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
1301  TF1 *g2 = (TF1*) gROOT->GetFunction("GF");
1302  if (! g2) {
1303  g2 = new TF1("GF",gfFunc, -5, 5, 9);
1304  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
1305  g2->SetParName(1,"mu"); g2->SetParLimits(1,-2.5,2.5);
1306  g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.2,0.8);
1307  g2->SetParName(3,"P"); g2->SetParLimits(3,0,0.5);
1308  g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.5);
1309  g2->SetParName(5,"e"); g2->SetParLimits(5,0.0,0.5);
1310  g2->SetParName(6,"d"); g2->SetParLimits(6,0.0,0.5);
1311  g2->SetParName(7,"Total");
1312  g2->SetParName(8,"Case");
1313  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
1314  }
1315  // Find pion peak
1316  Int_t nfound = fSpectrum->Search(proj);
1317  if (nfound < 1) return 0;
1318  Int_t npeaks = 0;
1319  Float_t *xpeaks = fSpectrum->GetPositionX();
1320  Float_t xp = 0;
1321  if (nfound > 2) nfound = 2;
1322  Double_t xpi = 9999;
1323  for (Int_t p = 0; p < nfound; p++) {
1324  xp = xpeaks[p];
1325  Int_t bin = proj->GetXaxis()->FindBin(xp);
1326  Double_t yp = proj->GetBinContent(bin);
1327  Double_t ep = proj->GetBinError(bin);
1328  if (yp-5*ep < 0) continue;
1329  if (xp < xpi) xpi = xp;
1330  }
1331  Double_t total = proj->Integral()*proj->GetBinWidth(5);
1332  // g2->SetParameters(0, proj->GetMean(), proj->GetRMS(), 0.1, 0.1, 0.1, 0.1,0.1,-1.);
1333 // Int_t binmax = proj->GetMaximumBin();
1334 // Double_t xmax = proj->GetXaxis()->GetBinCenter(binmax);
1335  g2->SetParameters(0, xpi, 0.35, 0.6, 0.1, 0.1, 0.1,0.1,-1.);
1336  // g2->FixParameter(3,2.86731e-01);
1337  g2->FixParameter(4,1e-6);
1338  g2->FixParameter(5,1e-6);
1339  g2->FixParameter(6,1e-6);
1340  g2->FixParameter(7,total);
1341  g2->FixParameter(8,-1);
1342  proj->Fit(g2,Opt.Data());
1343  g2->ReleaseParameter(3); g2->SetParLimits(3,0.0,TMath::Pi()/2);
1344  g2->ReleaseParameter(4); g2->SetParLimits(4,0.0,TMath::Pi()/2);
1345  g2->ReleaseParameter(5); g2->SetParLimits(5,0.0,TMath::Pi()/2);
1346  g2->ReleaseParameter(6); g2->SetParLimits(6,0.0,TMath::Pi()/2);
1347  Int_t iok = proj->Fit(g2,Opt.Data());
1348  if ( iok < 0) {
1349  cout << g2->GetName() << " fit has failed with " << iok << " for "
1350  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
1351  proj->Fit(g2,Opt.Data());
1352  }
1353  Opt += "m";
1354  iok = proj->Fit(g2,Opt.Data());
1355  if (iok < 0 ) return 0;
1356  if (! Opt.Contains("q",TString::kIgnoreCase)) {
1357  Double_t params[10];
1358  g2->GetParameters(params);
1359  Double_t X = params[1];
1360  Double_t Y = TMath::Exp(params[0]);
1361  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
1362  proj->GetListOfFunctions()->Add(pm);
1363  pm->SetMarkerStyle(23);
1364  pm->SetMarkerColor(kRed);
1365  pm->SetMarkerSize(1.3);
1366  for (int i = 0; i <= 4; i++) {
1367  TF1 *f = new TF1(*g2);
1368  f->SetName(Peaks[i].Name);
1369  f->FixParameter(8,i);
1370  f->SetLineColor(i+2);
1371  proj->GetListOfFunctions()->Add(f);
1372  }
1373  proj->Draw();
1374  }
1375  return g2;
1376 }
1377 //________________________________________________________________________________
1378 Double_t langaufun(Double_t *x, Double_t *par) {// $ROOTSYS/tutorials/fit/langaus.C
1379 
1380  //Fit parameters:
1381  //par[0]=Width (scale) parameter of Landau density
1382  //par[1]=Most Probable (MP, location) parameter of Landau density
1383  //par[2]=Total area (integral -inf to inf, normalization constant)
1384  //par[3]=Width (sigma) of convoluted Gaussian function
1385  //
1386  //In the Landau distribution (represented by the CERNLIB approximation),
1387  //the maximum is located at x=-0.22278298 with the location parameter=0.
1388  //This shift is corrected within this function, so that the actual
1389  //maximum is identical to the MP parameter.
1390 
1391  // Numeric constants
1392  Double_t invsq2pi = 0.3989422804014; // (2 pi)^(-1/2)
1393  Double_t mpshift = -0.22278298; // Landau maximum location
1394 
1395  // Control constants
1396  Double_t np = 100.0; // number of convolution steps
1397  Double_t sc = 5.0; // convolution extends to +-sc Gaussian sigmas
1398 
1399  // Variables
1400  Double_t xx, zz;
1401  Double_t mpc;
1402  Double_t fland;
1403  Double_t sum = 0.0;
1404  Double_t xlow,xupp;
1405  Double_t step;
1406  Double_t i;
1407 
1408 
1409  // MP shift correction
1410  mpc = TMath::Exp(par[1]) - mpshift * par[0];
1411 
1412  // Range of convolution integral
1413  xlow = x[0] - sc * par[3];
1414  xupp = x[0] + sc * par[3];
1415 
1416  step = (xupp-xlow) / np;
1417  // Convolution integral of Landau and Gaussian by sum
1418  for(i=1.0; i<=np/2; i++) {
1419  zz = xlow + (i-.5) * step;
1420  xx = TMath::Exp(zz);
1421  fland = xx*TMath::Landau(xx,mpc,par[0]) / par[0];
1422  sum += fland * TMath::Gaus(x[0],zz,par[3]);
1423 
1424  zz = xupp - (i-.5) * step;
1425  xx = TMath::Exp(zz);
1426  fland = xx*TMath::Landau(xx,mpc,par[0]) / par[0];
1427  sum += fland * TMath::Gaus(x[0],zz,par[3]);
1428  }
1429 
1430  return (par[2] * step * sum * invsq2pi / par[3]);
1431 }
1432 
1433 //________________________________________________________________________________
1434 Double_t gfL5Func(Double_t *x, Double_t *par) {
1435  // par[0] - norm
1436  // par[1] - pion position wrt Z_pion (Bichsel prediction)
1437  // par[2] - sigma
1438  // par[3] - proton signal
1439  // par[4] - Kaon -"-
1440  // par[5] - electorn -"-
1441  // par[6] - deuteron -"-
1442  // par[7] - Total
1443  // par[8] - width of Landau
1444  Double_t sigma = par[2];
1445  Double_t frac[5];
1446  Int_t i;
1447  frac[0] = 1;
1448  for (i = 1; i < 5; i++) {
1449  frac[i] = TMath::Sin(par[2+i]);
1450  frac[i] *= frac[i];
1451  frac[0] -= frac[i];
1452  }
1453  if (frac[0] < 0.4) return 0;
1454  Double_t Value = 0;
1455  Int_t i1 = 0;
1456  Int_t i2 = 4;
1457  Int_t icase = (Int_t ) par[9];
1458  if (icase >= 0) {i1 = i2 = icase;}
1459  // WidthL MP Total, SigmaG
1460  Double_t parLI[4] = { par[8], par[1], 1, 0.24};
1461  for (i = i1; i <= i2; i++) {
1462  Double_t Sigma = TMath::Sqrt(sigma*sigma + Peaks[i].sigma*Peaks[i].sigma);
1463  parLI[1] = par[1] + Peaks[i].peak;
1464  parLI[3] = Sigma;
1465  Value += frac[i]*langaufun(x,parLI);
1466  // cout << "i\t" << i << "\tx = " << x[0] << " frac " << frac[i] << "\t" << Value << endl;
1467  }
1468  return par[7]*TMath::Exp(par[0])*Value;
1469 }
1470 //________________________________________________________________________________
1471 TF1 *FitL5(TH1 *proj, Option_t *opt="", Int_t nhyps = 5) { // fit by 5 landau convoluted with gauss
1472  // fit in momentum range p = 0.45 - 0.50 GeV/c
1473  if (! proj) return 0;
1474  TString Opt(opt);
1475  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
1476  TF1 *g2 = (TF1*) gROOT->GetFunction("L5");
1477  if (! g2) {
1478  g2 = new TF1("L5",gfL5Func, -5, 5, 10);
1479  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
1480  g2->SetParName(1,"mu"); //g2->SetParLimits(1,-1.5,1.5);
1481  g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.2,0.8);
1482  g2->SetParName(3,"P");
1483  g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.5);
1484  g2->SetParName(5,"e"); g2->SetParLimits(5,0.0,0.5);
1485  g2->SetParName(6,"d"); g2->SetParLimits(6,0.0,0.5);
1486  g2->SetParName(7,"Total");
1487  g2->SetParName(8,"WidthL"); g2->SetParLimits(8,0.01,2.0);
1488  g2->SetParName(9,"case");
1489  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
1490  }
1491 
1492  Double_t total = proj->Integral()*proj->GetBinWidth(5);
1493  g2->SetParameters(0, proj->GetMean(), proj->GetRMS(), 0.0, 0.0, 0.0, 0.0,0.0,0.5,-1);
1494  g2->FixParameter(3,0);
1495  g2->FixParameter(4,0);
1496  g2->FixParameter(5,0);
1497  g2->FixParameter(6,0);
1498  g2->FixParameter(7,total);
1499  g2->FixParameter(9,-1);
1500  Int_t iok = proj->Fit(g2,Opt.Data());
1501  if (nhyps == 5) {
1502  g2->ReleaseParameter(3); g2->SetParLimits(3,0.0,TMath::Pi()/2);
1503  g2->ReleaseParameter(4); g2->SetParLimits(4,0.0,TMath::Pi()/2);
1504  g2->ReleaseParameter(5); g2->SetParLimits(5,0.0,TMath::Pi()/2);
1505  g2->ReleaseParameter(6); g2->SetParLimits(6,0.0,TMath::Pi()/2);
1506  iok = proj->Fit(g2,Opt.Data());
1507  }
1508  if ( iok < 0) {
1509  cout << g2->GetName() << " fit has failed with " << iok << " for "
1510  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
1511  proj->Fit(g2,Opt.Data());
1512  }
1513  Opt += "m";
1514  iok = proj->Fit(g2,Opt.Data());
1515  if (! Opt.Contains("q",TString::kIgnoreCase)) {
1516  Double_t params[10];
1517  g2->GetParameters(params);
1518  Double_t X = params[1];
1519  Double_t Y = TMath::Exp(params[0]);
1520  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
1521  proj->GetListOfFunctions()->Add(pm);
1522  pm->SetMarkerStyle(23);
1523  pm->SetMarkerColor(kRed);
1524  pm->SetMarkerSize(1.3);
1525  if (nhyps == 5) {
1526  for (int i = 0; i < nhyps; i++) {
1527  TF1 *f = new TF1(*g2);
1528  f->SetName(Form("L5%s",Peaks[i].Name));
1529  f->FixParameter(9,i);
1530  f->SetLineColor(i+2);
1531  proj->GetListOfFunctions()->Add(f);
1532  }
1533  proj->Draw();
1534  }
1535  }
1536  return g2;
1537 }
1538 //________________________________________________________________________________
1539 Double_t gbFunc(Double_t *x, Double_t *par) {
1540  // par[0] - norm
1541  // par[1] - pion position wrt Z_pion (Bichsel prediction)
1542  // par[2] - sigma
1543  // par[3] - proton signal
1544  // par[4] - Kaon -"-
1545  // par[5] - electorn -"-
1546  // par[6] - deuteron -"-
1547  // par[7] - Total
1548  // par[8] - <dX>
1549  if (! LandauF) Landau();
1550  static Double_t sigma_p[3] = {// sigma versus ::log(dX)
1551  5.31393e-01,// 1 p0 1.33485e-03 7.13072e-07 7.08416e-08
1552  -1.43277e-01,// 2 p1 3.36846e-03 6.62434e-07 -1.13681e-05
1553  2.43800e-02 // 3 p2 1.81240e-03 4.02492e-07 -2.08423e-05
1554  };
1555  Double_t sigmaC = par[2];
1556  Double_t frac[5];
1557  Int_t i;
1558  frac[0] = 1;
1559  for (i = 1; i < 5; i++) {
1560  frac[i] = TMath::Sin(par[2+i]);
1561  frac[i] *= frac[i];
1562  frac[0] -= frac[i];
1563  }
1564  Double_t Value = 0;
1565  static Double_t pMom = 0.475;
1566  static Double_t Xlog10bg[5];
1567  Double_t Ylog2dx = TMath::Log2(par[8]);
1568  Double_t Sigma = sigma_p[2];
1569  for (int n = 1; n>=0; n--) Sigma = Ylog2dx*Sigma + sigma_p[n];
1570  Double_t zMostProb[5];
1571  for (i = 0; i < 5; i++) {
1572  Xlog10bg[i] = TMath::Log10(pMom/Peaks[i].mass);
1573  zMostProb[i] = gBichsel->GetMostProbableZ(Xlog10bg[i],Ylog2dx);
1574  Double_t sigma = Sigma + sigmaC;
1575  // Double_t xi = (x[0] + zMostProb[0] - zMostProb[i])/sigma;
1576  // Double_t xi = (x[0] + par[1] + zMostProb[0] - zMostProb[i])/sigma;
1577  Double_t xi = (x[0] - par[1] - Peaks[i].peak)/sigma;
1578  Double_t Prob = LandauF->Eval(xi);
1579 
1580  Value += frac[i]*Prob;
1581  // cout << "i\t" << i << "\tx = " << x[0] << " frac " << frac[i] << "\t" << Value << endl;
1582  }
1583  return par[7]*TMath::Exp(par[0])*Value;
1584 }
1585 //________________________________________________________________________________
1586 TF1 *FitGB(TH1 *proj, Option_t *opt="", Double_t dX = 2.364) {
1587  if (!gBichsel) {
1588  gSystem->Load("StBichsel");
1589  gBichsel = Bichsel::Instance();
1590  }
1591  // fit in momentum range p = 0.45 - 0.50 GeV/c
1592  if (! proj) return 0;
1593  TString Opt(opt);
1594  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
1595  TF1 *g2 = (TF1*) gROOT->GetFunction("GB");
1596  if (! g2) {
1597  g2 = new TF1("GB",gbFunc, -5, 5, 9);
1598  g2->SetParName(0,"norm");
1599  g2->SetParName(1,"mu"); g2->SetParLimits(1,-1.5,1.5);
1600  g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.0,0.8);
1601  g2->SetParName(3,"P");
1602  g2->SetParName(4,"K");
1603  g2->SetParName(5,"e");
1604  g2->SetParName(6,"d");
1605  g2->SetParName(7,"Total");
1606  g2->SetParName(8,"dX");
1607  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
1608  }
1609 
1610  Double_t total = proj->Integral()*proj->GetBinWidth(5);
1611  g2->SetParameters(0, 1e-3, 0.01, 0.4, 0., 0., 0.,0.);
1612  g2->FixParameter(7,total);
1613  g2->FixParameter(8,dX);
1614  Int_t iok = proj->Fit(g2,Opt.Data());
1615  if ( iok < 0 ) {
1616  cout << g2->GetName() << " fit has failed with " << iok << " for "
1617  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
1618  proj->Fit(g2,Opt.Data());
1619  }
1620  Opt += "m";
1621  iok = proj->Fit(g2,Opt.Data());
1622  if (! Opt.Contains("q",TString::kIgnoreCase)) {
1623  Double_t params[10];
1624  g2->GetParameters(params);
1625  Double_t X = params[1];
1626  Double_t Y = TMath::Exp(params[0]);
1627  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
1628  proj->GetListOfFunctions()->Add(pm);
1629  pm->SetMarkerStyle(23);
1630  pm->SetMarkerColor(kRed);
1631  pm->SetMarkerSize(1.3);
1632  proj->Draw();
1633  }
1634  return g2;
1635 }
1636 //________________________________________________________________________________
1637 TF1 *FitG2(TH1 *proj, Option_t *opt="RQ") {
1638  if (! proj) return 0;
1639  Double_t params[9];
1640  TF1 *gaus = new TF1("gaus","gaus",-5.,5.);
1641  proj->Fit(gaus,opt);
1642  params[0] = gaus->GetParameter(0);
1643  params[1] = gaus->GetParameter(1);
1644  params[2] = gaus->GetParameter(2);
1645  params[3] = 1;
1646  params[4] = 2;
1647  params[5] = gaus->GetParameter(2);
1648  params[6] = 0;
1649  params[7] = 0;
1650  params[8] = 0;
1651  TF1 *g = new TF1("g","gaus(0)+gaus(3)",-5.,5.);
1652  g->SetParameters(params);
1653  g->SetParLimits(0,0,1.e10);
1654  g->SetParLimits(1,-5,5);
1655  g->SetParLimits(2,0.01,5);
1656  g->SetParLimits(3,0,1.e10);
1657  g->SetParLimits(4,-5,5);
1658  g->SetParLimits(5,0.01,5);
1659  proj->Fit(g,opt);
1660  g->GetParameters(params);
1661  if (TMath::Abs(params[1]) > TMath::Abs(params[4])) {
1662  g->SetParameter(0,params[3]);
1663  g->SetParameter(1,params[4]);
1664  g->SetParameter(2,params[5]);
1665  g->SetParameter(3,params[0]);
1666  g->SetParameter(4,params[1]);
1667  g->SetParameter(5,params[2]);
1668  proj->Fit(g,opt);
1669  }
1670  proj->Fit(g,opt);
1671  delete gaus;
1672  return g;
1673 }
1674 //________________________________________________________________________________
1675 TF1 *FitG3(TH1 *proj, Option_t *opt="RQ") {
1676  if (! proj) return 0;
1677  Double_t params[9];
1678  TF1 *gaus = new TF1("gaus","gaus",-5.,5.);
1679  TF1 *g = 0;
1680  proj->Fit(gaus,opt);
1681  params[0] = gaus->GetParameter(0);
1682  params[1] = gaus->GetParameter(1);
1683  params[2] = gaus->GetParameter(2);
1684  params[3] = 1;
1685  params[4] = 2;
1686  params[5] = gaus->GetParameter(2);
1687  params[6] = 0;
1688  params[7] = 0;
1689  params[8] = 0;
1690  delete gaus;
1691  g = new TF1("g","gaus(0)+gaus(3)",-5.,5.);
1692  g->SetParameters(params);
1693  g->SetParLimits(0,0,1.e10);
1694  g->SetParLimits(1,-5,5);
1695  g->SetParLimits(2,0.1,5);
1696  g->SetParLimits(3,0,1.e10);
1697  g->SetParLimits(4,-5,5);
1698  g->SetParLimits(5,0.1,5);
1699  proj->Fit(g,opt);
1700  g->GetParameters(params);
1701  if (TMath::Abs(params[1]) > TMath::Abs(params[4])) {
1702  g->SetParameter(0,params[3]);
1703  g->SetParameter(1,params[4]);
1704  g->SetParameter(2,params[5]);
1705  g->SetParameter(3,params[0]);
1706  g->SetParameter(4,params[1]);
1707  g->SetParameter(5,params[2]);
1708  }
1709  proj->Fit(g,opt);
1710  if (g->GetProb() > 0) return g;
1711  g->GetParameters(params);
1712  delete g;
1713  g = new TF1("g","gaus(0)+gaus(3)+gaus(6)",-5.,5.);
1714  g->SetParLimits(6,0,1.e10);
1715  g->SetParLimits(7,-5,5);
1716  g->SetParLimits(8,0.1,5);
1717  params[6] = 10.;
1718  params[7] = 0.5*(params[1]+params[4]);
1719  params[8] = params[2];
1720  g->SetParameters(params);
1721  proj->Fit(g,opt);
1722  return g;
1723 }
1724 //________________________________________________________________________________
1725 void FitB4G(Int_t icase = 0, Int_t hyp=-1, Int_t bin=0,
1726  Double_t xmin1=-1.0, Double_t xmax1 = 1.0,
1727  Double_t Mu2 = -.5,Double_t xmin2=-2., Double_t xmax2 = 4.,
1728  Double_t Mu3 = 0.1,Double_t xmin3=-2., Double_t xmax3 = 4.,
1729  Double_t Mu4 = 0.3
1730  )
1731 { // icase = 0 z; icase != 0 70
1732  Double_t sigmas[2] = {0.06,0.12};
1733  const Int_t nHYPS = NHYPS;
1734  TH2 *hists[nHYPS];
1735  TProfile *histp[nHYPS];
1736  TFile *fRootFile = (TFile *) gDirectory->GetFile();
1737  if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
1738  TFile *f = 0;
1739  if (bin == 0) {
1740  TString newfile("FitB4G");
1741  newfile += gSystem->BaseName(fRootFile->GetName());
1742  f = new TFile(newfile.Data(),"update");
1743  }
1744  Int_t i;
1745  for (i = 0; i< nHYPS; i++) {
1746  if (hyp >= 0 && hyp != i) continue;
1747  if (! icase) {
1748  hists[i] = (TH2 *) fRootFile->Get(HistNames[i]);
1749  if (!hists[i]) {printf("Cannot histogram %s\n",HistNames[i]); return;}
1750  }
1751  else {
1752  hists[i] = (TH2 *) fRootFile->Get(HistNames70[i]);
1753  if (!hists[i]) {printf("Cannot histogram %s\n",HistNames70[i]); return;}
1754  }
1755  histp[i] = (TProfile *) fRootFile->Get(HistNameP[i]);
1756  if (!histp[i]) {printf("Cannot histogram %s\n",HistNameP[i]); return;}
1757  }
1758  TF1 *g1 = new TF1("g1","gaus",xmin1,xmax1);
1759  TF1 *g = 0, *g2 = 0, *g3 = 0, *g4 = 0, *ga = 0;
1760  TCanvas *canvas = new TCanvas("canvas","canvas");
1761  if (Mu2 != 0) {
1762  if (Mu2 < -1.) {
1763  if (Mu2 <= -4.) ga = new TF1("ga","gaus(0)+exp(pol3(3))",xmin2,xmax2);
1764  else {
1765  if (Mu2 <= -3.) ga = new TF1("ga","gaus(0)+exp(pol2(3))",xmin2,xmax2);
1766  else ga = new TF1("ga","gaus(0)+exp(pol1(3))",xmin2,xmax2);
1767  }
1768  g = ga;
1769  g->SetParName(0,"Constant1");
1770  g->SetParName(1,"Mean1");
1771  g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
1772  g->SetParName(3,"Const");
1773  g->SetParName(4,"Slope1");
1774  g->SetParName(5,"Slope2");
1775  g->SetParName(6,"Slope3");
1776  g->SetParName(7,"Slope4");
1777  }
1778  else {
1779  g2 = new TF1("g2","gaus(0)+gaus(3)",xmin2,xmax2);
1780  g = g2;
1781  g->SetParName(0,"Constant1");
1782  g->SetParName(1,"Mean1");
1783  g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
1784  g->SetParName(3,"Constant2");
1785  g->SetParName(4,"Mean2");
1786  g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
1787  if (Mu3 != 0) {
1788  g3= new TF1("g3","gaus(0)+gaus(3)+gaus(6)",xmin3,xmax3);
1789  g = g3;
1790  g->SetParName(0,"Constant1");
1791  g->SetParName(1,"Mean1");
1792  g->SetParName(2,"Sigma1"); g->SetParLimits(2,sigmas[0],sigmas[1]);
1793  g->SetParName(3,"Constant2");
1794  g->SetParName(4,"Mean2");
1795  g->SetParName(5,"Sigma2"); g->SetParLimits(5,sigmas[0],sigmas[1]);
1796  g->SetParName(6,"Constant3");
1797  g->SetParName(7,"Mean3");
1798  g->SetParName(8,"Sigma3"); g->SetParLimits(8,sigmas[0],sigmas[1]);
1799  if (Mu4 != 0) {
1800  g4= new TF1("g4","gaus(0)+gaus(3)+gaus(6)+gaus(9)",xmin3,xmax3);
1801  g4->SetParName(0,"Constant1");
1802  g4->SetParName(1,"Mean1");
1803  g4->SetParName(2,"Sigma1"); g4->SetParLimits(2,sigmas[0],sigmas[1]);
1804  g4->SetParName(3,"Constant2");
1805  g4->SetParName(4,"Mean2");
1806  g4->SetParName(5,"Sigma2"); g4->SetParLimits(5,sigmas[0],sigmas[1]);
1807  g4->SetParName(6,"Constant3");
1808  g4->SetParName(7,"Mean3");
1809  g4->SetParName(8,"Sigma3"); g4->SetParLimits(8,sigmas[0],sigmas[1]);
1810  g4->SetParName(9,"Constant4");
1811  g4->SetParName(10,"Mean4");
1812  g4->SetParName(11,"Sigma4"); g4->SetParLimits(11,sigmas[0],sigmas[1]);
1813  }
1814  }
1815  }
1816  }
1817  TH1 *proj = 0;
1818  Double_t params[12];
1819  Int_t k;//,ibin;
1820  Int_t Bin = TMath::Abs(bin);
1821  for (k = 0; k<nHYPS; k++) {
1822  if (hyp != -1 && k != hyp) continue;
1823  TH2 *hist = hists[k];
1824  Int_t nx = hist->GetNbinsX();
1825  for (i=1; i<=nx; i++) {
1826  if (bin != 0 && Bin != i) continue;
1827  // if (proj) delete proj;
1828  Char_t line[40];
1829  sprintf(line,"%s_%i",hist->GetName(),i);
1830  TString name(line);
1831  proj = hist->ProjectionY(name.Data(),i,i);
1832  Int_t ix1=proj->GetXaxis()->FindBin(-.1);
1833  Int_t ix2=proj->GetXaxis()->FindBin(0.1);
1834  if (proj->Integral(ix1,ix2) < 100.) {
1835  printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hist->GetName(),Bin,k,proj->Integral());
1836  delete proj;
1837  continue;
1838  }
1839  Int_t NFitPoints = 0;
1840  Double_t chisq;
1841  Double_t xrange = 0.6;
1842  g = g1;
1843  g->SetParLimits(0,0,1.e7);
1844  g->SetParLimits(1,-xrange,xrange);
1845  g->SetParLimits(2,sigmas[0],sigmas[1]);
1846  proj->Fit(g->GetName(),"r");
1847  g->GetParameters(params);
1848  Int_t kcase = -1;
1849  if (xmin1 < -0.5 && xmax1 > 0.5 && g->GetChisquare() < 1.e3 // g->GetProb() > 1.e-3
1850  ) goto Done;
1851  // && TMath::Abs(g->GetParameter(1)) < 0.05) goto Done;
1852  if (ga) {
1853  g = ga;
1854  g1->GetParameters(params);
1855  params[1] = 0.;
1856  params[3] = 0.;
1857  params[4] = 0.;
1858  params[5] = 0.;
1859  params[6] = 0.;
1860  params[7] = 0.;
1861  params[8] = 0.;
1862  g->SetParameters(params);
1863  kcase = 0;
1864  if (g->GetChisquare() < 1.e3 // g->GetProb() > 1.e-3
1865  ) goto Done;
1866  }
1867  else {
1868  if (g2) {
1869  params[1] = 0.;
1870  params[3] = params[0];
1871  params[4] = Mu2;
1872  params[5] = 2.0*params[2];
1873  g = g2;
1874  g->SetParameters(params);
1875  g->SetParLimits(0,0,1.e7);
1876  g->SetParLimits(1,-xrange,xrange);
1877  // g->SetParLimits(1,-.3,.3);
1878  g->SetParLimits(2,sigmas[0],sigmas[1]);
1879  g->SetParLimits(3,0,1.e7);
1880  g->SetParLimits(5,sigmas[0],sigmas[1]);
1881  proj->Fit(g->GetName(),"r");
1882  g->GetParameters(params);
1883  kcase = 2;
1884  if (g->GetChisquare() < 1.e3 // g->GetProb() > 1.e-3
1885  ) goto Done;
1886  if (g3) {
1887  params[1] = 0.;
1888  params[6] = params[0];
1889  params[7] = Mu3;
1890  params[8] = 2.0*params[2];
1891  g = g3;
1892  g->SetParameters(params);
1893  g->SetParLimits(0,0,1.e7);
1894  g->SetParLimits(1,-xrange,xrange);
1895  // g->SetParLimits(1,-.1,.1);
1896  g->SetParLimits(2,sigmas[0],sigmas[1]);
1897  g->SetParLimits(3,0,1.e7);
1898  g->SetParLimits(5,sigmas[0],sigmas[1]);
1899  g->SetParLimits(6,0,1.e7);
1900  g->SetParLimits(8,sigmas[0],sigmas[1]);
1901  proj->Fit(g->GetName(),"r");
1902  g->GetParameters(params);
1903  kcase = 3;
1904  if (g->GetChisquare() < 1.e3 // g->GetProb() > 1.e-3
1905  ) goto Done;
1906  if (g4) {
1907  params[1] = 0.;
1908  params[9] = params[0];
1909  params[10] = Mu4;
1910  params[11] = 2.0*params[2];
1911  g = g4;
1912  g->SetParameters(params);
1913  g->SetParLimits(0,0,1.e7);
1914  g->SetParLimits(1,-xrange,xrange);
1915  // g->SetParLimits(1,-.1,.1);
1916  g->SetParLimits(2,sigmas[0],sigmas[1]);
1917  g->SetParLimits(3,0,1.e7);
1918  g->SetParLimits(5,sigmas[0],sigmas[1]);
1919  g->SetParLimits(6,0,1.e7);
1920  g->SetParLimits(8,sigmas[0],sigmas[1]);
1921  g->SetParLimits(9,0,1.e7);
1922  g->SetParLimits(11,sigmas[0],sigmas[1]);
1923  proj->Fit(g->GetName(),"r");
1924  kcase = 4;
1925  g->GetParameters(params);
1926  }
1927  }
1928  }
1929  }
1930  Done:
1931  if (g) {
1932  proj->Fit(g->GetName(),"RM");
1933  canvas->Update();
1934  Int_t l = 1;
1935  Double_t mu = g->GetParameter(l);
1936  for (int m=2;m<=kcase;m++) {
1937  if (TMath::Abs(mu) > TMath::Abs(g->GetParameter(3*m-2))) {
1938  l = 3*m-2; mu = g->GetParameter(l);
1939  // printf("l=%i\n",l);
1940  }
1941  }
1942  Nu[N] = g->GetParameter(l);// printf("l=%i Nu=%f\n",l,Nu[N]);
1943  Mu[N] = Nu[N] + histp[k]->GetBinContent(i);
1944  dMu[N] = g->GetParError(l);//Mu[N];
1945  // Mu[N] = TMath::Log(Mu[N]);
1946  Sigma[N] = g->GetParameter(2);
1947  dSigma[N] = g->GetParError(2);
1948  NFitPoints = g->GetNumberFitPoints();
1949  Int_t NDF = g->GetNDF();
1950  Double_t prob = g->GetProb();//TMath::Prob(chisq, NDF);
1951  chisq = g->GetChisquare();
1952  X[N] = hist->GetXaxis()->GetBinCenter(i);
1953  dX[N] = hist->GetXaxis()->GetBinWidth(i);
1954  Double_t pionM = Masses[k]*pow(10.,X[N]);
1955  printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
1956  "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
1957  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob);
1958  // if (bin < 0 && prob > 1.e-7 && TMath::Abs(Nu[N]) < 0.05) {
1959  // if ( bin >= 0 && prob > 1.e-7 && TMath::Abs(Nu[N]) < 0.10 || bin < 0) {
1960  // if ( bin >= 0 && TMath::Abs(Nu[N]) < 0.10 || bin < 0) {
1961  if ( bin >= 0) {
1962  printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
1963  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
1964  TString FileN("FitPars");
1965  if (hyp > -1) FileN += hist->GetName();
1966  FileN += ".h";
1967  FILE *fp = fopen(FileN.Data(),"a");
1968  if (fp) {
1969  if (N == 0) {
1970  TDatime time;
1971  fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
1972  fprintf(fp,
1973  "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
1974  }
1975  fprintf(fp,
1976  "{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
1977  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
1978  fclose(fp);
1979  }
1980  proj->Write();
1981  }
1982  else printf ("================== Skip it\n");
1983  }
1984  N++;
1985  proj->Draw();// cnew->Update();
1986  }
1987  }
1988  if (f) {delete f;}
1989 }
1990 //________________________________________________________________________________
1991 Int_t FitG(TH1 *proj, TF1 *g, TF1 *ga){//, Double_t scaleM=-2., Double_t scaleP=2.) {
1992  Double_t params[10];
1993  proj->Fit("g","R");
1994  Double_t chisq = g->GetChisquare();
1995  if (chisq <= 0. || chisq > 1.e10) return -1;
1996  g->GetParameters(params);
1997  params[3] = 0;
1998  params[4] = 0;
1999  params[5] = 0;
2000  params[6] = 0;
2001  params[7] = 0;
2002  params[8] = 0;
2003  params[9] = 0;
2004  ga->SetParameters(params);
2005  proj->Fit("ga","r");
2006  // proj->Fit("ga","rIM");
2007  chisq = g->GetChisquare();
2008  if (chisq <= 0. || chisq > 1.e10) return -1;
2009  return 0;
2010 }
2011 //________________________________________________________________________________
2012 Int_t FitGG(TH1 *proj, TF1 *g1, TF1 *g2=0, TF1 *ga2=0, Double_t scaleM=-2., Double_t scaleP=2.) {
2013  Double_t params[9];
2014  proj->Fit("g1","R");
2015  g1->GetParameters(params);
2016  Double_t chisq = g1->GetChisquare();
2017  if (chisq <= 0. || chisq > 1.e10) return -1;
2018  params[3] = 1;
2019  params[4] = params[1]+0.1;
2020  params[5] = params[2];
2021  g2->SetParameters(params);
2022  proj->Fit("g2");
2023  chisq = g2->GetChisquare();
2024  if (chisq <= 0. || chisq > 1.e10) return -1;
2025  g2->GetParameters(params);
2026  if (ga2) {
2027  params[6] = 0;
2028  params[7] = 0;
2029  params[8] = 0;
2030  ga2->SetParameters(params);
2031  proj->Fit("ga2","R");
2032  ga2->GetParameters(params);
2033  ga2->SetRange(params[1]+scaleM*params[2],params[1]+scaleP*params[2]);
2034  proj->Fit("ga2","R");
2035  }
2036 // if (gPad) {
2037 // gPad->Update();
2038 // gPad->WaitPrimitive();
2039 // }
2040  return 0;
2041 }
2042 //________________________________________________________________________________
2043 void FitX(TH2 *hist=0, Double_t range=1, Int_t Ibin = 0) {
2044  TFile *fRootFile = (TFile *) gDirectory->GetFile();
2045  if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
2046  if (!hist) {
2047  const Char_t *HistName = "FPoints";
2048  hist = (TH2D *) fRootFile->Get(HistName);
2049  if (!hist) {printf("Cannot histogram %s\n",HistName); return;}
2050  }
2051  // const Int_t nx = hist->GetNbinsX();
2052  TF1 *g = new TF1("g","gaus",-range,range);
2053  // TF1 *ga = g;
2054  TString fitN("gaus(0)+exp(pol3(3))");
2055  if (range > 2.) fitN = "gaus(0)+exp(pol5(3))";
2056  TF1 *ga= new TF1("ga",fitN.Data(),-range,range);
2057  ga->SetParName(0,"Area Pi");
2058  ga->SetParName(1,"Mean Pi");
2059  ga->SetParName(2,"Sigma Pi");
2060  ga->SetParName(3,"A0");
2061  ga->SetParName(4,"A1");
2062  ga->SetParName(5,"A2");
2063  ga->SetParName(6,"A3");
2064  ga->SetParName(7,"A4");
2065  ga->SetParName(8,"A5");
2066 
2067  TCanvas *c = new TCanvas("Fit","Fit results");
2068  int i;
2069  TString name(hist->GetName());
2070  name += "MuFG";
2071  Int_t nBins = hist->GetXaxis()->GetNbins();
2072  Double_t xlow = hist->GetXaxis()->GetXmin();
2073  Double_t xup = hist->GetXaxis()->GetXmax();
2074  TH1D *MuF = new TH1D(name.Data(),"Avarage shift versus no. of measurement points",
2075  nBins, xlow, xup);
2076  name = hist->GetName();
2077  name += "SigmaFG";
2078  TH1D *SigmaF = new TH1D(name.Data(),"Sigma of z versus no. of measurement points",
2079  nBins, xlow, xup);
2080  Int_t i1 = 1, i2 = nBins;
2081  if (Ibin > 0 && Ibin <= nBins) {i1 = Ibin; i2 = Ibin;}
2082  Double_t XFitP, dXFitP, MuFitP,dMuFitP,SigmaFitP,dSigmaFitP;
2083  TH1 *proj = 0;
2084  for (i=i1; i<=i2; i++) {
2085  if (proj) delete proj;
2086  proj = hist->ProjectionY("proj",i,i);
2087  XFitP = hist->GetXaxis()->GetBinCenter(i);
2088  dXFitP = hist->GetXaxis()->GetBinWidth(i);
2089  Double_t chisq = -999;
2090  if (proj->Integral() < 1000) {
2091  if (N>0) goto NEXT;
2092  continue;
2093  }
2094  if (FitG(proj,g,ga)) goto NEXT;
2095  // proj->Fit("g","R");
2096  if (ga->GetParError(1) <= 0 || ga->GetParError(1) > 0.01 ) {
2097  printf("============= REJECT ================\n");
2098  goto NEXT;
2099  }
2100  MuFitP = ga->GetParameter(1);
2101  dMuFitP = ga->GetParError(1);
2102  // MuFitP = TMath::Log(MuFitP);
2103  SigmaFitP = ga->GetParameter(2);
2104  dSigmaFitP = ga->GetParError(2);
2105  chisq = ga->GetChisquare();
2106  if (chisq > 1.e4) goto NEXT;
2107  MuF->SetCellContent(i,0,MuFitP);
2108  MuF->SetCellError(i,0,dMuFitP);
2109  SigmaF->SetCellContent(i,0,SigmaFitP);
2110  SigmaF->SetCellError(i,0,dSigmaFitP);
2111  proj->Draw(); c->Update();
2112  printf("Bin: %i x: %f +/- %f MuF: %f+/-%f Sigma: %f+/-%f chisq: %f \n",
2113  i,XFitP,dXFitP,MuFitP,dMuFitP,SigmaFitP,dSigmaFitP,chisq);
2114  NEXT:
2115  N++;
2116  // delete proj;
2117  }
2118  if (! Ibin ) {
2119  c->Clear();
2120  // c->Divide(2,1);
2121  // c->cd(1);
2122  // MuF->Draw();
2123  // c->cd(2);
2124  SigmaF->SetMarkerStyle(20);
2125  TF1 *p = new TF1("p","[0]/pow(x,[1])");
2126  p->SetParameter(0,0.64);
2127  p->SetParameter(1,0.50);
2128  SigmaF->Fit("p","R");
2129  SigmaF->SetMinimum(0.06);
2130  SigmaF->Draw();
2131  }
2132  TString NewRootFile(hist->GetName());
2133  NewRootFile += fRootFile->GetName();
2134  TFile *f = new TFile(NewRootFile.Data(),"update");
2135  MuF->Write();
2136  SigmaF->Write();
2137  delete f;
2138  fRootFile->cd();
2139 }
2140 //________________________________________________________________________________
2141 void Fit4G(Int_t ng=2, Int_t hyp=-1, Int_t bin=0,
2142  Double_t xmin1=-1.0, Double_t xmax1 = 1.0,
2143  Double_t Mu2 = -.5,Double_t xmin2=-1., Double_t xmax2 = 1.,
2144  Double_t Mu3 = 0.0,Double_t xmin3=-1., Double_t xmax3 = 1.,
2145  Double_t Mu4 = 0.0){
2146  Double_t sigmas[2] = {0.06,0.12};
2147  const Int_t nHYPS = 4;
2148  TH2 *hists[nHYPS];
2149  TProfile *histp[nHYPS];
2150  TFile *fRootFile = (TFile *) gDirectory->GetFile();
2151  if (! fRootFile ) {printf("Cannot find/open %s",fRootFile->GetName()); return;}
2152  TFile *f = 0;
2153  if (bin == 0) {
2154  TString newfile("BBFit");
2155  newfile += gSystem->BaseName(fRootFile->GetName());
2156  f = new TFile(newfile.Data(),"update");
2157  }
2158  for (int i = 0; i< nHYPS; i++) {
2159  if (hyp >= 0 && hyp != i) continue;
2160  hists[i] = (TH2 *) fRootFile->Get(HistNames[i]);
2161  if (!hists[i]) {printf("Cannot histogram %s\n",HistNames[i]); return;}
2162  histp[i] = (TProfile *) fRootFile->Get(HistNameP[i]);
2163  if (!histp[i]) {printf("Cannot histogram %s\n",HistNameP[i]); return;}
2164  }
2165  TF1 *g1 = new TF1("g1","gaus",xmin1,xmax1);
2166  TF1 *g = 0, *g2 = 0, *g3 = 0, *g4 = 0, *ga = 0;
2167  TCanvas *canvas = new TCanvas("canvas","canvas");
2168  if (ng<0) {
2169  ga = new TF1("ga",Form("gaus(0)+exp(pol%i(3))",-ng),xmin2,xmax2);
2170  ga->SetParName(0,"Constant1");
2171  ga->SetParName(1,"Mean1");
2172  ga->SetParName(2,"Sigma1");
2173  ga->SetParName(3,"Const");
2174  for (int m = 0; m <= -ng; m++) ga->SetParName(m+3,Form("Slope%i",m));
2175  }
2176  else {
2177  for (int m = 2; m <= ng; m++) {
2178  if (m == 2) {g2 = new TF1("g2","gaus(0)+gaus(3)",xmin2,xmax2); g = g2;}
2179  if (m == 3) {g3 = new TF1("g3","gaus(0)+gaus(3)+gaus(6)",xmin3,xmax3); g = g3;}
2180  if (m == 4) {g4 = new TF1("g4","gaus(0)+gaus(3)+gaus(6)+gaus(9)",-1.,1.); g = g4;}
2181  for (int k = 0; k < 3*m; k += 3) {
2182  g->SetParName(k ,Form("Constant%i",m));
2183  g->SetParName(k+1,Form("Mean%i",m));
2184  g->SetParName(k+2,Form("Sigma%i",m));
2185  }
2186  }
2187  }
2188  TH1 *proj = 0;
2189  Double_t params[12];
2190  Int_t k,i;//,ibin;
2191  Int_t Bin = TMath::Abs(bin);
2192  for (k = 0; k<nHYPS; k++) {
2193  if (hyp != -1 && k != hyp) continue;
2194  TH2 *hist = hists[k];
2195  Int_t nx = hist->GetNbinsX();
2196  for (i=1; i<=nx; i++) {
2197  if (bin != 0 && Bin != i) continue;
2198  // if (proj) delete proj;
2199  Char_t line[40];
2200  sprintf(line,"%s_%i",hist->GetName(),i);
2201  TString name(line);
2202  proj = hist->ProjectionY(name.Data(),i,i);
2203  Int_t ix1=proj->GetXaxis()->FindBin(-.1);
2204  Int_t ix2=proj->GetXaxis()->FindBin(0.1);
2205  if (proj->Integral(ix1,ix2) < 100) {
2206  printf("hist:%s bin %i for hyp %i has only %10.0f entries\n",hist->GetName(),Bin,k,proj->Integral());
2207  delete proj;
2208  continue;
2209  }
2210  Int_t NFitPoints = 0;
2211  Double_t xrange = 0.6;
2212  Double_t chisq;
2213  g = g1;
2214  g->SetParLimits(0,0,1.e7);
2215  g->SetParLimits(1,-xrange,xrange);
2216  // g->SetParLimits(1,-.1,.1);
2217  g->SetParLimits(2,sigmas[0],sigmas[1]);
2218  proj->Fit(g->GetName(),"ri");
2219  g->GetParameters(params);
2220  Int_t kcase = -1;
2221  if (xmin1 < -0.5 && xmax1 > 0.5 && g->GetProb() > 1.e-7
2222  && TMath::Abs(g->GetParameter(1)) < 0.05) goto Done;
2223  if (ga) {
2224  g = ga;
2225  g1->GetParameters(params);
2226  params[1] = 0.;
2227  params[3] = 0.;
2228  params[4] = 0.;
2229  params[5] = 0.;
2230  params[6] = 0.;
2231  params[7] = 0.;
2232  params[8] = 0.;
2233  g->SetParameters(params);
2234  kcase = 0;
2235  if (g->GetProb() > 1.e-7) goto Done;
2236  }
2237  else {
2238  if (g2) {
2239  params[1] = 0.;
2240  params[3] = params[0];
2241  params[4] = Mu2;
2242  params[5] = 2.0*params[2];
2243  g = g2;
2244  g->SetParameters(params);
2245  g->SetParLimits(0,0,1.e7);
2246  g->SetParLimits(1,-xrange,xrange);
2247  // g->SetParLimits(1,-.1,.1);
2248  g->SetParLimits(2,sigmas[0],sigmas[1]);
2249  g->SetParLimits(3,0,1.e7);
2250  g->SetParLimits(5,sigmas[0],sigmas[1]);
2251  proj->Fit(g->GetName(),"ri");
2252  g->GetParameters(params);
2253  kcase = 2;
2254  if (g->GetProb() > 1.e-7) goto Done;
2255  if (g3) {
2256  params[1] = 0.;
2257  params[6] = params[0];
2258  params[7] = Mu3;
2259  params[8] = 2.0*params[2];
2260  g = g3;
2261  g->SetParameters(params);
2262  g->SetParLimits(0,0,1.e7);
2263  g->SetParLimits(1,-xrange,xrange);
2264  // g->SetParLimits(1,-.1,.1);
2265  g->SetParLimits(2,sigmas[0],sigmas[1]);
2266  g->SetParLimits(3,0,1.e7);
2267  g->SetParLimits(5,sigmas[0],sigmas[1]);
2268  g->SetParLimits(6,0,1.e7);
2269  g->SetParLimits(8,sigmas[0],sigmas[1]);
2270  proj->Fit(g->GetName(),"ri");
2271  g->GetParameters(params);
2272  kcase = 3;
2273  if (g->GetProb() > 1.e-7) goto Done;
2274  if (g4) {
2275  params[1] = 0.;
2276  params[9] = params[0];
2277  params[10] = Mu4;
2278  params[11] = 2.0*params[2];
2279  g = g4;
2280  g->SetParameters(params);
2281  g->SetParLimits(0,0,1.e7);
2282  g->SetParLimits(1,-xrange,xrange);
2283  // g->SetParLimits(1,-.1,.1);
2284  g->SetParLimits(2,sigmas[0],sigmas[1]);
2285  g->SetParLimits(3,0,1.e7);
2286  g->SetParLimits(5,sigmas[0],sigmas[1]);
2287  g->SetParLimits(6,0,1.e7);
2288  g->SetParLimits(8,sigmas[0],sigmas[1]);
2289  g->SetParLimits(9,0,1.e7);
2290  g->SetParLimits(11,sigmas[0],sigmas[1]);
2291  proj->Fit(g->GetName(),"ri");
2292  kcase = 4;
2293  g->GetParameters(params);
2294  }
2295  }
2296  }
2297  }
2298  Done:
2299  if (g) {
2300  proj->Fit(g->GetName(),"RIM");
2301  canvas->Update();
2302  Int_t l = 1;
2303  Double_t mu = g->GetParameter(l);
2304  for (int m=2;m<=kcase;m++) {
2305  if (TMath::Abs(mu) > TMath::Abs(g->GetParameter(3*m-2))) {
2306  l = 3*m-2; mu = g->GetParameter(l);
2307  // printf("l=%i\n",l);
2308  }
2309  }
2310  Nu[N] = g->GetParameter(l);// printf("l=%i Nu=%f\n",l,Nu[N]);
2311  Mu[N] = Nu[N] + histp[k]->GetBinContent(i);
2312  dMu[N] = g->GetParError(l);
2313  // Mu[N] = TMath::Log(Mu[N]);
2314  Sigma[N] = g->GetParameter(2);
2315  dSigma[N] = g->GetParError(2);
2316  NFitPoints = g->GetNumberFitPoints();
2317  Int_t NDF = g->GetNDF();
2318  Double_t prob = g->GetProb();//TMath::Prob(chisq, NDF);
2319  chisq = g->GetChisquare();
2320  X[N] = hist->GetXaxis()->GetBinCenter(i);
2321  dX[N] = hist->GetXaxis()->GetBinWidth(i);
2322  Double_t pionM = Masses[k]*pow(10.,X[N]);
2323  printf ("%s :hyp = %i bin=%i, Point=%i, x=%f, p=%f, Delta_I=%f, I=%f, Sigma_I=%f,\n"
2324  "chisq=%f, NoPoints=%i,ndf=%i, prob=%f\n",
2325  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob);
2326  // if (bin < 0 && prob > 1.e-7 && TMath::Abs(Nu[N]) < 0.05) {
2327  // if ( bin >= 0 && prob > 1.e-7 && TMath::Abs(Nu[N]) < 0.10 || bin < 0) {
2328  if ( (bin >= 0 && TMath::Abs(Nu[N]) < 0.10) || bin < 0) {
2329  printf("{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
2330  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
2331  TString FileN("FitPars");
2332  if (hyp > -1) FileN += HistNames[hyp];
2333  FileN += ".h";
2334  FILE *fp = fopen(FileN.Data(),"a");
2335  if (fp) {
2336  if (N == 0) {
2337  TDatime time;
2338  fprintf(fp,"// Date: Time = %i : %i\n",time.GetDate(), time.GetTime());
2339  fprintf(fp,
2340 "// bin, Point, x, p, Delta_I, I, Sigma_I, chisq, NoPoints,ndf, prob\n");
2341  }
2342  fprintf(fp,
2343  "{\"%-4s\",%2i,%4i,%6i,%6.3f,%7.3f,%10.6f,%10.6f,%8.5f,%10.3f},//%3i,%3i,%5.3f -- %s\n",
2344  Names[k],k,i,N,X[N],pionM,Nu[N],Mu[N],dMu[N],chisq,NFitPoints,NDF,prob,g->GetName());
2345  fclose(fp);
2346  }
2347  proj->Write();
2348  }
2349  else printf ("================== Skip it\n");
2350  }
2351  N++;
2352  proj->Draw();// cnew->Update();
2353  }
2354  }
2355  if (f) {delete f;}
2356 }
2357 //#define DEBUG
2358 //________________________________________________________________________________
2359 TList *ListOfKeys() {
2360  TList *list = 0;
2361  TSeqCollection *files = gROOT->GetListOfFiles();
2362  TIter next(files);
2363  TFile *f = 0;
2364  while ((f = (TFile *) next())) {
2365  f->cd();
2366  list = f->GetListOfKeys();
2367  break;
2368  }
2369  return list;
2370 }
2371 //________________________________________________________________________________
2372 TH1 *FindHistograms(const Char_t *Name, const Char_t *fit = "") {
2373  TH1 *hist = 0;
2374 #ifdef DEBUG
2375  cout << "FindHistograms(\"" << Name << "\",\"" << fit << "\")" << endl;
2376 #endif
2377  TSeqCollection *files = gROOT->GetListOfFiles();
2378  TIter next(files);
2379  TFile *f = 0;
2380  TString Fit(fit);
2381  if (Fit != "") {
2382  while ((f = (TFile *) next())) {
2383 #ifdef DEBUG
2384  cout << "look in " << f->GetName() << endl;
2385 #endif
2386  TString fName(gSystem->BaseName(f->GetName()));
2387  if (fName.BeginsWith(Name)) {
2388  hist = (TH1 *) f->Get(fit);
2389 #ifdef DEBUG
2390  if (hist) cout << "Found histogram " << hist->GetName() << " in file " << f->GetName() << endl;
2391 #endif
2392  return hist;
2393  }
2394  }
2395  } else {
2396  while ((f = (TFile *) next())) {
2397  f->cd();
2398 #ifdef DEBUG
2399  cout << "look in " << f->GetName() << endl;
2400 #endif
2401  hist = (TH1 *) f->Get(Name);
2402  if (hist) {
2403 #ifdef DEBUG
2404  cout << "Found histogram " << hist->GetName() << " in file " << f->GetName() << endl;
2405 #endif
2406  return hist;
2407  }
2408  }
2409  }
2410  return hist;
2411 }
2412 //________________________________________________________________________________
2413 void DrawSummary0(TFile *f, const Char_t *opt) {
2414  if (! f) return;
2415  gStyle->SetTimeOffset(788936400);
2416  f->cd();
2417  TList *keys = f->GetListOfKeys();
2418  if (! keys) return;
2419  keys->Sort();
2420  TIter next(keys);
2421  TKey *key = 0;
2422  TObjArray hists;
2423  TString Opt(opt);
2424  TString Title;
2425  TF1 *powfit = new TF1("powfit","[0]*pow(x,[1])",40,80);
2426  powfit->SetParameters(0.5,-0.5);
2427  while ((key = (TKey *) next())) {
2428  TH1 *hist = FindHistograms(key->GetName());
2429  if (! hist) continue;
2430  if (hist->GetEntries() < 1) continue;
2431  if (Opt != "") {
2432  TString name(hist->GetName());
2433  if (! name.Contains(Opt.Data())) continue;
2434  }
2435  hists.AddLast(hist);
2436  TString hName(hist->GetName());
2437  if ((( hist->IsA()->InheritsFrom( "TH3C" ) ||
2438  hist->IsA()->InheritsFrom( "TH3S" ) ||
2439  hist->IsA()->InheritsFrom( "TH3I" ) ||
2440  hist->IsA()->InheritsFrom( "TH3F" ) ||
2441  hist->IsA()->InheritsFrom( "TH3D" ) ) ||
2442  ( hist->IsA()->InheritsFrom( "TH2C" ) ||
2443  hist->IsA()->InheritsFrom( "TH2S" ) ||
2444  hist->IsA()->InheritsFrom( "TH2I" ) ||
2445  hist->IsA()->InheritsFrom( "TH2F" ) ||
2446  hist->IsA()->InheritsFrom( "TH2D" ) ))&&
2447  ( ! (hName.BeginsWith("Points") || hName.BeginsWith("TPoints") || hName.BeginsWith("MPoints")) ) ) {
2448  TH1 *mu = FindHistograms(key->GetName(),"mu");
2449  if (! mu) continue;
2450  mu->SetName(Form("%s_mu",hist->GetName()));
2451  hists.AddLast(mu);
2452  TH1 *sigma = FindHistograms(key->GetName(),"sigma");
2453  if (! sigma) continue;
2454  sigma->SetName(Form("%s_sigma",hist->GetName()));
2455  hists.AddLast(sigma);
2456  }
2457  }
2458  Int_t N = hists.GetEntriesFast();
2459  Int_t nx = (Int_t) (TMath::Sqrt(N));
2460  if (nx*nx < N) nx++;
2461  Int_t ny = N/nx;
2462  if (nx*ny < N) ny++;
2463  cout << "N " << N << " nx " << nx << " ny " << ny << endl;
2464  TString Tag(gSystem->BaseName(f->GetName()));
2465  Tag.ReplaceAll(".root","");
2466  Tag += opt;
2467  // TCanvas *c1 = new TCanvas(Tag,Tag,200,10,700,780);
2468  Double_t dx = 0.98/nx;
2469  Double_t dy = 0.98/ny;
2470  TObjArray pads(N);
2471  for (Int_t iy = 0; iy < ny; iy++) {
2472  for (Int_t ix = 0; ix < nx; ix++) {
2473  Int_t i = ix + nx*iy;
2474  if (i >= N) continue;
2475  TH1 *hist = (TH1 *) hists[i];
2476  if (! hist) continue;
2477  Double_t x1 = 0.01 + dx* ix;
2478  Double_t x2 = 0.01 + dx*(ix+1);
2479  Double_t y1 = 0.99 - dy*(iy+1);
2480  Double_t y2 = 0.99 - dy* iy;
2481  TPad *pad = new TPad(Form("pad_%s",hist->GetName()),hist->GetTitle(),x1,y1,x2,y2);
2482  // cout << "Create pad " << pad->GetName() << "/" << pad->GetTitle() << endl;
2483  pad->Draw();
2484  pads.AddLast(pad);
2485  }
2486  }
2487  for (Int_t iy = 0; iy < ny; iy++) {
2488  for (Int_t ix = 0; ix < nx; ix++) {
2489  Int_t i = ix + nx*iy;
2490  if (i >= N) continue;
2491  TH1 *hist = (TH1 *) hists[i];
2492  if (! hist) continue;
2493  TPad *pad = (TPad *) pads[i];
2494  if (! pad) continue;
2495  pad->cd();
2496  TProfile *prof = 0;
2497  if ( hist->IsA()->InheritsFrom( "TProfile" ) ) prof = (TProfile *) hist;
2498  Int_t NbinsX = hist->GetNbinsX();
2499  Int_t xmin = NbinsX;
2500  Int_t xmax = 0;
2501  Double_t ymin = 1e9;
2502  Double_t ymax = -1e9;
2503  // cout << "Draw pad " << pad->GetName() << "/" << pad->GetTitle() << "\t x " << ix << " y " << iy << endl;
2504  //3D
2505  if ( hist->IsA()->InheritsFrom( "TH3C" ) ||
2506  hist->IsA()->InheritsFrom( "TH3S" ) ||
2507  hist->IsA()->InheritsFrom( "TH3I" ) ||
2508  hist->IsA()->InheritsFrom( "TH3F" ) ||
2509  hist->IsA()->InheritsFrom( "TH3D" ) ) {
2510  ((TH3 *)hist)->Project3DProfile("xy")->Draw("colz");
2511  goto TIMEAxis;
2512  }
2513  //2D
2514  if ( hist->IsA()->InheritsFrom( "TH2C" ) ||
2515  hist->IsA()->InheritsFrom( "TH2S" ) ||
2516  hist->IsA()->InheritsFrom( "TH2I" ) ||
2517  hist->IsA()->InheritsFrom( "TH2F" ) ||
2518  hist->IsA()->InheritsFrom( "TH2D" ) ) {
2519  if (hist->GetMaximum() > 0 && hist->GetMinimum() >= 0) pad->SetLogz(1);
2520  TString hName(hist->GetName());
2521  if (hName.BeginsWith("Points") || hName.BeginsWith("TPoints") || hName.BeginsWith("MPoints")) {
2522  Title = hist->GetTitle();
2523  Title.ReplaceAll("/sigma","");
2524  hist->SetTitle(Title);
2525  TAxis *y = hist->GetYaxis();
2526  Int_t iy1 = y->FindBin(-0.15);
2527  Int_t iy2 = y->FindBin( 0.50);
2528  y->SetRange(iy1,iy2);
2529  ((TH2 *)hist)->Draw("colz");
2530  TH1 *mu = FindHistograms(hist->GetName(),"mu");
2531  if (! mu) goto TIMEAxis;
2532  mu->SetMarkerColor(1);
2533  mu->SetMarkerStyle(20);
2534  mu->Draw("same");
2535  mu->Fit("pol0","e0r","",10,120);
2536  TLegend *leg = new TLegend(0.25,0.6,0.9,0.9,"");
2537  TF1 *f = mu->GetFunction("pol0");
2538  if (f) {
2539  f->Draw("same");
2540  Title = Form("#mu = %5.2f +/- %5.2f %\%",100*f->GetParameter(0),100*f->GetParError(0));
2541  cout << Title << endl;
2542  leg->AddEntry(mu,Title.Data());
2543  }
2544  TH1 *sigma = FindHistograms(hist->GetName(),"sigma");
2545  if (! sigma) goto TIMEAxis;
2546  sigma->SetMarkerColor(1);
2547  sigma->SetMarkerStyle(21);
2548  sigma->Draw("same");
2549  sigma->Fit(powfit,"r0");
2550  f = sigma->GetFunction("powfit");
2551  if (f) {
2552  f->Draw("same");
2553  Title = Form("#sigma(@76cm) = %5.2f%\%",100*f->Eval(76));
2554  cout << Title << endl;
2555  leg->AddEntry(sigma,Title.Data());
2556  }
2557  leg->Draw();
2558  } else hist->Draw("colz");
2559  goto TIMEAxis;
2560  }
2561  //1D + Prof
2562  NbinsX = hist->GetNbinsX();
2563  xmin = NbinsX;
2564  xmax = 0;
2565  ymin = 1e9;
2566  ymax = -1e9;
2567  for (Int_t i = 1; i <= NbinsX; i++) {
2568  Double_t y = hist->GetBinContent(i);
2569  if ((prof && prof->GetBinEntries(i)) || y > 0) {
2570  Int_t x = i;
2571  xmin = TMath::Min(xmin,x);
2572  xmax = TMath::Max(xmax,x);
2573  ymin = TMath::Min(ymin,y);
2574  ymax = TMath::Max(ymax,y);
2575  }
2576  }
2577  if (ymin < ymax) {
2578  hist->SetMaximum(1.1*ymax);
2579  hist->SetMinimum(0.9*ymin);
2580  // cout << "Set min/max for " << hist->GetName() << "\t" << hist->GetMinimum() << "/" << hist->GetMaximum() << endl;
2581  }
2582  if (xmin < xmax) hist->GetXaxis()->SetRange(xmin,xmax);
2583  hist->Draw();
2584  TIMEAxis:
2585  TAxis *xax = hist->GetXaxis();
2586  if (xax->GetBinLowEdge(1) > 1e6) {
2587  xax->SetTimeDisplay(1);
2588  gPad->Modified();
2589  }
2590  }
2591  }
2592 }
2593 //________________________________________________________________________________
2594 void DrawSummary(const Char_t *opt="") {
2595  TSeqCollection *files = gROOT->GetListOfFiles();
2596  TIter next(files);
2597  TFile *f = 0;
2598  TString FName("");
2599  while ((f = (TFile *) next())) {
2600  FName = gSystem->BaseName(f->GetName());
2601  if (FName.BeginsWith("Hist")) break;
2602  }
2603  if (! f) {
2604  cout << "Hist* root file has not been found" << endl;
2605  return;
2606  }
2607  DrawSummary0(f,opt);
2608 }
2609 //________________________________________________________________________________
2610 Double_t gNFunc(Double_t *x=0, Double_t *par=0) {
2611  static Double_t sigmaOLD = -1;
2612  static Double_t fracpiOld = -1;
2613  if (! x || ! par) {
2614  sigmaOLD = -1;
2615  fracpiOld = -1;
2616  return 0;
2617  }
2618  // par[0] - norm
2619  // par[1] - pion position wrt Z_pion (Bichsel prediction)
2620  // par[2] - sigma
2621  // par[3] - proton signal
2622  // par[4] - Kaon -"-
2623  // par[5] - electorn -"-
2624  // par[6] - deuteron -"-
2625  // par[7] - Total
2626 #ifdef __HEED_MODEL__
2627  // par[8] - case
2628  // par[9] - scale
2629  // par[10]- row
2630  Int_t row = par[10];
2631  StdEdxModel::ESector kTpcOuterInner = StdEdxModel::kTpcOuter;
2632  if (row > 0 && row <= 13) kTpcOuterInner = StdEdxModel::kTpcInner;
2633 #endif /* __HEED_MODEL__ */
2634  // ratio dN/dx_h / dN/dx_pi: P K e d
2635  static Double_t dNdxR[4] = {1.97273, 1.32040, 1.16313, 3.42753};
2636  static Int_t _debug = 0;
2637  static TCanvas *c1 = 0;
2638  if (_debug) {
2639  c1 = (TCanvas *) gROOT->GetListOfCanvases()->FindObject("c1");
2640  if (! c1) c1 = new TCanvas();
2641  else c1->Clear();
2642  c1->cd();
2643  }
2644  enum {NFIT_HYP = 5, NT=6};
2645  static TH1D *hists[6] = {0};
2646  static Int_t icaseOLD = -1;
2647  static Double_t meanPion = -1, RMSPion = -1, mpvPion = -1;
2648  static Double_t ln10 = TMath::Log(10.);
2649  Double_t sigma = par[2];
2650  if (sigma != sigmaOLD) {
2651 #ifndef __HEED_MODEL__
2652  TF1 *zdE = StdEdxModel::instance()->zdEdx();
2653 #else /* __HEED_MODEL__ */
2654  TF1 *zdE = StdEdxModel::instance()->zdEdx(kTpcOuterInner);
2655 #endif /* __HEED_MODEL__ */
2656  sigmaOLD = sigma;
2657  fracpiOld = -1;
2658  static const Char_t *names[5] = {"pi","P","K","e","d"};
2659  for (Int_t i = 0; i < NFIT_HYP; i++) {
2660  if (hists[i]) hists[i]->Reset();
2661  else {
2662  hists[i] = new TH1D(Form("dE%s",names[i]),Form("Expected dE for %s",names[i]),100,-5,5);
2663  hists[i]->SetMarkerColor(i+1);
2664  hists[i]->SetLineColor(i+1);
2665  }
2666  Double_t entries = projNs[i]->GetEntries();
2667  Double_t mean = projNs[i]->GetMean();
2668 #ifndef __HEED_MODEL__
2669  Double_t mpv = StdEdxModel::instance()->zMPV()->Eval(mean,sigma);
2670 #else /* __HEED_MODEL__ */
2671  Double_t mpv = StdEdxModel::instance()->zMPV(kTpcOuterInner)->Eval(mean,sigma);
2672 #endif /* __HEED_MODEL__ */
2673  Double_t RMS = projNs[i]->GetRMS();
2674  if (i == 0) {
2675  meanPion = mean;
2676  mpvPion = mpv;
2677  RMSPion = RMS;
2678  }
2679  Int_t nx = projNs[i]->GetNbinsX();
2680  for (Int_t ix = 1; ix <= nx; ix++) {
2681  Double_t v = projNs[i]->GetBinContent(ix);
2682  if (v > 0.0) {
2683 #ifndef __HEED_MODEL__
2684  Double_t n_PL10 = projNs[i]->GetBinCenter(ix);
2685  Double_t n_P = TMath::Exp(n_PL10*ln10);
2686 #else /* __HEED_MODEL__ */
2687  Double_t n_PL = projNs[i]->GetBinCenter(ix);
2688  Double_t n_P = TMath::Exp(n_PL);
2689 #endif /* __HEED_MODEL__ */
2690  Double_t Sigma = TMath::Sqrt(sigma*sigma + 1./n_P);
2691  zdE->SetParameter(1, n_P);
2692  zdE->SetParameter(2,mpv - mpvPion);
2693  zdE->SetParameter(3,Sigma);
2694  hists[i]->Add(zdE,v/entries);
2695  }
2696  }
2697  if (c1) {
2698  hists[i]->Draw();
2699  c1->Update();
2700  }
2701  }
2702  if (! hists[NFIT_HYP]) hists[NFIT_HYP] = new TH1D("dEAll","Expected dE for All",100,-5,5);
2703  }
2704  Double_t frac[NFIT_HYP] = {0};
2705  static Double_t fracOld[6] = {-1};
2706  Int_t i;
2707  frac[0] = 1;
2708  Bool_t updatedFractions = kFALSE;
2709  for (i = 1; i < NFIT_HYP; i++) {
2710  frac[i] = TMath::Sin(par[2+i]);
2711  frac[i] *= frac[i];
2712  if (TMath::Abs(frac[i]-fracOld[i]) > 1e-7) {
2713  fracOld[i] = frac[i];
2714  updatedFractions = kTRUE;
2715  }
2716  frac[0] -= frac[i];
2717  }
2718  if (fracOld[0] != frac[0]) updatedFractions = kTRUE;
2719  Int_t icase = (Int_t) par[8];
2720  Int_t i1 = 0;
2721  Int_t i2 = NFIT_HYP - 1;
2722  if (icase >= 0) {i1 = i2 = icase;}
2723  if (icase != icaseOLD || updatedFractions) {
2724  icaseOLD = icase;
2725  fracpiOld = frac[0];
2726  hists[NFIT_HYP]->Reset();
2727  for (i = i1; i <= i2; i++) {
2728  if (frac[i] > 1e-7) {
2729  hists[NFIT_HYP]->Add(hists[i], frac[i]);
2730  }
2731  }
2732  if (c1) {
2733  hists[NFIT_HYP]->Draw();
2734  c1->Update();
2735  }
2736  }
2737 #ifndef __HEED_MODEL__
2738  Double_t Value = par[7]*TMath::Exp(par[0])*hists[NFIT_HYP]->Interpolate(x[0]-par[1]);
2739 #else /* __HEED_MODEL__ */
2740  Double_t Value = par[7]*TMath::Exp(par[0])*hists[NFIT_HYP]->Interpolate(par[9]*(x[0]-par[1]));
2741 #endif /* __HEED_MODEL__ */
2742  return Value;
2743 }
2744 //________________________________________________________________________________
2745 #ifndef __HEED_MODEL__
2746 TF1 *FitNF(TH1 *proj, Option_t *opt) {// fit with no. of primary clusters
2747 #else /* __HEED_MODEL__ */
2748 TF1 *FitNF(TH1 *proj, Option_t *opt, Int_t row) {// fit with no. of primary clusters
2749 #endif /* __HEED_MODEL__ */
2750  for (Int_t ih = 0; ih < 5; ih++) {
2751  if (! projNs[ih]) {
2752  cout << "projNs[" << ih << "] is not defined. Abort." << endl;
2753  return 0;
2754  }
2755  }
2756  static TSpectrum *fSpectrum = 0;
2757  if (! fSpectrum) {
2758  fSpectrum = new TSpectrum(6);
2759  }
2760  // fit in momentum range p = 0.45 - 0.50 GeV/c
2761  if (! proj) return 0;
2762  TString Opt(opt);
2763  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
2764  TF1 *g2 = (TF1*) gROOT->GetFunction("GN");
2765  enum {NFIT_HYP = 5}; // ignore e and d
2766  if (! g2) {
2767 #ifndef __HEED_MODEL__
2768  g2 = new TF1("GN",gNFunc, -5, 5, 9);
2769  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
2770 #else /* __HEED_MODEL__ */
2771  g2 = new TF1("GN",gNFunc, -5, 5, 11);
2772  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
2773 #endif /* __HEED_MODEL__ */
2774  g2->SetParName(1,"mu"); g2->SetParLimits(1,-2.5,2.5);
2775  g2->SetParName(2,"Sigma"); g2->FixParameter(2, 0); // g2->SetParLimits(2,0.,0.5);
2776  g2->SetParName(3,"P"); g2->SetParLimits(3,0.0,TMath::Pi()/2);
2777  g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.30);
2778  g2->SetParName(5,"e"); g2->SetParLimits(5,0.0,0.30); g2->FixParameter(5,0);
2779  g2->SetParName(6,"d"); g2->SetParLimits(6,0.0,0.30); g2->FixParameter(6,0);
2780  // g2->SetParName(6,"ScaleL"); g2->SetParLimits(6,-2.5,2.5);
2781  g2->SetParName(7,"Total");
2782  g2->SetParName(8,"Case");
2783 #ifdef __HEED_MODEL__
2784  g2->SetParName(9,"Scale"); g2->SetParameter(9,1); g2->SetParLimits(9,0.5,1.5);
2785  g2->SetParName(10,"row");
2786 #endif /* __HEED_MODEL__ */
2787  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
2788  }
2789  // Find pion peak
2790  Int_t nfound = fSpectrum->Search(proj);
2791  if (nfound < 1) return 0;
2792  Int_t npeaks = 0;
2793  Float_t *xpeaks = fSpectrum->GetPositionX();
2794  Float_t xp = 0;
2795  if (nfound > 2) nfound = 2;
2796  Double_t xpi = 9999;
2797  for (Int_t p = 0; p < nfound; p++) {
2798  xp = xpeaks[p];
2799  Int_t bin = proj->GetXaxis()->FindBin(xp);
2800  Double_t yp = proj->GetBinContent(bin);
2801  Double_t ep = proj->GetBinError(bin);
2802  if (yp-5*ep < 0) continue;
2803  if (xp < xpi) xpi = xp;
2804  }
2805  Double_t total = proj->Integral()*proj->GetBinWidth(5);
2806 #ifndef __HEED_MODEL__
2807  g2->SetParameters(0, xpi, 0.0, 0.0, 0.1, 0.1, 0.0,0.0,-1.);
2808 #else /* __HEED_MODEL__ */
2809  g2->SetParameters(0, xpi, 0.0, 0.0, 0.1, 0.1, 0.0,0.0,-1., 1., row);
2810 #endif /* __HEED_MODEL__ */
2811  g2->FixParameter(2, 0);
2812 #ifdef __HEED_MODEL__
2813  g2->FixParameter(3, 0);
2814  g2->FixParameter(4, 0);
2815 #endif /* __HEED_MODEL__ */
2816  g2->FixParameter(5,0);
2817  g2->FixParameter(6,0);
2818  g2->FixParameter(7,total);
2819  g2->FixParameter(8,-1);
2820 #ifdef __HEED_MODEL__
2821  g2->SetParameter(9,1); //TMath::Log(10.));
2822  g2->FixParameter(10, row);
2823 #endif /* __HEED_MODEL__ */
2824  gNFunc();
2825  TFitResultPtr res = proj->Fit(g2,Opt.Data());
2826  Int_t iok = res->Status();
2827  if ( iok < 0) {
2828  cout << g2->GetName() << " fit has failed with " << iok << " for "
2829  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
2830  proj->Fit(g2,Opt.Data());
2831  }
2832  // g2->ReleaseParameter(5); g2->SetParameter(5,0.01); g2->SetParLimits(5,0.0,0.30);
2833  // g2->ReleaseParameter(6); g2->SetParameter(6,0.01); g2->SetParLimits(6,0.0,0.30);
2834  Opt += "m";
2835  gNFunc();
2836  res = proj->Fit(g2,Opt.Data());
2837  iok = res->Status();
2838  if (iok < 0 ) return 0;
2839  if (! Opt.Contains("q",TString::kIgnoreCase)) {
2840  gNFunc();
2841  proj->Draw();
2842  Double_t params[10];
2843  g2->GetParameters(params);
2844  Double_t X = params[1];
2845  Double_t Y = TMath::Exp(params[0]);
2846  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
2847  proj->GetListOfFunctions()->Add(pm);
2848  pm->SetMarkerStyle(23);
2849  pm->SetMarkerColor(kRed);
2850  pm->SetMarkerSize(1.3);
2851  for (int i = 0; i < NFIT_HYP; i++) {
2852  TF1 *f = new TF1(*g2);
2853  f->SetName(Peaks[i].Name);
2854  f->FixParameter(8,i);
2855  f->SetLineColor(i+2);
2856  gNFunc();
2857  f->Draw("same");
2858  proj->GetListOfFunctions()->Add(f);
2859  }
2860  }
2861  return g2;
2862 }
2863 //________________________________________________________________________________
2864 Double_t gXFunc(Double_t *x=0, Double_t *par=0) {
2865  // par[0] - norm
2866  // par[1] - pion position wrt Z_pion (Bichsel prediction)
2867  // par[2] - sigma
2868  // par[3] - proton signal
2869  // par[4] - Kaon -"-
2870  // par[5] - electorn -"-
2871  // par[6] - ScaleL
2872  // par[7] - Total
2873  // par[8] - Case
2874  // par[9] - dX
2875  // par[10]- ddX
2876  // ratio dN/dx_h / dN/dx_pi: P K e d
2877  // static Double_t dNdxR[4] = {1.97273, 1.32040, 1.16313, 3.42753};
2878  // static Double_t dNdxL10[4] = {0.295068, 0.120707, 0.0656301, 0.534982};
2879  static Double_t dNdxL10[4] = {0.590131, 0.241411, 0.131261, 1.06996};
2880  static Int_t _debug = 0;
2881  static TCanvas *c1 = 0;
2882  if (_debug) {
2883  c1 = (TCanvas *) gROOT->GetListOfCanvases()->FindObject("c1");
2884  if (! c1) c1 = new TCanvas();
2885  else c1->Clear();
2886  c1->cd();
2887  }
2888  enum {NFIT_HYP = 3, NT};
2889  Double_t frac[NT];
2890  static Double_t sigmaOLD = -1;
2891  static Double_t fracpiOld = -1;
2892  static Int_t icaseOLD = -1;
2893  static Double_t dNPion = -1, mpvPion = -1;
2894  static Double_t ln10 = TMath::Log(10.);
2895  static TH1D *hists[5] = {0};
2896  if (! x || ! par) {
2897  sigmaOLD = -1;
2898  return 0;
2899  }
2900  Double_t norm = par[0];
2901  Double_t mu = par[1];
2902  Double_t sigma = par[2];
2903  Double_t ScaleL= par[6];
2904  Double_t dX = par[9];
2905  Double_t ddX = par[10];
2906  Int_t i;
2907  frac[0] = 1;
2908  for (i = 1; i <= NFIT_HYP; i++) {
2909  frac[i] = TMath::Sin(par[2+i]);
2910  frac[i] *= frac[i];
2911  frac[0] -= frac[i];
2912  }
2913  TF1 *zdE = StdEdxModel::instance()->zdEdx();
2914  if (sigma != sigmaOLD) {
2915  sigmaOLD = sigma;
2916  fracpiOld = -1;
2917  static const Char_t *names[5] = {"pi","P","K","e","d"};
2918  static Double_t dNdx_pion = 29.22; // 1/cm
2919  Double_t dNPion = dNdx_pion*dX;
2920  for (i = 0; i < NFIT_HYP; i++) {
2921  if (hists[i]) hists[i]->Reset();
2922  else {
2923  hists[i] = new TH1D(Form("dE%s",names[i]),Form("Expected dE for %s",names[i]),100,-5,5);
2924  hists[i]->SetMarkerColor(i+1);
2925  hists[i]->SetLineColor(i+1);
2926  }
2927  Double_t dN = dNPion;
2928  if (i) dN *= TMath::Power(10.,dNdxL10[i-1]);
2929  Double_t Sigma = TMath::Sqrt(sigma*sigma + 1./dN + (ddX/dX)*(ddX/dX));
2930  // if (i) Sigma = TMath::Sqrt(Sigma*Sigma + 0.15*0.15);
2931  Double_t mpv = StdEdxModel::instance()->zMPV()->Eval(TMath::Log10(dN),Sigma);
2932  if (i == 0) {
2933  mpvPion = mpv;
2934  }
2935  zdE->SetParameter(1, dN);
2936  zdE->SetParameter(2,mpv - mpvPion);
2937  zdE->SetParameter(3,Sigma);
2938  hists[i]->Add(zdE);
2939  if (c1) {
2940  hists[i]->Draw();
2941  c1->Update();
2942  }
2943  }
2944  if (! hists[NFIT_HYP]) hists[NFIT_HYP] = new TH1D("dEAll","Expected dE for All",100,-5,5);
2945  }
2946  Int_t icase = (Int_t) par[8];
2947  Int_t i1 = 0;
2948  Int_t i2 = NFIT_HYP;
2949  if (icase >= 0) {i1 = i2 = icase;}
2950  if (icase != icaseOLD || icase >= 0 || TMath::Abs(fracpiOld - frac[0]) > 1e-7) {
2951  icaseOLD = icase;
2952  fracpiOld = frac[0];
2953  hists[NFIT_HYP]->Reset();
2954  for (i = i1; i <= i2; i++) {
2955  if (frac[i] > 1e-7) {
2956  hists[NFIT_HYP]->Add(hists[i], frac[i]);
2957  }
2958  }
2959  if (c1) {
2960  hists[NFIT_HYP]->Draw();
2961  c1->Update();
2962  }
2963  }
2964  Double_t Value = par[7]*TMath::Exp(par[0])*hists[NFIT_HYP]->Interpolate(x[0]-par[1]);
2965  return Value;
2966 }
2967 //________________________________________________________________________________
2968 TF1 *FitXF(TH1 *proj, Option_t *opt, Double_t dX = 1.25, Double_t ddX = 0.01) {// fit with no. of primary clusters
2969  static TSpectrum *fSpectrum = 0;
2970  if (! fSpectrum) {
2971  fSpectrum = new TSpectrum(6);
2972  }
2973  // fit in momentum range p = 0.45 - 0.50 GeV/c
2974  if (! proj) return 0;
2975  TString Opt(opt);
2976  // Bool_t quet = Opt.Contains("Q",TString::kIgnoreCase);
2977  TF1 *g2 = (TF1*) gROOT->GetFunction("GX");
2978  enum {NFIT_HYP = 3}; // ignore e and d
2979  if (! g2) {
2980  g2 = new TF1("GX",gXFunc, -5, 5, 11);
2981  g2->SetParName(0,"norm"); g2->SetParLimits(0,-80,80);
2982  g2->SetParName(1,"mu"); g2->SetParLimits(1,-2.5,2.5);
2983  g2->SetParName(2,"Sigma"); g2->SetParLimits(2,0.,0.5);
2984  g2->SetParName(3,"P"); g2->SetParLimits(3,0.0,TMath::Pi()/2);
2985  g2->SetParName(4,"K"); g2->SetParLimits(4,0.0,0.30);
2986  g2->SetParName(5,"e"); g2->FixParameter(5,0);
2987  g2->SetParName(6,"ScaleL"); g2->FixParameter(6,0);
2988  g2->SetParName(7,"Total");
2989  g2->SetParName(8,"Case");
2990  g2->SetParName(9,"dX"); g2->FixParameter(9,dX);
2991  g2->SetParName(10,"ddX"); g2->FixParameter(10,ddX);
2992  // g2->SetParName(7,"factor"); g2->SetParLimits(7,-.1,0.1);
2993  }
2994  // Find pion peak
2995  Int_t nfound = fSpectrum->Search(proj);
2996  if (nfound < 1) return 0;
2997  Int_t npeaks = 0;
2998  Float_t *xpeaks = fSpectrum->GetPositionX();
2999  Float_t xp = 0;
3000  if (nfound > 2) nfound = 2;
3001  Double_t xpi = 9999;
3002  for (Int_t p = 0; p < nfound; p++) {
3003  xp = xpeaks[p];
3004  Int_t bin = proj->GetXaxis()->FindBin(xp);
3005  Double_t yp = proj->GetBinContent(bin);
3006  Double_t ep = proj->GetBinError(bin);
3007  if (yp-5*ep < 0) continue;
3008  if (xp < xpi) xpi = xp;
3009  }
3010  Double_t total = proj->Integral()*proj->GetBinWidth(5);
3011  g2->SetParameters(0, xpi, 0.0, 0.6, 0.1, 0.1, 0.0,0.0,-1., dX, ddX);
3012  g2->FixParameter(5,0);
3013  g2->FixParameter(6,0);
3014  g2->FixParameter(7,total);
3015  g2->FixParameter(8,-1);
3016  g2->FixParameter(9,dX);
3017  g2->FixParameter(10,ddX);
3018  gXFunc();
3019  TFitResultPtr res = proj->Fit(g2,Opt.Data());
3020  Int_t iok = res->Status();
3021  if ( iok < 0) {
3022  cout << g2->GetName() << " fit has failed with " << iok << " for "
3023  << proj->GetName() << "/" << proj->GetTitle() << " Try one again" << endl;
3024  proj->Fit(g2,Opt.Data());
3025  }
3026  // g2->ReleaseParameter(5); g2->SetParLimits(5,0.0,TMath::Pi()/2);
3027  Opt += "m";
3028  gXFunc();
3029  res = proj->Fit(g2,Opt.Data());
3030  iok = res->Status();
3031  if (iok < 0 ) return 0;
3032  if (! Opt.Contains("q",TString::kIgnoreCase)) {
3033  gXFunc();
3034  proj->Draw();
3035  Double_t params[11];
3036  g2->GetParameters(params);
3037  Double_t X = params[1];
3038  Double_t Y = TMath::Exp(params[0]);
3039  TPolyMarker *pm = new TPolyMarker(1, &X, &Y);
3040  proj->GetListOfFunctions()->Add(pm);
3041  pm->SetMarkerStyle(23);
3042  pm->SetMarkerColor(kRed);
3043  pm->SetMarkerSize(1.3);
3044  for (int i = 0; i < NFIT_HYP; i++) {
3045  TF1 *f = new TF1(*g2);
3046  f->SetName(Peaks[i].Name);
3047  f->FixParameter(8,i);
3048  f->SetLineColor(i+2);
3049  gXFunc();
3050  f->Draw("same");
3051  proj->GetListOfFunctions()->Add(f);
3052  }
3053  }
3054  return g2;
3055 }
3056 //________________________________________________________________________________
3057 void dEdxFit() {}
3058 //________________________________________________________________________________
3059 void dEdxFit(const Char_t *HistName,const Char_t *FitName = "GP",
3060  Option_t *opt="R",
3061  Int_t ix = -1, Int_t jy = -1,
3062  Int_t mergeX=1, Int_t mergeY=1,
3063  Double_t nSigma=3, Int_t pow=1) {
3064  TCanvas *canvas = 0;
3065  TString Opt(opt);
3066  if (! Opt.Contains("Q",TString::kIgnoreCase)) {
3067  canvas = (TCanvas *) gROOT->GetListOfCanvases()->FindObject("Fit");
3068  if (! canvas) canvas = new TCanvas("Fit","Fit results");
3069  else canvas->Clear();
3070  }
3071  TList *list = (TList *) gROOT->GetListOfFiles();
3072  if (! list) {printf("File list is empty\n"); return;}
3073  TIter next(list);
3074  TFile *fRootFile = (TFile *) next(); // gDirectory->GetFile();
3075  if (! fRootFile ) {printf("There is no opened file\n"); return;}
3076  TH1 *hist = 0;// = (TH1 *) fRootFile->Get(HistName);
3077  // fRootFile->GetObject(HistName,hist);
3078  TObject *obj = fRootFile->Get(HistName);
3079  if (!obj) {printf("Cannot find %s\n", HistName); return;}
3080  if (obj->IsA()->InheritsFrom( "TH1" )) {
3081  hist = (TH1 *) obj;
3082  } else if (obj->IsA()->InheritsFrom( "THnSparse" )) {
3083  hist = ((THnSparse *) obj)->Projection(1,0);
3084  }
3085  if (! hist) return;
3086  TAxis *xax = hist->GetXaxis();
3087  Int_t nx = xax->GetNbins(); printf ("nx = %i",nx);
3088  Axis_t xmin = xax->GetXmin(); printf (" xmin = %f",xmin);
3089  Axis_t xmax = xax->GetXmax(); printf (" xmax = %f\n",xmax);
3090  TAxis *yax = hist->GetYaxis();
3091  Int_t dim = hist->GetDimension();
3092  Int_t ny = yax->GetNbins();
3093  Int_t NX = nx;
3094  Int_t NY = ny;
3095  if (dim < 3) ny = 1; printf ("ny = %i",ny);
3096  Axis_t ymin = yax->GetXmin(); printf (" ymin = %f",ymin);
3097  Axis_t ymax = yax->GetXmax(); printf (" ymax = %f\n",ymax);
3098  struct Fit_t {
3099  Float_t i;
3100  Float_t j;
3101  Float_t x;
3102  Float_t y;
3103  Float_t mean;
3104  Float_t rms;
3105  Float_t peak;
3106  Float_t mu;
3107  Float_t sigma;
3108  Float_t entries;
3109  Float_t chisq;
3110  Float_t prob;
3111  Float_t a0;
3112  Float_t a1;
3113  Float_t a2;
3114  Float_t a3;
3115  Float_t a4;
3116  Float_t a5;
3117  Float_t Npar;
3118  Float_t dpeak;
3119  Float_t dmu;
3120  Float_t dsigma;
3121  Float_t da0;
3122  Float_t da1;
3123  Float_t da2;
3124  Float_t da3;
3125  Float_t da4;
3126  Float_t da5;
3127 #ifdef __HEED_MODEL__
3128  Float_t Scale;
3129  Float_t dScale;
3130 #endif /* __HEED_MODEL__ */
3131  };
3132  Fit_t Fit;
3133  // TString NewRootFile(gSystem->DirName(fRootFile->GetName()));
3134  TString NewRootFile(gSystem->DirName(gSystem->BaseName(fRootFile->GetName())));
3135  NewRootFile += "/";
3136  NewRootFile += HistName;
3137  NewRootFile += FitName;
3138  if (ix >= 0) NewRootFile += Form("_X%i",ix);
3139  if (jy >= 0) NewRootFile += Form("_Y%i",jy);
3140  if (mergeX != 1) NewRootFile += Form("_x%i",mergeX);
3141  if (mergeY != 1) NewRootFile += Form("_y%i",mergeY);
3142  // NewRootFile += "_2_";
3143  NewRootFile += gSystem->BaseName(fRootFile->GetName());
3144  if (! FitP) {
3145  if (! fOut) {
3146  fOut = new TFile(NewRootFile.Data(),"update");
3147  if (! fOut) fOut = new TFile(NewRootFile.Data(),"new");
3148  if (fOut) cout << NewRootFile << " has been opened." << endl;
3149  else {cout << "Failed to open " << NewRootFile << endl; return;}
3150  }
3151  FitP = (TNtuple *) fOut->Get("FitP");
3152  }
3153  if (! FitP) {
3154  FitP = new TNtuple("FitP","Fit results",
3155 #ifndef __HEED_MODEL__
3156  "i:j:x:y:mean:rms:peak:mu:sigma:entries:chisq:prob:a0:a1:a2:a3:a4:a5:Npar:dpeak:dmu:dsigma:da0:da1:da2:da3:da4:da5");
3157 #else /* __HEED_MODEL__ */
3158  "i:j:x:y:mean:rms:peak:mu:sigma:entries:chisq:prob:a0:a1:a2:a3:a4:a5:Npar:dpeak:dmu:dsigma:da0:da1:da2:da3:da4:da5:Scale:dScale");
3159 #endif /* __HEED_MODEL__ */
3160  FitP->SetMarkerStyle(20);
3161  FitP->SetLineWidth(2);
3162  }
3163  TH1 *mean = (TH1 *) fOut->Get("mean");
3164  TH1 *rms = (TH1 *) fOut->Get("rms");
3165  TH1 *entries = (TH1 *) fOut->Get("entries");
3166  TH1 *mu = (TH1 *) fOut->Get("mu");
3167  TH1 *sigma= (TH1 *) fOut->Get("sigma");
3168  TH1 *chisq= (TH1 *) fOut->Get("chisq");
3169  if (! mu) {
3170  if (dim == 3) {
3171  mean = new TH2D("mean",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3172  rms = new TH2D("rms",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3173  entries = new TH2D("entries",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3174  mu = new TH2D("mu",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3175  sigma = new TH2D("sigma",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3176  chisq = new TH2D("chisq",hist->GetTitle(),nx,xmin,xmax,ny,ymin,ymax);
3177  }
3178  else {
3179  if (dim == 2 || dim == 1) {
3180  mean = new TH1D("mean",hist->GetTitle(),nx,xmin,xmax);
3181  rms = new TH1D("rms",hist->GetTitle(),nx,xmin,xmax);
3182  entries = new TH1D("entries",hist->GetTitle(),nx,xmin,xmax);
3183  mu = new TH1D("mu",hist->GetTitle(),nx,xmin,xmax);
3184  sigma = new TH1D("sigma",hist->GetTitle(),nx,xmin,xmax);
3185  chisq = new TH1D("chisq",hist->GetTitle(),nx,xmin,xmax);
3186  }
3187  else {
3188  printf("Histogram %s has wrong dimension %i\n", hist->GetName(),dim);
3189  return;
3190  }
3191  }
3192  }
3193  Double_t params[20] = {0};
3194  TH1 *proj = 0;
3195  TF1 *g = 0;
3196  Int_t ix1 = ix, jy1 = jy;
3197  if (ix >= 0) nx = ix;
3198  if (jy >= 0) ny = jy;
3199  if (ix1 < 0) ix1 = 0;
3200  if (jy1 < 0) jy1 = 0;
3201  for (int i=ix1;i<=nx-mergeX+1;i++){
3202  Int_t ir0 = i;
3203  Int_t ir1=i+mergeX-1;
3204  if (i == 0) {ir0 = 1; ir1 = nx;}
3205  for (int j=jy1;j<=ny-mergeY+1;j++){
3206  Int_t jr0 = j;
3207  Int_t jr1 = j+mergeY-1;
3208  if (j == 0) {jr0 = 1; jr1 = ny;}
3209  if (dim == 3) {
3210  if (ir0 == ir1 && jr0 == jr1) {
3211  proj = ((TH3 *) hist)->ProjectionZ(Form("f%i_%i", ir0, jr0 ),ir0,ir1,jr0,jr1);
3212  cout<<"Histogram "<<proj->GetName()<<" was created"<<endl;
3213  } else {
3214  proj = ((TH3 *) hist)->ProjectionZ(Form("f%i_%i_%i_%i",ir0,ir1,jr0,jr1),ir0,ir1,jr0,jr1);
3215  cout<<"Histogram "<<proj->GetName()<<" was created"<<endl;
3216  }
3217  if (! proj) continue;
3218  TString title(proj->GetTitle());
3219  title += Form(" in x [%5.1f,%5.1f] and y [%5.1f,%5.1f] range",
3220  xax->GetBinLowEdge(ir0), xax->GetBinUpEdge(ir1),
3221  yax->GetBinLowEdge(jr0), yax->GetBinUpEdge(jr1));
3222  proj->SetTitle(title.Data());
3223  }
3224  else {
3225  if (ir0 == ir1)
3226  { proj = ((TH2 *) hist)->ProjectionY(Form("f%i", ir0), ir0,ir1);cout<<"Histogram "<<proj->GetName()<<" was created"<<endl;}
3227  else
3228  { proj = ((TH2 *) hist)->ProjectionY(Form("f%i_%i",ir0,ir1),ir0,ir1);cout<<"Histogram "<<proj->GetName()<<" was created"<<endl;}
3229  if (! proj) continue;
3230  TString title(proj->GetTitle());
3231  title += Form("in x [%5.1f,%5.1f] range",
3232  xax->GetBinLowEdge(ir0), xax->GetBinUpEdge(ir1));
3233  proj->SetTitle(title.Data());
3234  }
3235  cout << "i/j\t" << i << "/" << j << "\t" << proj->GetName() << "\t" << proj->GetTitle() << "\t" << proj->Integral() << endl;
3236  // continue;
3237  memset (&Fit, 0, sizeof(Fit_t));
3238  Fit.i = (2.*i+mergeX-1.)/2;
3239  Fit.j = (2.*j+mergeY-1.)/2;
3240  Fit.x = 0.5*(xax->GetBinLowEdge(i) + xax->GetBinUpEdge(i+mergeX-1));
3241  Fit.y = 0.5*(yax->GetBinLowEdge(j) + yax->GetBinUpEdge(j+mergeY-1));
3242  Fit.mean = proj->GetMean();
3243  Fit.rms = proj->GetRMS();
3244  Fit.chisq = -100;
3245  Fit.prob = 0;
3246  Fit.entries = proj->Integral();
3247  if (Fit.entries < 100) {delete proj; continue;}
3248  if (TString(FitName) == "GP") g = FitGP(proj,opt,nSigma,pow);
3249  else if (TString(FitName) == "G2") g = FitG2(proj,opt);
3250  else if (TString(FitName) == "NF" && dim == 3) {
3251  TH3 *hists[5] = {0};
3252  static const Char_t *names[5] = {"pi","P","K","e","d"};
3253  for (Int_t ih = 0; ih < 5; ih++) {
3254  hists[ih] = (TH3 *) fRootFile->Get(Form("%s%s",HistName,names[ih]));
3255  SafeDelete(projNs[ih]);
3256  if (hists[ih]) {
3257  if (ir0 == ir1 && jr0 == jr1) {
3258  projNs[ih] = ((TH3 *) hists[ih])->ProjectionZ(Form("f%s%i_%i",names[ih], ir0, jr0 ),ir0,ir1,jr0,jr1);
3259  cout<<"Histogram "<<projNs[ih]->GetName()<<" was created"<<endl;
3260  } else {
3261  projNs[ih] = ((TH3 *) hists[ih])->ProjectionZ(Form("f%s%i_%i_%i_%i",names[ih],ir0,ir1,jr0,jr1),ir0,ir1,jr0,jr1);
3262  cout<<"Histogram "<<projNs[ih]->GetName()<<" was created"<<endl;
3263  }
3264  }
3265  }
3266  Opt = opt;
3267  Opt += "S";
3268 #ifndef __HEED_MODEL__
3269  g = FitNF(proj,Opt);
3270 #else /* __HEED_MODEL__ */
3271  g = FitNF(proj,Opt,Fit.j);
3272 #endif /* __HEED_MODEL__ */
3273  }
3274  else if (TString(FitName) == "XF" && dim == 3) {
3275  Opt = opt;
3276  Opt += "S";
3277  TProfile2D *pdX = (TProfile2D *) fRootFile->Get(Form("%sdX",HistName));
3278  if (pdX) {
3279  Double_t dX = 0, ddX = 0;
3280  Int_t nn = 0;
3281  for (Int_t ii = ir0; ii <= ir1; ii++) {
3282  for (Int_t jj = jr0; jj <= jr1; jj++) {
3283  Double_t d = pdX->GetBinContent(ii,jj);
3284  Double_t e = pdX->GetBinError(ii,jj);
3285  nn++;
3286  dX += d;
3287  ddX += e*e;
3288  }
3289  }
3290  dX = dX/nn;
3291  ddX = TMath::Sqrt(ddX/nn);
3292  g = FitXF(proj,Opt,dX,ddX);
3293  } else {
3294  g = FitXF(proj,Opt);
3295  }
3296  }
3297  else if (TString(FitName) == "GF") g = FitGF(proj,opt);
3298  else if (TString(FitName) == "L5") g = FitL5(proj,opt,5);
3299  else if (TString(FitName) == "L1") g = FitL5(proj,opt,0);
3300 #ifdef __USE_ROOFIT__
3301  else if (TString(FitName) == "R1") g = FitR5(proj,opt,0);
3302  else if (TString(FitName) == "R5") g = FitR5(proj,opt,5);
3303  else if (TString(FitName) == "RL5")
3304  {
3305  Bool_t outer = kFALSE;
3306  if (NX == 45) {outer = i > 13;}
3307  else if (NY == 45) {outer = j > 13;}
3308  g = FitRL5(proj,outer);
3309  // g->Print();
3310  // cout<<"TEST (SELEKCJA) : i = "<<i<<", j = "<<j<<", outer = "<<outer<<endl;
3311  }
3312  else if (TString(FitName) == "RL1") g = FitRL1(proj);
3313 #endif /* __USE_ROOFIT__ */
3314  else if (TString(FitName) == "GB") {
3315  Double_t dX = 2.0; // <dX> Outer
3316  g = FitGB(proj,opt,dX);
3317  }
3318  else {cout << FitName << " has not been definded" << endl; break;}
3319  if ( g ) {
3320  Int_t kpeak = 0;
3321  if (TString(FitName) == "RL5" || TString(FitName) == "RL1") kpeak = 10;
3322  g->GetParameters(params);
3323  Fit.Npar = g->GetNpar();
3324  Fit.chisq = g->GetChisquare();
3325  Fit.prob = g->GetProb();
3326  Fit.peak = params[kpeak]; // norm, Mu for RL5
3327  Fit.mu = params[1];
3328  Fit.sigma = TMath::Abs(params[2]);
3329  Fit.a0 = params[3]; // FitGF "P"
3330  Fit.a1 = params[4]; // "K"
3331  Fit.a2 = params[5]; // "e"
3332  Fit.a3 = params[6]; // "d"
3333  Fit.a4 = params[7]; // "Total"
3334  Fit.a5 = params[8]; // "Case", sigma of Landau for L5
3335  Fit.dpeak = g->GetParError(kpeak);
3336  Fit.dmu = g->GetParError(1);
3337  Fit.dsigma = g->GetParError(2);
3338  Fit.da0 = g->GetParError(3);
3339  Fit.da1 = g->GetParError(4);
3340  Fit.da2 = g->GetParError(5);
3341  Fit.da3 = g->GetParError(6);
3342  Fit.da4 = g->GetParError(7);
3343  Fit.da5 = g->GetParError(8);
3344 #ifdef __HEED_MODEL__
3345  if (Fit.Npar > 9) {
3346  Fit.Scale = params[9];
3347  Fit.dScale = g->GetParError(9);
3348  }
3349 #endif /* __HEED_MODEL__ */
3350  } else {
3351  delete proj; continue;
3352  }
3353  // Fit.chisq = g3->GetChisquare();
3354  // if (Fit.prob > 0) {
3355  if (dim == 3) {
3356  mean->SetBinContent(i,j,Fit.mean);
3357  rms->SetBinContent(i,j,Fit.rms);
3358  entries->SetBinContent(i,j,Fit.entries);
3359  mu->SetBinContent(i,j,Fit.mu);
3360  mu->SetBinError(i,j,g->GetParError(1));
3361  sigma->SetBinContent(i,j,Fit.sigma);
3362  sigma->SetBinError(i,j,g->GetParError(2));
3363  chisq->SetBinContent(i,j,Fit.chisq);
3364  }
3365  else {
3366  mean->SetBinContent(i,Fit.mean);
3367  rms->SetBinContent(i,Fit.rms);
3368  entries->SetBinContent(i,Fit.entries);
3369  mu->SetBinContent(i,Fit.mu);
3370  mu->SetBinError(i,g->GetParError(1));
3371  sigma->SetBinContent(i,Fit.sigma);
3372  sigma->SetBinError(i,g->GetParError(2));
3373  chisq->SetBinContent(i,Fit.chisq);
3374  }
3375  // }
3376  printf("%i/%i %f/%f mean %f rms = %f entries = %f mu = %f sigma = %f chisq = %f prob = %f\n",
3377  i,j,Fit.x,Fit.y,Fit.mean,Fit.rms,Fit.entries,Fit.mu,Fit.sigma,Fit.chisq,Fit.prob);
3378  if (FitP) FitP->Fill(&Fit.i);
3379  if (canvas) {
3380  canvas->Update();
3381  }
3382  fOut->cd();
3383  proj->Write();
3384  }
3385  }
3386  fOut->cd();
3387  FitP->Write();
3388  mean->Write();
3389  rms->Write();
3390  entries->Write();
3391  mu->Write();
3392  sigma->Write();
3393  chisq->Write();
3394 }
peak_t
&lt;peak postion&gt;=&quot;&quot;&gt; at p = [0.45,0.50] GeV/c wrt pion
Definition: dEdxFit.C:110
Generated by   doxygen 1.8.5