StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
PIDFitter.cxx
1 //
3 // $Id: PIDFitter.cxx,v 1.11 2007/07/12 19:47:31 fisyak Exp $
4 //
5 // Authors: Aihong Tang
6 //
8 //
9 // Description: fit histo. produced by StPidAmpMaker.
10 //
12 
13 #include "PIDFitter.h"
14 #include "Stsstream.h"
15 #include "Stiostream.h"
16 #include <math.h>
17 //#include "TH1.h"
18 #include "TH2.h"
19 #include "TLine.h"
20 #include "TF1.h"
21 #include <float.h>
22 #include "TStyle.h"
23 #include "TCanvas.h"
24 #include "TMath.h"
25 #include "TFile.h"
26 #include "TGraphErrors.h"
27 
28 
29 #include "BetheBlochFunction.hh"
30 //#include "MaxllBoltz.hh"
31 
32 
33 #include <Stiostream.h>
34 
35 ClassImp(PIDFitter);
36 Double_t sigmaNSampleFitFcn(Double_t* x, Double_t *par);
37 
38 //----------------------------------------------------------------------
39 PIDFitter::PIDFitter(){
40  Init();
41 }
42 //----------------------------------------------------------------------
43 void PIDFitter::Process( Char_t* sigmaOfSigmTrialInputName,
44  Char_t* sigmaOfSigmTrialOutputName,
45  Char_t* phaseSpaceCalibInputName,
46  Char_t* phaseSpaceCalibOutputName,
47  Char_t* gausFitInputName,
48  Char_t* gausFitOutputName,
49  Char_t* ampFitOutputName ){
50 
51  GetSigmaOfSingleTrail(sigmaOfSigmTrialInputName,sigmaOfSigmTrialOutputName);
52  DoPhaseSpaceCalibration(phaseSpaceCalibInputName,phaseSpaceCalibOutputName);
53  FitMultiGaus(gausFitInputName,gausFitOutputName);
54  ExtrapAmp(gausFitOutputName,ampFitOutputName);
55 
56 }
57 //----------------------------------------------------------------------
58 void PIDFitter::GetSigmaOfSingleTrail(Char_t* fileName4SigmaOfSingleTrail,Char_t* sigmaFileName){//use pion
59 
60  cout<<" getting sigma of single trail..."<<endl;
61  TFile* fitResoFile;
62 
63  if (mWriteSigmaNSampleGraph) {
64  fitResoFile = new TFile("sigmaNSample.root","UPDATE");
65  }
66 
67 
68  ofstream SigmaFileOut;
69  SigmaFileOut.open(sigmaFileName);
70 
71  float pionPosition4Calib=0.4;
72 
73  int i = int((pionPosition4Calib-mPStart)*mNPBins/(mPEnd-mPStart) );
74 
75  TH1F* theHist=0;
76 
77  TFile histo4CalibFile(fileName4SigmaOfSingleTrail,"READ");
78 
79  TH1F* theDummyHisto = (TH1F *)histo4CalibFile.Get("h2000");
80 
81 
82  for (int j=0; j<mNEtaBins; j++) {
83 
84 
85  TGraphErrors* sigmaNSampleGraph = new TGraphErrors();
86 
87  TString grName;
88  grName.Append(fileName4SigmaOfSingleTrail);
89  grName.Append("SigmaNSample");
90  grName.Append((i<10) ? "0" : "");
91  grName +=i;
92  grName +=j;
93  grName.ReplaceAll(".root","");
94 
95  sigmaNSampleGraph->SetTitle(grName.Data());
96  sigmaNSampleGraph->SetName(grName.Data());
97 
98 
99  for (int k=0; k<mNNHitsBins;k++){
100 
101  TH1F* tempHisto = new TH1F(*theDummyHisto);
102  tempHisto->Reset();
103 
104  char *theName = new char[80];
105  if (i<10)
106  sprintf(theName,"h0%d%d%d",i,j,k);
107  else sprintf(theName,"h%d%d%d",i,j,k);
108 
109  theHist= (TH1F *)histo4CalibFile.Get(theName);
110  delete theName;
111 
112 
113  for (int mm=0; mm<tempHisto->GetNbinsX(); mm++)
114  tempHisto->SetBinContent(mm,(theHist->GetBinContent(mm)));
115 
116 
117  Float_t theError=0;
118 
119  Float_t theFitGausHalfRange = 0.03e-05*0.95;
120 
121  Float_t mean=FitResoGaus(tempHisto,theFitGausHalfRange,theError,
122  0.005e-4,0.04e-4,1,j,k,pionPosition4Calib);
123  Float_t sigma=FitResoGaus(tempHisto,theFitGausHalfRange,theError,
124  0.005e-4,0.04e-4,2,j,k,pionPosition4Calib);
125 
126  sigmaNSampleGraph->SetPoint(sigmaNSampleGraph->GetN(),(k*mNNHitsBins+(5./2.)),sigma/mean);
127  sigmaNSampleGraph->SetPointError((sigmaNSampleGraph->GetN()-1),0,theError);
128 
129  if (tempHisto) delete tempHisto;
130 
131  }
132 
133 
134  TF1* f1 = new TF1("f1",sigmaNSampleFitFcn,12.,40.,1);
135  sigmaNSampleGraph->Fit("f1","R");
136 
137  cout<<"sigNSamplePion1Graph fit parameter: "<<f1->GetParameter(0);
138  double theSigmaOfSingleTrail= f1->GetParameter(0);
139  double theErrorOfSigmaOfSingleTrail = f1->GetParError(0);
140  if(theErrorOfSigmaOfSingleTrail) {/*warnOff*/}
141  if (mWriteSigmaNSampleGraph) {
142  fitResoFile->cd();
143  sigmaNSampleGraph->SetMarkerColor(2);
144  sigmaNSampleGraph->SetMarkerStyle(20);
145  sigmaNSampleGraph->Write(grName.Data(),TObject::kOverwrite | TObject::kSingleKey);
146  fitResoFile->Write();
147  }
148 
149 
150  if (f1) delete f1;
151  if (sigmaNSampleGraph) delete sigmaNSampleGraph;
152  mSigmaOfSingleTrail[j]=theSigmaOfSingleTrail;
153 
154 
155 
156  SigmaFileOut<<j<<endl;
157  SigmaFileOut<<theSigmaOfSingleTrail<<endl;//" +/- "<<theErrorOfSigmaOfSingleTrail<<endl;
158 
159  } //integrated over all eta bins
160 
161  if (mWriteSigmaNSampleGraph) fitResoFile->Close();
162  SigmaFileOut.close();
163 
164 
165  cout<<" end of getting sigma of single trail..."<<endl;
166 
167 }
168 
169 //----------------------------------------------------------------------
170 float PIDFitter::FitResoGaus(TH1F* resoHist,float fitRange,float& er,float theStart, float theEnd, int ParIndex, int j, int k, float thePPosition){
171  //get the par from fitting of sigma _ Nsample.
172 
173  TFile* fitResoFile;
174 
175  if (mWriteGaus4SigmaNSampleHist) {
176  fitResoFile = new TFile("Gaus4SigmaNSample.root","UPDATE"); //no name change!
177  }
178 
179  Double_t tempSigmaOfSingleTrial =0.42;
180 
181  RefreshPresettings(resoHist, tempSigmaOfSingleTrial, k, thePPosition);
182 
183  //setup vary range
184 
185  Double_t ECenterVary=0.3;
186  Double_t EWidthVary=0.4;
187 
188  Double_t PiCenterVary=0.3;
189  Double_t PiHeightVary=0.2;
190  Double_t PiWidthVary=0.4;
191 
192  Double_t KCenterVary=0.3;
193  Double_t KWidthVary=0.4;
194 
195  TF1* sumGs
196  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)", 0., 0.08e-4);
197 
198  Double_t pars_b[9] = {mPiGausHeight, mPiGausCenter, mPiSigma,
199  mEGausHeight, mEGausCenter, mESigma,
200  mKGausHeight, mKGausCenter, mKSigma };
201 
202  sumGs->SetParameters(&pars_b[0]);
203 
204  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
205  mPiGausHeight*(1.+PiHeightVary) );
206  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
207  mPiGausCenter*(1.+PiCenterVary) );
208  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
209  mPiSigma*(1.+PiWidthVary) );
210 
211  sumGs->SetParLimits( 3, 0.,mEGausHeight);
212  sumGs->SetParLimits( 4, mEGausCenter*(1.-ECenterVary),
213  mEGausCenter*(1.+ECenterVary) );
214  sumGs->SetParLimits( 5, mESigma*(1.-EWidthVary),
215  mESigma*(1.+EWidthVary) );
216 
217  sumGs->SetParLimits( 6, 0.,mKGausHeight);
218  sumGs->SetParLimits( 7, mKGausCenter*(1.-KCenterVary),
219  mKGausCenter*(1.+KCenterVary) );
220  sumGs->SetParLimits( 8, mKSigma*(1.-KWidthVary),
221  mKSigma*(1.+KWidthVary) );
222 
223 
224  resoHist->Fit("sumGs","R");
225 
226  if (mWriteGaus4SigmaNSampleHist){
227  fitResoFile->cd();
228  resoHist->Write();
229  }
230 
231  if (mWriteGaus4SigmaNSampleHist) fitResoFile->Close();
232  float thePar=sumGs->GetParameter(ParIndex);
233  if (mWriteGaus4SigmaNSampleHist) { delete fitResoFile; fitResoFile=0;}
234  if (sumGs) { delete sumGs; }
235 
236  return thePar;
237 }
238 
239 //----------------------------------------------------------------------
240 void PIDFitter::RefreshPresettings(TH1F* resoHist, double tempSigmaOfSingleTrial, int k, float thePPosition){
241 
242  mEGausHeight =0; mESigma=0.;
243  mPiGausHeight =0; mPiSigma=0.;
244  mKGausHeight =0; mKSigma=0.;
245  mPGausHeight =0; mPSigma=0.;
246 
247  mEGausCenter = electronBandCenter ->Eval(thePPosition,0,0);
248  mPiGausCenter = pionBandCenter ->Eval(thePPosition,0,0);
249  mKGausCenter = kaonBandCenter ->Eval(thePPosition,0,0);
250  mPGausCenter = antiprotonBandCenter->Eval(thePPosition,0,0);
251 
252  PresetHeightAndSigma(mEGausCenter, mEGausHeight, mESigma,
253  resoHist, tempSigmaOfSingleTrial, k);
254  PresetHeightAndSigma(mPiGausCenter, mPiGausHeight, mPiSigma,
255  resoHist, tempSigmaOfSingleTrial, k);
256  PresetHeightAndSigma(mKGausCenter, mKGausHeight, mKSigma,
257  resoHist, tempSigmaOfSingleTrial, k);
258  PresetHeightAndSigma(mPGausCenter, mPGausHeight, mPSigma,
259  resoHist, tempSigmaOfSingleTrial, k);
260 
261 }
262 
263 //----------------------------------------------------------------------
264 void PIDFitter::PresetHeightAndSigma(double center, double& height, double& sigma, TH1F* resoHist, double tempSigmaOfSingleTrial, int k){
265 
266  int nbins = resoHist->GetNbinsX();
267  int centBin =
268  int((center-mDedxStart)*float(nbins)/((mDedxEnd-mDedxStart)));
269 
270  height = (centBin<nbins) ? resoHist->GetBinContent(centBin) : 0.;
271  sigma = tempSigmaOfSingleTrial*center/TMath::Sqrt((k+0.5)*((double(mNNHitsEnd)-double(mNNHitsStart))/double(mNNHitsBins)));
272 
273 }
274 
275 
276 //----------------------------------------------------------------------
277 void PIDFitter::DoPhaseSpaceCalibration(Char_t* fileName4Calibration,Char_t* phaseSpaceCalibFileName){
278 
279  cout<<" doing phase space calibration.......... "<<endl;
280 
281  double protonRef[mNEtaBins][mNNHitsBins];
282  double pionRef[mNEtaBins][mNNHitsBins];
283 
284  double fitEnd = mDedxEnd;
285  double fitStart = mDedxStart;
286 
287  TFile histo4CalibFile(fileName4Calibration,"READ");
288 
289  TH1F* theHist=0;
290  TF1* sumGs=0;
291 
292  int pionPosition4Calib =int((0.5-mPStart)*mNPBins/(mPEnd-mPStart) );
293  int protonPosition4Calib =int((0.7-mPStart)*mNPBins/(mPEnd-mPStart) );
294 
295 
296  //get protonRef
297  for (int j=0; j<mNEtaBins; j++)
298  for (int k=2; k<mNNHitsBins;k++){
299 
300  int i = protonPosition4Calib; //for protonRef
301 
302  //if dE/dx does not shift regarding to nhits,
303  //set kname=4.(just use long tracks)
304  int kname = k;
305  int jname = j;
306 
307  char *theName = new char[80];
308  if (i<10)
309  sprintf(theName,"h0%d%d%d",i,jname,kname);
310  else sprintf(theName,"h%d%d%d",i,jname,kname);
311 
312  theHist= (TH1F *)histo4CalibFile.Get(theName);
313  delete theName;
314 
315  float pPosition =
316  float(i)*((mPEnd-mPStart)/float(mNPBins)) +
317  (((mPEnd-mPStart)/mNPBins)/2.);
318 
319  //setup initial values.
320  RefreshPresettings(theHist, mSigmaOfSingleTrail[j] , k, pPosition);
321 
322  //setup vary range
323  Double_t KCenterVary=0.15;
324  Double_t KWidthVary=0.03;
325 
326  Double_t PiCenterVary=0.1;
327  Double_t PiHeightVary=0.1;
328  Double_t PiWidthVary=0.05;
329 
330  Double_t PCenterVary=0.1;
331  Double_t PWidthVary=0.03;
332 
333 
334  sumGs
335  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)", fitStart, fitEnd);
336 
337  Double_t pars_c[9] = { mPiGausHeight, mPiGausCenter, mPiSigma,
338  mKGausHeight, mKGausCenter, mKSigma,
339  mPGausHeight, mPGausCenter, mPSigma };
340 
341  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
342  mPiGausHeight*(1.+PiHeightVary) );
343  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
344  mPiGausCenter*(1.+PiCenterVary) );
345  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
346  mPiSigma*(1.+PiWidthVary) );
347 
348 
349  sumGs->SetParLimits( 3, 0.,mKGausHeight);
350  sumGs->SetParLimits( 4, mKGausCenter*(1.-KCenterVary),
351  mKGausCenter*(1.+KCenterVary) );
352  sumGs->SetParLimits( 5, mKSigma*(1.-KWidthVary),
353  mKSigma*(1.+KWidthVary) );
354 
355 
356  sumGs->SetParLimits( 6, 0.,mPGausHeight);
357  sumGs->SetParLimits( 7, mPGausCenter*(1.-PCenterVary),
358  mPGausCenter*(1.+PCenterVary) );
359  sumGs->SetParLimits( 8, mPSigma*(1.-PWidthVary),
360  mPSigma*(1.+PWidthVary) );
361 
362  sumGs->SetParameters(&pars_c[0]);
363 
364  theHist->Fit("sumGs","NR");
365 
366  protonRef[j][k]=sumGs->GetParameter(7);
367 
368  if (sumGs) {delete sumGs;}
369  }
370 
371 
372  //get pionRef
373  for (int j=0; j<mNEtaBins; j++)
374  for (int k=2; k<mNNHitsBins;k++){
375 
376  int i = pionPosition4Calib; //for pionRef
377 
378  //if dE/dx does not shift regarding to nhits,
379  //set kname=4.(just use long tracks)
380  int kname = k;
381  int jname = j;
382 
383  char *theName = new char[80];
384  if (i<10)
385  sprintf(theName,"h0%d%d%d",i,jname,kname);
386  else sprintf(theName,"h%d%d%d",i,jname,kname);
387 
388  theHist= (TH1F *)histo4CalibFile.Get(theName);
389 
390  float pPosition =
391  float(i)*((mPEnd-mPStart)/float(mNPBins)) +
392  (((mPEnd-mPStart)/mNPBins)/2.);
393 
394  //setup initial values.
395  RefreshPresettings(theHist, mSigmaOfSingleTrail[j] , k, pPosition);
396 
397  //setup vary range
398  Double_t KCenterVary=0.15;
399  Double_t KWidthVary=0.03;
400 
401  Double_t PiCenterVary=0.1;
402  Double_t PiHeightVary=0.1;
403  Double_t PiWidthVary=0.05;
404 
405  Double_t PCenterVary=0.1;
406  Double_t PWidthVary=0.03;
407 
408  sumGs
409  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)", fitStart, fitEnd);
410 
411  Double_t pars_c[9] = { mPiGausHeight, mPiGausCenter, mPiSigma,
412  mKGausHeight, mKGausCenter, mKSigma,
413  mPGausHeight, mPGausCenter, mPSigma };
414 
415  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
416  mPiGausHeight*(1.+PiHeightVary) );
417  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
418  mPiGausCenter*(1.+PiCenterVary) );
419  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
420  mPiSigma*(1.+PiWidthVary) );
421 
422 
423  sumGs->SetParLimits( 3, 0.,mKGausHeight );
424  sumGs->SetParLimits( 4, mKGausCenter*(1.-KCenterVary),
425  mKGausCenter*(1.+KCenterVary) );
426  sumGs->SetParLimits( 5, mKSigma*(1.-KWidthVary),
427  mKSigma*(1.+KWidthVary) );
428 
429 
430  sumGs->SetParLimits( 6, 0.,mPGausHeight );
431  sumGs->SetParLimits( 7, mPGausCenter*(1.-PCenterVary),
432  mPGausCenter*(1.+PCenterVary) );
433  sumGs->SetParLimits( 8, mPSigma*(1.-PWidthVary),
434  mPSigma*(1.+PWidthVary) );
435 
436  sumGs->SetParameters(&pars_c[0]);
437 
438  theHist->Fit("sumGs","NR");
439 
440  pionRef[j][k]=sumGs->GetParameter(1);
441 
442 
443  if (sumGs) {delete sumGs; }
444  }
445 
446  for (int j=0; j<mNEtaBins; j++)
447  for (int k=2; k<mNNHitsBins;k++){
448  double theCalib=4.68971e-07;
449  double theOffSet=3.01022e-07;
450  double deltaReference=protonRef[j][k]-pionRef[j][k];
451  double mimimumIonizingPionPosition=
452  float(pionPosition4Calib)*((mPEnd-mPStart)/float(mNPBins)) +
453  (((mPEnd-mPStart)/mNPBins)/2.);
454  double protonTestPosition=
455  float(protonPosition4Calib)*((mPEnd-mPStart)/float(mNPBins)) +
456  (((mPEnd-mPStart)/mNPBins)/2.);
457 
458  BBScalePar[j][k]=look4MinDeltaDiff(theCalib*0.7, theCalib*1.3, 100,mimimumIonizingPionPosition ,protonTestPosition , deltaReference);
459  BBOffSetPar[j][k]= theOffSet+minimumIonizingdEdx(BBScalePar[j][k],mimimumIonizingPionPosition )-pionRef[j][k];
460 
461  cout<<" j "<<j<<" k "<<k<<" BBScalePar["<<j<<"]["<<k<<"]="<<BBScalePar[j][k]<<"; BBOffSetPar["<<j<<"]["<<k<<"]="<<BBOffSetPar[j][k]<<"; "<<endl;
462 
463  }
464 
465  //now deal with k<2
466  for (int j=0; j<mNEtaBins; j++)
467  for (int k=0; k<2; k++) {
468  BBScalePar[j][k]=BBScalePar[j][2];
469  BBOffSetPar[j][k]=BBOffSetPar[j][2];
470  } //nhits<10, hard to calibrate. use nhits 10-15 instead.
471 
472 
473  ofstream PhaseSpaceFileOut;
474  PhaseSpaceFileOut.open(phaseSpaceCalibFileName);
475 
476  for (int j=0; j<mNEtaBins; j++)
477  for (int k=0; k<mNNHitsBins; k++) {
478  PhaseSpaceFileOut<<j<<endl<<k<<endl;
479  PhaseSpaceFileOut<<electronBandCenter->GetParameter(0)<<endl;
480  PhaseSpaceFileOut<<electronBandCenter->GetParameter(1)<<endl;
481  PhaseSpaceFileOut<<BBOffSetPar[j][k]<<endl;
482  PhaseSpaceFileOut<<BBScalePar[j][k]<<endl;
483  PhaseSpaceFileOut<<electronBandCenter->GetParameter(6)<<endl;
484 
485  }
486 
487  PhaseSpaceFileOut.close();
488 
489  cout<<" phase space calibration finished .......... "<<endl;
490 
491 }
492 
493 
494 //--------------------------------------------------------------------------
495 void PIDFitter::Init(){
496 
497  mWriteGaus4SigmaNSampleHist=kFALSE;
498 
499  //create a prototype for update in FitResoGaus()
500  if (mWriteGaus4SigmaNSampleHist){
501  TFile* fitResoFile = new TFile("Gaus4SigmaNSample.root","RECREATE");
502  fitResoFile->Write();
503  fitResoFile->Close();
504  delete fitResoFile;
505  }
506 
507  mWriteSigmaNSampleGraph=kFALSE;
508 
509  if (mWriteSigmaNSampleGraph) {
510  TFile* fitResoFile = new TFile("sigmaNSample.root","RECREATE");
511  delete fitResoFile;
512  }
513 
514 
515 
516  mSigmaOfSingleTrail = new double[mNEtaBins];
517 
518  BBOffSetPar = new double*[mNEtaBins];
519  BBScalePar = new double*[mNEtaBins];
520 
521  for (int ii=0; ii<mNEtaBins; ii++) {
522  BBOffSetPar[ii] = new double[mNNHitsBins];
523  BBScalePar[ii] = new double[mNNHitsBins];
524  }
525 
526 
527 
528  //allocate functions...
529  electronBandCenter
530  =new TF1("electronBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
531  pionBandCenter
532  =new TF1("pionBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
533  kaonBandCenter
534  =new TF1("kaonBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
535  antiprotonBandCenter
536  =new TF1("antiprotonBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
537 
538  pionKaonBandCenter //for drawing line between pion and kaon bands
539  =new TF1("pionKaonBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
540  kaonAntiprotonBandCenter //for drawing line between kaon and antiproton bands
541  =new TF1("kaonAntiprotonBandCenter",BetheBlochFunction, mPStart,mPEnd, NParameters);
542 
543 
544 
545 
546 
547  Double_t electronPars[7]
548  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.511e-3, 4.68971e-07, 0.0005 };
549  Double_t pionPars[7]
550  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.13957, 4.68971e-07, 0.0005 };
551  Double_t kaonPars[7]
552  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.49368, 4.68971e-07, 0.0005 };
553  Double_t antiprotonPars[7]
554  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.93827, 4.68971e-07, 0.0005 };
555 
556  Double_t pionKaonPars[7]
557  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.387447*1.1, 4.68971e-07, 0.0005 };
558 
559  Double_t kaonAntiprotonPars[7]
560  ={ 1.09344, 0.0199, 3.01022e-07, 1, 0.804893*1.05, 4.68971e-07, 0.0005 };
561 
562  electronBandCenter->SetParameters(&electronPars[0]);
563  pionBandCenter->SetParameters(&pionPars[0]);
564  kaonBandCenter->SetParameters(&kaonPars[0]);
565  antiprotonBandCenter->SetParameters(&antiprotonPars[0]);
566  pionKaonBandCenter->SetParameters(&pionKaonPars[0]);
567  kaonAntiprotonBandCenter->SetParameters(&kaonAntiprotonPars[0]);
568 
569 
570 
571 }
572 
573 
574 //--------------------------------------------------------------------------
575 double PIDFitter::delta(double calib, double pionPosition, double protonPosition){ //for phase space calibration use only.
576 
577  pionBandCenter->SetParameter(5,calib);
578  antiprotonBandCenter->SetParameter(5,calib);
579 
580  return (antiprotonBandCenter->Eval(protonPosition,0,0)-
581  pionBandCenter->Eval(pionPosition,0,0));
582 }
583 
584 //--------------------------------------------------------------------------
585 double PIDFitter::look4MinDeltaDiff(double calibStart, double calibEnd, int calibSteps, double pionPosition, double protonPosition, double DeltaRef){
586 
587  double calibSeg=(calibEnd-calibStart)/double(calibSteps);
588  double thisCalib=calibStart;
589  double minDeltaDiffCalib=5000;//calib associated with minDeltaDiff.
590  double minDeltaDiff=5000;
591 
592  do {
593 
594  double myDelta=delta(thisCalib,pionPosition,protonPosition);
595  double diff=TMath::Abs(myDelta-DeltaRef);
596  if (diff<minDeltaDiff) {
597  minDeltaDiff=diff;
598  minDeltaDiffCalib=thisCalib;
599  }
600 
601  thisCalib=thisCalib+calibSeg;
602  }while(thisCalib<calibEnd);
603 
604  return minDeltaDiffCalib;
605 }
606 
607 //--------------------------------------------------------------------------
608 double PIDFitter::minimumIonizingdEdx(double calib, double pionPosition){
609  pionBandCenter->SetParameter(5,calib);
610 
611  return pionBandCenter->Eval(pionPosition,0,0);
612 }
613 
614 
615 //--------------------------------------------------------------------------
616 void PIDFitter::FitMultiGaus(Char_t* fileNameOfInput, Char_t* fileNameOfOutput){
617  //ATTN: DoPhaseSpaceCalibration() and GetSigmaOfSingleTrail have to be called
618  //befor calling this function, otherwise it brokes.
619 
620  TFile histoFile(fileNameOfInput,"READ");
621 
622  TFile fittedHistoFile(fileNameOfOutput,"RECREATE");
623 
624 
625  TH1F* theHist =0;
626  TF1* sumGs =0;
627 
628  for (int i=0; i<mNPBins; i++)
629  for (int j=0; j<mNEtaBins; j++)
630  for (int k=0; k<mNNHitsBins;k++){
631 
632  cout<<" p bin # "<<i<<", eta bin # "<<j<<", ndedx bin # "<<k<<endl;
633 
634  char *theName = new char[80];
635  if (i<10)
636  sprintf(theName,"h0%d%d%d",i,j,k);
637  else sprintf(theName,"h%d%d%d",i,j,k);
638 
639  theHist= (TH1F *)histoFile.Get(theName);
640 
641  float pPosition =
642  float(i)*((mPEnd-mPStart)/float(mNPBins)) +
643  (((mPEnd-mPStart)/mNPBins)/2.);
644 
645  //setup initial values.
646 
647  electronBandCenter->SetParameter(5,BBScalePar[j][k]);
648  pionBandCenter->SetParameter(5,BBScalePar[j][k]);
649  kaonBandCenter->SetParameter(5,BBScalePar[j][k]);
650  antiprotonBandCenter->SetParameter(5,BBScalePar[j][k]);
651  pionKaonBandCenter->SetParameter(5,BBScalePar[j][k]);
652  kaonAntiprotonBandCenter->SetParameter(5,BBScalePar[j][k]);
653 
654  electronBandCenter->SetParameter(2,BBOffSetPar[j][k]);
655  pionBandCenter->SetParameter(2,BBOffSetPar[j][k]);
656  kaonBandCenter->SetParameter(2,BBOffSetPar[j][k]);
657  antiprotonBandCenter->SetParameter(2,BBOffSetPar[j][k]);
658  pionKaonBandCenter->SetParameter(2,BBOffSetPar[j][k]);
659  kaonAntiprotonBandCenter->SetParameter(2,BBOffSetPar[j][k]);
660 
661  //setup initial values.
662  RefreshPresettings(theHist, mSigmaOfSingleTrail[j] , k, pPosition);
663 
664  //setup vary range
665  Double_t KCenterVary=0.1;
666  Double_t KWidthVary=0.1;
667 
668  Double_t PiCenterVary=0.1;
669  Double_t PiHeightVary=0.1;
670  Double_t PiWidthVary=0.1;
671 
672  Double_t PCenterVary=0.1;
673  Double_t PWidthVary=0.1;
674 
675  Double_t ECenterVary=0.;
676  Double_t EWidthVary=0.1;
677 
678  //for low p range, the width for pion and e is big, allow them vary for a better fitting.
679 
680  if (pPosition<0.13) {
681  PiWidthVary=0.5;
682  PiCenterVary=0.2;
683  EWidthVary=0.5;
684  ECenterVary=0.2;
685  }
686 
687 
689  // p<0.2 g=pion+e //
690  // 0.2 <p<0.44 g=pion+e+kaon //
691  // 0.44<p<0.85 g=pion+e+kaon+antiproton //
692  // p>0.85 g=pion+kaon+antiproton //
694 
695  Double_t fitStart=0;
696  Double_t fitEnd=0;
697 
698  if (pPosition<junction1){//p<0.2 g=pion+e
699 
700  Double_t kaonBandLowEdge=pionKaonBandCenter->Eval(pPosition,0.,0.);
701  fitEnd = TMath::Min(kaonBandLowEdge,mDedxEnd);
702 
703  sumGs
704  = new TF1("sumGs","gaus(0)+gaus(3)", fitStart, fitEnd);
705 
706  Double_t pars_a[6] = { mPiGausHeight, mPiGausCenter, mPiSigma,
707  mEGausHeight, mEGausCenter, mESigma };
708 
709  sumGs->SetParameters(&pars_a[0]);
710 
711  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
712  mPiGausHeight*(1.+PiHeightVary) );
713  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
714  mPiGausCenter*(1.+PiCenterVary) );
715  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
716  mPiSigma*(1.+PiWidthVary) );
717 
718  sumGs->SetParLimits( 3, 0., mEGausHeight );
719  sumGs->SetParLimits( 4, mEGausCenter*(1.-ECenterVary),
720  mEGausCenter*(1.+ECenterVary));
721  sumGs->SetParLimits( 5, mESigma*(1.-EWidthVary),
722  mESigma*(1.+EWidthVary) );
723 
724 
725  if (theHist->GetEntries()>1.) theHist->Fit("sumGs","R");
726 
727  fittedHistoFile.cd();
728  theHist->Write();
729  if (sumGs) delete sumGs;
730 
731  }
732 
733 
734 
735  if (pPosition>=junction1 && pPosition<junction2){//0.2<p<0.44 g=pion+e+kaon
736 
737  Double_t antiprotonBandLowEdge=kaonAntiprotonBandCenter->Eval(pPosition,0.,0.);
738  fitEnd = TMath::Min(antiprotonBandLowEdge,mDedxEnd);
739 
740  sumGs
741  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)", fitStart, fitEnd);
742 
743  Double_t pars_b[9] = { mPiGausHeight, mPiGausCenter, mPiSigma,
744  mEGausHeight, mEGausCenter, mESigma,
745  mKGausHeight, mKGausCenter, mKSigma };
746 
747  sumGs->SetParameters(&pars_b[0]);
748 
749 
750  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
751  mPiGausHeight*(1.+PiHeightVary));
752  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
753  mPiGausCenter*(1.+PiCenterVary) );
754  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
755  mPiSigma*(1.+PiWidthVary) );
756 
757  sumGs->SetParLimits( 3, 0., mEGausHeight );
758  sumGs->SetParLimits( 4, mEGausCenter*(1.-ECenterVary),
759  mEGausCenter*(1.+ECenterVary) );
760  sumGs->SetParLimits( 5, mESigma*(1.-EWidthVary),
761  mESigma*(1.+EWidthVary) );
762 
763  sumGs->SetParLimits( 6, 0., mKGausHeight );
764  sumGs->SetParLimits( 7, mKGausCenter*(1.-KCenterVary),
765  mKGausCenter*(1.+KCenterVary) );
766  sumGs->SetParLimits( 8, mKSigma*(1.-KWidthVary),
767  mKSigma*(1.+KWidthVary) );
768 
769 
770  theHist->Fit("sumGs","R");
771 
772  fittedHistoFile.cd();
773  theHist->Write();
774  if (sumGs) delete sumGs;
775  }
776 
777 
778 
779  if (pPosition>=junction2 && pPosition<junction3){//0.44<p<0.85 g=pion+e+kaon+antiproton
780 
781  fitEnd = mDedxEnd;
782 
783  sumGs
784  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)+gaus(9)", fitStart, fitEnd);
785 
786  Double_t pars_c[12] = { mPiGausHeight, mPiGausCenter, mPiSigma,
787  mEGausHeight, mEGausCenter, mESigma,
788  mKGausHeight, mKGausCenter, mKSigma,
789  mPGausHeight, mPGausCenter, mPSigma };
790 
791  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
792  mPiGausHeight*(1.+PiHeightVary) );
793  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
794  mPiGausCenter*(1.+PiCenterVary) );
795  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
796  mPiSigma*(1.+PiWidthVary) );
797 
798 
799  sumGs->SetParLimits( 3, 0., mEGausHeight );
800  sumGs->SetParLimits( 4, mEGausCenter*(1.-ECenterVary),
801  mEGausCenter*(1.+ECenterVary) );
802  sumGs->SetParLimits( 5, mESigma*(1.-EWidthVary),
803  mESigma*(1.+EWidthVary) );
804 
805 
806  sumGs->SetParLimits( 6, 0., mKGausHeight );
807  sumGs->SetParLimits( 7, mKGausCenter*(1.-KCenterVary),
808  mKGausCenter*(1.+KCenterVary) );
809  sumGs->SetParLimits( 8, mKSigma*(1.-KWidthVary),
810  mKSigma*(1.+KWidthVary) );
811 
812 
813  sumGs->SetParLimits( 9, 0., mPGausHeight );
814  sumGs->SetParLimits(10, mPGausCenter*(1.-PCenterVary),
815  mPGausCenter*(1.+PCenterVary) );
816  sumGs->SetParLimits(11, mPSigma*(1.-PWidthVary),
817  mPSigma*(1.+PWidthVary) );
818 
819  sumGs->SetParameters(&pars_c[0]);
820 
821  theHist->Fit("sumGs","R");
822 
823  fittedHistoFile.cd();
824  theHist->Write();
825  if (sumGs) delete sumGs;
826  }
827 
828 
829  if (pPosition>=junction3 && pPosition<2.){//p>0.85 g=pion+kaon+antiproton
830 
831  fitEnd = mDedxEnd;
832 
833  sumGs
834  = new TF1("sumGs","gaus(0)+gaus(3)+gaus(6)", fitStart, fitEnd);
835 
836  Double_t pars_d[9] = { mPiGausHeight, mPiGausCenter, mPiSigma,
837  mKGausHeight, mKGausCenter, mKSigma,
838  mPGausHeight, mPGausCenter, mPSigma };
839 
840  sumGs->SetParLimits( 0, mPiGausHeight*(1.-PiHeightVary),
841  mPiGausHeight*(1.+PiHeightVary));
842  sumGs->SetParLimits( 1, mPiGausCenter*(1.-PiCenterVary),
843  mPiGausCenter*(1.+PiCenterVary));
844  sumGs->SetParLimits( 2, mPiSigma*(1.-PiWidthVary),
845  mPiSigma*(1.+PiWidthVary));
846 
847 
848 
849  sumGs->SetParLimits( 3, 0., mKGausHeight );
850  sumGs->SetParLimits( 4, mKGausCenter*(1.-KCenterVary),
851  mKGausCenter*(1.+KCenterVary));
852  sumGs->SetParLimits( 5, mKSigma*(1.-KWidthVary),
853  mKSigma*(1.+KWidthVary));
854 
855 
856  sumGs->SetParLimits( 6, 0., mPGausHeight );
857  sumGs->SetParLimits( 7, mPGausCenter*(1.-PCenterVary),
858  mPGausCenter*(1.+PCenterVary));
859  sumGs->SetParLimits( 8, mPSigma*(1.-PWidthVary),
860  mPSigma*(1.+PWidthVary));
861 
862  sumGs->SetParameters(&pars_d[0]);
863 
864  theHist->Fit("sumGs","R");
865 
866  fittedHistoFile.cd();
867  theHist->Write();
868  if (sumGs) delete sumGs;
869  }
870 
871 
872  }
873 
874 
875  fittedHistoFile.Write();
876  fittedHistoFile.Close();
877 }
878 
879 //---------------------------------------------------------------------------
880 void PIDFitter::ExtrapAmp(Char_t* fileNameOfInput, Char_t* fileNameOfOutput){
881 
882  cout<<" **********extraping pion amp. begin********* "<<endl;
883 
884  Double_t pBegin = mPStart;
885  Double_t pEnd = mPEnd;
886 
887  TFile* gausHistFile = new TFile(fileNameOfInput,"READ");
888 
889  TFile* ampHistFile = new TFile(fileNameOfOutput,"RECREATE");
890 
891 
892  for (int j=0; j<mNEtaBins; j++)
893  for (int k=0; k<mNNHitsBins;k++){
894 
895  char *PiAmpName = new char[80];
896  sprintf(PiAmpName,"PiAmp%d%d",j,k);
897  char *EAmpName = new char[80];
898  sprintf(EAmpName,"EAmp%d%d",j,k);
899  char *KAmpName = new char[80];
900  sprintf(KAmpName,"KAmp%d%d",j,k);
901  char *PAmpName = new char[80];
902  sprintf(PAmpName,"PAmp%d%d",j,k);
903 
904  TH1F* PiAmpHisto = new TH1F(PiAmpName, PiAmpName, mNPBins,mPStart,mPEnd);
905  TH1F* EAmpHisto = new TH1F(EAmpName, EAmpName, mNPBins,mPStart,mPEnd);
906  TH1F* KAmpHisto = new TH1F(KAmpName, KAmpName, mNPBins,mPStart,mPEnd);
907  TH1F* PAmpHisto = new TH1F(PAmpName, PAmpName, mNPBins,mPStart,mPEnd);
908 
909  PiAmpHisto->SetLineColor(4);
910  EAmpHisto->SetLineColor(6);
911  KAmpHisto->SetLineColor(2);
912  PAmpHisto->SetLineColor(3);
913 
914 
915 
916  double nhitsPosition=double((k+0.5)*(float(mNNHitsEnd-mNNHitsStart)/float(mNNHitsBins)));
917 
918 
919  for (int i=0; i<mNPBins; i++){
920 
921  float pPosition =
922  float(i)*((pEnd-pBegin)/float(mNPBins)) +
923  (((pEnd-pBegin)/mNPBins)/2.);
924 
925  char *theHistName = new char[80];
926  if (i<10)
927  sprintf(theHistName,"h0%d%d%d",i,j,k);
928  else sprintf(theHistName,"h%d%d%d",i,j,k);
929 
930  TH1F* theHist= (TH1F *)gausHistFile->Get(theHistName);
931 
932 
933  TF1* sumGs=theHist->GetFunction("sumGs");
934 
935  if (!sumGs) continue; //due to bin 0-3 might not be fitted (no entries)
936 
937  if (pPosition<junction1){//p<0.2 g=pion+e
938  PiAmpHisto->SetBinContent(i+1,sumGs->GetParameter(0));
939  PiAmpHisto->SetBinError(i+1, sumGs->GetParError(0));
940  EAmpHisto->SetBinContent(i+1,sumGs->GetParameter(3));
941  EAmpHisto->SetBinError(i+1, sumGs->GetParError(3));
942  }
943 
944  if (pPosition>=junction1 && pPosition<junction2){//0.2<p<0.38 g=pion+e+kaon
945  PiAmpHisto->SetBinContent(i+1,sumGs->GetParameter(0));
946  PiAmpHisto->SetBinError(i+1, sumGs->GetParError(0));
947  EAmpHisto->SetBinContent(i+1,sumGs->GetParameter(3));
948  EAmpHisto->SetBinError(i+1, sumGs->GetParError(3));
949  KAmpHisto->SetBinContent(i+1,sumGs->GetParameter(6));
950  KAmpHisto->SetBinError(i+1, sumGs->GetParError(6));
951  }
952  if (pPosition>=junction2 && pPosition<junction3 ){//0.85>p>0.38 g=pion+e+kaon+antiproton
953  PiAmpHisto->SetBinContent(i+1,sumGs->GetParameter(0));
954  PiAmpHisto->SetBinError(i+1, sumGs->GetParError(0));
955  EAmpHisto->SetBinContent(i+1,sumGs->GetParameter(3));
956  EAmpHisto->SetBinError(i+1, sumGs->GetParError(3));
957  KAmpHisto->SetBinContent(i+1,sumGs->GetParameter(6));
958  KAmpHisto->SetBinError(i+1, sumGs->GetParError(6));
959  PAmpHisto->SetBinContent(i+1,sumGs->GetParameter(9));
960  PAmpHisto->SetBinError(i+1, sumGs->GetParError(9));
961  }
962 
963  if (pPosition>=junction3 && pPosition<2.){//p>0.85 g=pion+kaon+antiproton
964  PiAmpHisto->SetBinContent(i+1,sumGs->GetParameter(0));
965  PiAmpHisto->SetBinError(i+1, sumGs->GetParError(0));
966  KAmpHisto->SetBinContent(i+1,sumGs->GetParameter(3));
967  KAmpHisto->SetBinError(i+1, sumGs->GetParError(3));
968  PAmpHisto->SetBinContent(i+1,sumGs->GetParameter(6));
969  PAmpHisto->SetBinError(i+1, sumGs->GetParError(6));
970  }
971 
972 
973  }
974 
975  //fit Pi Amp
976  double low = 0.15;
977  double high = 0.6;
978  double ampPeakPos = 0.45;
979 
980 
981  //find the peak position of Pi amp.
982  for (int t=0; t<4; t++){
983  TF1* gs = new TF1("gs","gaus",low,high);
984  PiAmpHisto->Fit("gs","NWR");
985  low = gs->GetParameter(1)-0.13;
986  high = gs->GetParameter(1)+0.13;
987  ampPeakPos= gs->GetParameter(1);
988  delete gs;
989  }
990 
991  // if (nhitsPosition <= 15.) {//ampPeakPos found fails for low ndedx, assign a number.
992  ampPeakPos=0.6;
993  // }
994 
995 
996  TF1* PiFcnCenter = new TF1("PiFcnCenter","pol3",0.4,0.7);
997 
998  TF1* PiFcnLeft = new TF1("PiFcnLeft","pol3",0.1,ampPeakPos-0.05);
999  PiFcnLeft->SetLineColor(4);
1000 
1001  TF1* PiFcnRight = new TF1("PiFcnRight","expo",0.7,1.4);
1002  PiFcnRight->SetLineColor(3);
1003 
1004  PiAmpHisto->Fit("PiFcnCenter","WR+");
1005  // PiAmpHisto->Fit("PiFcnLeft", "WR+");
1006 
1007  // blow error between 0.85-1.2 GeV/c, because K pi crosssing, put less weight on fitting for that range.
1008 
1009  float crossBandBegin = 0.8;
1010  float brossBandEnd = 1.2;
1011 
1012 
1013  TGraphErrors* PiAmpGraphRight = new TGraphErrors();
1014 
1015  for (int ampbin = 0; ampbin<PiAmpHisto->GetNbinsX(); ampbin++){
1016  if (PiAmpHisto->GetBinCenter(ampbin) < crossBandBegin ||
1017  PiAmpHisto->GetBinCenter(ampbin) > brossBandEnd){
1018  PiAmpGraphRight->SetPoint(PiAmpGraphRight->GetN(), PiAmpHisto->GetBinCenter(ampbin), PiAmpHisto->GetBinContent(ampbin));
1019  PiAmpGraphRight->SetPointError(PiAmpGraphRight->GetN()-1, PiAmpHisto->GetBinCenter(ampbin), PiAmpHisto->GetBinContent(ampbin));
1020  }
1021  }
1022 
1023 
1024  PiAmpGraphRight->Fit("PiFcnRight", "WR+");
1025 
1026  PiAmpHisto->GetListOfFunctions()->Add(PiAmpGraphRight->GetFunction("PiFcnRight"));
1027 
1028 
1029  if (PiFcnCenter) delete PiFcnCenter;
1030  if (PiFcnLeft) delete PiFcnLeft;
1031  if (PiFcnRight) delete PiFcnRight;
1032 
1033  //now fit E tail
1034  float EAmpFitBegin = EAmpHisto->GetBinCenter(EAmpHisto->GetMaximumBin())+0.03;//was +0.04
1035 
1036 
1037  if (nhitsPosition <= 15. && nhitsPosition >10) EAmpFitBegin +=0.05; //peak shift.
1038  if (nhitsPosition <= 10. && nhitsPosition >5) EAmpFitBegin +=0.09; //peak shift.
1039 
1040 
1041  TF1* EFcnLeft = new TF1("EFcnLeft","expo",EAmpFitBegin,EAmpFitBegin+0.22);
1042  EFcnLeft->SetLineColor(2);
1043  TF1* EFcnRight = new TF1("EFcnRight","expo",0.35,0.8);
1044  EFcnRight->SetLineColor(4);
1045 
1046  // EAmpHisto->Fit("EFcnLeft","RW+");
1047 
1048  crossBandBegin = 0.45;
1049  brossBandEnd = 0.6;
1050 
1051  TGraphErrors* EAmpGraphRight = new TGraphErrors();
1052 
1053  for (int ampbin = 0; ampbin<EAmpHisto->GetNbinsX(); ampbin++){
1054  if (EAmpHisto->GetBinCenter(ampbin) < crossBandBegin ||
1055  EAmpHisto->GetBinCenter(ampbin) > brossBandEnd){
1056  EAmpGraphRight->SetPoint(EAmpGraphRight->GetN(), EAmpHisto->GetBinCenter(ampbin), EAmpHisto->GetBinContent(ampbin));
1057  EAmpGraphRight->SetPointError(EAmpGraphRight->GetN()-1, EAmpHisto->GetBinCenter(ampbin), EAmpHisto->GetBinContent(ampbin));
1058  }
1059  }
1060 
1061 
1062  EAmpGraphRight->Fit("EFcnRight", "WR+");
1063 
1064  EAmpHisto->GetListOfFunctions()->Add(EAmpGraphRight->GetFunction("EFcnRight"));
1065 
1066 
1067 
1068 
1069 
1070  if (EFcnLeft) delete EFcnLeft;
1071  if (EFcnRight) delete EFcnRight;
1072 
1073 
1074  //now fit K
1075  low = 0.4;
1076  high = 0.9;
1077  ampPeakPos = 0.68;
1078 
1079  //find the peak position of K amp.
1080 
1081  /*
1082  for (int t=0; t<4; t++){
1083  TF1* gs = new TF1("gs","gaus",low,high);
1084  KAmpHisto->Fit("gs","NWR");
1085  low = gs->GetParameter(1)-0.1;
1086  high = gs->GetParameter(1)+0.1;
1087  ampPeakPos= gs->GetParameter(1);
1088  delete gs;
1089  }
1090  */
1091 
1092 
1093  TF1* KFcnCenter = new TF1("KFcnCenter","gaus",low,high);
1094 
1095  TF1* KFcnLeft = new TF1("KFcnLeft","pol3", (nhitsPosition <= 25.) ? 0.28 : 0.35, 0.65);
1096  KFcnLeft->SetLineColor(4);
1097 
1098  TF1* KFcnRight = new TF1("KFcnRight","expo",0.8,2.);
1099  KFcnRight->SetLineColor(3);
1100 
1101 
1102  crossBandBegin = 0.48;
1103  brossBandEnd = 0.58;
1104 
1105 
1106  TGraphErrors* KAmpGraphLeft = new TGraphErrors();
1107 
1108  for (int ampbin = 0; ampbin<KAmpHisto->GetNbinsX(); ampbin++){
1109  if (KAmpHisto->GetBinCenter(ampbin) < crossBandBegin ||
1110  KAmpHisto->GetBinCenter(ampbin) > brossBandEnd){
1111  KAmpGraphLeft->SetPoint(KAmpGraphLeft->GetN(), KAmpHisto->GetBinCenter(ampbin), KAmpHisto->GetBinContent(ampbin));
1112  KAmpGraphLeft->SetPointError(KAmpGraphLeft->GetN()-1, KAmpHisto->GetBinCenter(ampbin), KAmpHisto->GetBinContent(ampbin));
1113  }
1114  }
1115 
1116  KAmpGraphLeft->Fit("KFcnLeft", "WR+");
1117  KAmpHisto->GetListOfFunctions()->Add(KAmpGraphLeft->GetFunction("KFcnLeft"));
1118  crossBandBegin = 0.9;
1119  brossBandEnd = 1.5;
1120 
1121 
1122  TGraphErrors* KAmpGraphRight = new TGraphErrors();
1123 
1124  for (int ampbin = 0; ampbin<KAmpHisto->GetNbinsX(); ampbin++){
1125  if (KAmpHisto->GetBinCenter(ampbin) < crossBandBegin ||
1126  KAmpHisto->GetBinCenter(ampbin) > brossBandEnd){
1127  KAmpGraphRight->SetPoint(KAmpGraphRight->GetN(), KAmpHisto->GetBinCenter(ampbin), KAmpHisto->GetBinContent(ampbin));
1128  KAmpGraphRight->SetPointError(KAmpGraphRight->GetN()-1, KAmpHisto->GetBinCenter(ampbin), KAmpHisto->GetBinContent(ampbin));
1129  }
1130  }
1131 
1132 
1133  KAmpGraphRight->Fit("KFcnRight", "WR+");
1134 
1135  KAmpHisto->GetListOfFunctions()->Add(KAmpGraphRight->GetFunction("KFcnRight"));
1136 
1137  KAmpHisto->Fit("KFcnCenter","WR+");
1138 
1139  //for the part overlay with FcnLeft, use FcnLeft
1140  if ( KFcnCenter->GetXmin()<KFcnLeft->GetXmax() && KFcnLeft->GetXmax() < KFcnCenter->GetXmax() )
1141  KAmpHisto->GetFunction("KFcnCenter")->SetRange(KFcnLeft->GetXmax()-0.05, KFcnCenter->GetXmax());
1142 
1143  // KAmpHisto->GetListOfFunctions()->Add(KAmpHisto->GetFunction("KFcnCenter"));
1144 
1145  if (KFcnCenter) delete KFcnCenter;
1146  if (KFcnLeft) delete KFcnLeft;
1147  if (KFcnRight) delete KFcnRight;
1148 
1149 
1150  //now fit P
1151 
1152  TF1* PFcnLeft = new TF1("PFcnLeft","pol4",0.6,1.4);
1153  TF1* PFcnRight = new TF1("PFcnRight","expo",1.4,2.);
1154 
1155  crossBandBegin = 1.5;
1156  brossBandEnd = 1.8;
1157 
1158 
1159  TGraphErrors* PAmpGraphRight = new TGraphErrors();
1160 
1161  for (int ampbin = 0; ampbin<PAmpHisto->GetNbinsX(); ampbin++){
1162  if (PAmpHisto->GetBinCenter(ampbin) < crossBandBegin ||
1163  PAmpHisto->GetBinCenter(ampbin) > brossBandEnd){
1164  PAmpGraphRight->SetPoint(PAmpGraphRight->GetN(), PAmpHisto->GetBinCenter(ampbin), PAmpHisto->GetBinContent(ampbin));
1165  PAmpGraphRight->SetPointError(PAmpGraphRight->GetN()-1, PAmpHisto->GetBinCenter(ampbin), PAmpHisto->GetBinContent(ampbin));
1166  }
1167  }
1168 
1169 
1170  PAmpGraphRight->Fit("PFcnRight", "WR+");
1171 
1172  PAmpHisto->GetListOfFunctions()->Add(PAmpGraphRight->GetFunction("PFcnRight"));
1173 
1174 
1175  PAmpHisto->Fit("PFcnLeft","RW+");
1176 
1177  if (PFcnRight) delete PFcnRight;
1178  if (PFcnLeft) delete PFcnLeft;
1179 
1180  //write hito.
1181 
1182  ampHistFile->cd();
1183  PiAmpHisto->Write(PiAmpHisto->GetName(),TObject::kOverwrite | TObject::kSingleKey);
1184  EAmpHisto->Write(EAmpHisto->GetName(),TObject::kOverwrite | TObject::kSingleKey);
1185  KAmpHisto->Write(KAmpHisto->GetName(),TObject::kOverwrite | TObject::kSingleKey);
1186  PAmpHisto->Write(PAmpHisto->GetName(),TObject::kOverwrite | TObject::kSingleKey);
1187 
1188 
1189  }
1190 
1191 ampHistFile->Write();
1192 ampHistFile->Close();
1193 gausHistFile->Close();
1194  cout<<" **********extraping pion amp. end********* "<<endl;
1195 }
1196 
1197 
1198 Double_t sigmaNSampleFitFcn(Double_t* x, Double_t *par){
1199 
1200  if (x[0]) return par[0]/::sqrt(x[0]);
1201  else return 0.;
1202 }
1203 
1205 //
1206 // $Log: PIDFitter.cxx,v $
1207 // Revision 1.11 2007/07/12 19:47:31 fisyak
1208 // Add includes for ROOT 5.16
1209 //
1210 // Revision 1.10 2003/09/07 03:49:04 perev
1211 // gcc 3.2 + WarnOff
1212 //
1213 // Revision 1.9 2003/09/02 17:58:46 perev
1214 // gcc 3.2 updates + WarnOff
1215 //
1216 // Revision 1.8 2003/08/04 18:46:26 perev
1217 // warnOff
1218 //
1219 // Revision 1.7 2003/07/25 15:04:17 aihong
1220 // change fit range etc.
1221 //
1222 // Revision 1.4 2002/12/17 17:05:48 aihong
1223 // write out sigmaNsample graph
1224 //
1225 // Revision 1.3 2002/03/27 14:19:21 aihong
1226 // use two functions to fit e+/- amp. instead of one
1227 //
1228 // Revision 1.2 2002/02/25 18:31:24 jeromel
1229 // More SetFormat() stripped.
1230 //
1231 // Revision 1.1 2002/02/14 21:25:55 aihong
1232 // re-install the new version
1233 //
1234 
1235 //
PIDFitter::kaonAntiprotonBandCenter
TF1 * kaonAntiprotonBandCenter
//for drawing line between pion and kaon bands
Definition: PIDFitter.h:58
PIDFitter::BBOffSetPar
double ** BBOffSetPar
//for drawing line between kaon and antiproton bands
Definition: PIDFitter.h:60
Generated by   doxygen 1.8.5