PWGLF/STRANGENESS/LambdaK0PbPb/FitLambdaSpectrum.C

   1 static  int      myAliceBlue    = TColor::GetColor(0,0,192);
   2 static  int      myAliceRed    = TColor::GetColor(192,0,0);
   3 void FitLambdaSpectrum(const char* name, const char * listName = "clambdak0Histo_00", const char* HistName = "h2PtVsMassLambda", const char* SaveName = "SPECTRA_MK_Lambda_7TeV.root")
   4
   5 ////////////////////////zaciatok FitLambdaSpectrum//////////////////////////////////////////////////////////
   6 {
   7
   8   gROOT->SetStyle("Plain");
   9   gROOT->LoadMacro("run.C");
  10   InitAndLoadLibs();
  11
  12   TFile *myFile = new TFile(name);
  13   myFile->ls();
  14   TList *myList = myFile->Get(listName); // FIXME: find suffix?
  15
  16   if(!myList) {
  17     cout << "Cannot get list " << listName << " from file " << name << endl;
  18     return;
  19   }
  20
  21   TCanvas *myCan = new TCanvas("spectra","lambda");
  22   myCan->Draw();
  23   myCan->cd();
  24
  25   TLatex *myKind = new TLatex(0.15,0.92,"ALICE: LHC10b Pass1 7 TeV"); // FIXME
  26   myKind->SetNDC();
  27   myKind->SetTextSize(0.03);
  28   myKind->SetTextColor(myAliceBlue);
  29   myKind->Draw();
  30
  31
  32   TText *textTitle = new TText(0.6,0.86,"");
  33   textTitle->SetNDC();
  34   textTitle->SetTextSize(0.04);
  35   textTitle->Draw();
  36
  37   TPad *myPad1 = new TPad("myPad1","myPad1",0,0,1,0.85);
  38   myPadSetUp(myPad1,0.12,0.02,0.02,0.15);
  39   myPad1->Draw();
  40   myPad1->cd();
  41
  42   //    gStyle->SetOptTitle(0);
  43   //    gStyle->SetOptStat(0);
  44
  45   TH1F * SaveParameters= new TH1F("SaveParameters","parameters",6,0,6);
  46
  47   ///zaciatok binovania/////////////
  48   const Int_t range=33;
  49
  50   Double_t b[range];
  51
  52   b[0]=0;
  53   b[1]=0.5;
  54   for (Int_t i = 2;i<17;i++)
  55     {
  56       b[i]=b[i-1]+0.1;
  57       cout<<b[i]<<endl;
  58     }
  59
  60   for (Int_t i = 17;i<27;i++)
  61     {
  62       b[i]=b[i-1]+0.2;
  63       cout<<b[i]<<endl;
  64     }
  65   b[27]=4.5;
  66   b[28]=5.0;
  67   b[29]=5.5;
  68   b[30]=6.5;
  69   b[31]=8.0;
  70   b[32]=12.0;
  71
  72
  73   //koniec binovania
  74   TH1F *DrawSpectrumLambda = new TH1F("DrawSpectrumLambda","#Lambda Spectrum;p_{t} [GeV/c]; N",32,b);
  75
  76   ///////////////////////////////////////////////run over all bins//////////////////////////////////////////////////////////////
  77   int iLoBin = 15, iHiBin = 16;
  78   //    int iLoBin = 71, iHiBin = 72, hMax = 250;
  79
  80   for (Int_t rBin = 4;rBin<33;rBin++){
  81     //for (Int_t rBin = 36;rBin<51;rBin++){
  82
  83
  84     myBinCounting(myList, DrawSpectrumLambda, rBin,iLoBin ,iHiBin,HistName,SaveParameters);
  85
  86     if ((rBin>=1)&&(rBin<16)) {iLoBin = iLoBin+2; iHiBin = iHiBin+2;}
  87
  88     if  (rBin==16) {iLoBin = iLoBin+2; iHiBin = iHiBin+4;}
  89
  90     if ((rBin>=17)&&(rBin<26)){iLoBin = iLoBin+4; iHiBin = iHiBin+4;}
  91
  92     if  (rBin==26) {iLoBin = iLoBin+4; iHiBin = iHiBin+10;}
  93
  94     if  (rBin==27) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;}
  95
  96     if  (rBin==28) {iLoBin = iLoBin+10; iHiBin = iHiBin+10;}
  97
  98     if  (rBin==29) {iLoBin = iLoBin+10; iHiBin = iHiBin+20;}
  99
 100     if  (rBin==30) {iLoBin = iLoBin+20; iHiBin = iHiBin+30;}
 101
 102     if  (rBin==31) {iLoBin = iLoBin+30; iHiBin = iHiBin+80;}
 103
 104
 105
 106     char saveFileName[60];
 107     char saveFileName1[60];
 108     sprintf(saveFileName,"hm_InvariantMassLambdafit_bin%d.gif",rBin);
 109     sprintf(saveFileName1,"hm_InvariantMassLambdafit_bin%d.root",rBin);
 110     myCan->SaveAs(saveFileName);
 111     myCan->SaveAs(saveFileName1);
 112     myCan->Clear();
 113
 114   }
 115   //////////////////////////koniec cyklu cez jednotlive biny/////////////////////////////
 116
 117
 118
 119   //////////////////////////vykreslenie a ulozenie histogramu///////////////////////////
 120   DrawSpectrumLambda->Draw();
 121
 122   TFile *SPECTRA_MK_Lambda = new TFile(SaveName,"RECREATE");
 123   SPECTRA_MK_Lambda->WriteObject(DrawSpectrumLambda,"DrawSpectrumLambda");
 124
 125   SPECTRA_MK_Lambda->Close();
 126   //////////////////////////koniec ulozenia a vykreslenia/////////////////////////////////
 127
 128
 129 }
 130 /////////end FitLambdaSpectrum///////////////////////////////////////////////////////////////////////////////
 131
 132
 133
 134
 135
 136
 137 //////////////////////funkcia pre fitovanie////////////////////////////////////////////////////////
 138
 139
 140 void myBinCounting(TList *myList = 0, TH1F *DrawSpectrumLambda = 0, Int_t rBin = 0, Int_t iLoBin =0, Int_t iHiBin =0, const char* HistName =0,TH1F *SaveParameters = 0){
 141
 142
 143   TH2F *h2PtVsMassLambda = (TH2F*)myList->FindObject(HistName);
 144   TH1D *h1MassLambdaPtBin = h2PtVsMassLambda->ProjectionX("h1MassLambdaPtBin",iLoBin,iHiBin);
 145   h1MassLambdaPtBin->SetFillColor(0);
 146   //    h1MassLambdaPtBin->SetTitle("h1MassLambdaPtBin;invariant mass (GeV/c^{2});Counts");
 147   //h1MassLambdaPtBin->Draw("SAME");
 148
 149
 150   ///////////////////////////////////////Gauss+2nd pol fit////////////////////////////////////////////
 151
 152   TH1D *myHistGauss = h1MassLambdaPtBin;
 153   //Double_t hMax=myHistGauss->GetMaximum();
 154   //myHistGauss->SetMaximum(hMax);
 155   myHistGauss->Draw();
 156   //   double myLevel = 500, mySlope = 0, mySlope2 = 0 ;
 157   //   double myNorm  = 100, myMean = 0.497, myWidth = 0.003;
 158   double myRange[2]={0,0};
 159
 160   // Gaussain Fit
 161
 162   double    myLevel = 0, mySlope = 0, mySlope2 = 0;
 163   double     myNorm  = 0, myMean  = 1.116, myWidth = 0.003;
 164   //    myNorm=(myHistGauss->GetMaximum()*TMath::Sqrt(2*TMath::Pi()));
 165   myNorm=myHistGauss->GetMaximum();
 166   //     myRange[0]=1.09; myRange[1]=1.14;
 167   myRange[0]=1.095; myRange[1]=1.135;
 168
 169   if(rBin==4)
 170     {
 171       //myLevel =-2.46E5; mySlope = 3.39E5; mySlope2 = -1.03E5; myNorm  = 5.12E3;
 172       myRange[0]=1.10; myRange[1]=1.14;
 173     }
 174
 175   if(rBin==3)
 176     {
 177       myLevel =-2.91E6; mySlope = 5.21E6; mySlope2 = -2.33E6; myNorm  = 9.32E3;
 178       myRange[0]=1.104; myRange[1]=1.14;
 179     }
 180
 181   char cRange[30],cInterval[30],cIntegration[30],cFitInfo[100], cFitInfo2[100], cFitInfo3[100];
 182
 183   TF1 *fPolGauss = new TF1("fpolgauss",myPolGauss,myRange[0],myRange[1],6);
 184   fPolGauss->SetTitle("The gaussian function");
 185   fPolGauss->SetParNames("level","slope","slope2","norm","mean","width");
 186   fPolGauss->SetParameters(myLevel,mySlope,mySlope2,myNorm,myMean,myWidth);
 187   fPolGauss->SetLineStyle(2);
 188   fPolGauss->SetLineWidth(2);
 189   fPolGauss->SetLineColor(602);
 190   ///    fPolGauss->Draw("SAME");
 191   myHistGauss->Fit(fPolGauss,"IREM");
 192   fPolGauss->Draw("SAME");
 193
 194   //myHistGauss->Fit(fPolGauss,"LIREM");
 195
 196   myNorm  = fPolGauss->GetParameter(3);
 197   myMean  = fPolGauss->GetParameter(4);
 198   myWidth = fPolGauss->GetParameter(5);
 199
 200
 201
 202   /////store parameters for next fit if own will fall////
 203   //if ((myMean-(4*myWidth))<(myMean+(4*myWidth)))
 204   //if (myWidth>0.0005)
 205   if ((myWidth>0.0005)&&(myWidth<0.004))
 206     {
 207       for (Int_t i=1;i<7;i++) SaveParameters->SetBinContent(i,0);
 208       SaveParameters->SetBinContent(1,fPolGauss->GetParameter(0));
 209       SaveParameters->SetBinContent(2,fPolGauss->GetParameter(1));
 210       SaveParameters->SetBinContent(3,fPolGauss->GetParameter(2));
 211       SaveParameters->SetBinContent(4,fPolGauss->GetParameter(3));
 212       SaveParameters->SetBinContent(5,fPolGauss->GetParameter(4));
 213       SaveParameters->SetBinContent(6,fPolGauss->GetParameter(5));
 214
 215
 216     }
 217   ///////////////////////////////////////////////////////
 218
 219   ///////////if parameters are wrong use previous/////
 220   //if ((myMean-(4*myWidth))>=(myMean+(4*myWidth)))
 221   //if (myWidth<=0.0005)
 222   if ((myWidth<=0.0005)||(myWidth>=0.004))
 223     {
 224       cout<<"Previous parameters"<<endl;
 225       for (Int_t i=1;i<7;i++) cout<<SaveParameters->GetBinContent(i)<<endl;
 226       myLevel=SaveParameters->GetBinContent(1);
 227       mySlope=SaveParameters->GetBinContent(2);
 228       mySlope2=SaveParameters->GetBinContent(3);
 229       myNorm=SaveParameters->GetBinContent(4);
 230       myMean=SaveParameters->GetBinContent(5);
 231       myWidth=SaveParameters->GetBinContent(6);
 232
 233       cout<<"nova sirka"<<myWidth<<endl;
 234     }
 235   ////////////////////////////////////////////////////
 236
 237
 238   TF1 *fGauss = new TF1("fgauss",myGauss,myRange[0],myRange[1],3);
 239   fGauss->SetTitle("The gaussian function");
 240   fGauss->SetParNames("norm","mean","width");
 241   fGauss->SetParameters(myNorm,myMean,myWidth);
 242   fGauss->SetLineColor(myAliceRed);
 243   fGauss->SetFillColor(myAliceRed);
 244   fGauss->SetFillStyle(3005);
 245   fGauss->SetLineWidth(1);
 246   fGauss->Draw("SAME");
 247
 248   int    NumSigmaSignal      = 4;
 249   int    NumSigmaBackground   =6;
 250   double integralWidth = NumSigmaSignal*myWidth;
 251
 252   double binWidth      = myHistGauss->GetBinWidth(1);
 253   double integralGauss = fGauss->Integral(myMean-integralWidth,myMean+integralWidth);
 254   sprintf(cRange,"%d#sigma range (GeV/c^{2})",NumSigmaSignal);
 255   sprintf(cInterval,"[%.3f;%.3f]",myMean-integralWidth,myMean+integralWidth);
 256   sprintf(cIntegration,"integral: %.0f",integralGauss/binWidth);
 257   sprintf(cFitInfo,"#Chi^{2}/ndf = %.1f/%d",fPolGauss->GetChisquare(),fPolGauss->GetNDF());
 258   //     sprintf(cFitInfo2,"mean = %.4f #pm %.4f",fPolGauss->GetParameter(4),fPolGauss->GetParError(4));
 259   //     sprintf(cFitInfo3,"width = %.5f #pm %.5f",fPolGauss->GetParameter(5),fPolGauss->GetParError(5));
 260
 261   cout<<cRange<<endl;
 262   cout<<cInterval<<endl;
 263   cout<<cIntegration<<endl;
 264   cout<<cFitInfo<<endl;
 265   //cout<<cFitInfo2<<endl;
 266   //cout<<cFitInfo3<<endl;
 267
 268
 269   //break;
 270
 271   ////////////////////////////////////////koniec Gauss+2nd pol fit////////////////////////////////////
 272
 273
 274   Double_t GaussMean =myMean; //fPolGauss->GetParameter(4); //stredna hodnota gausovskeho fitu
 275   Double_t GaussSigma =myWidth; //fPolGauss->GetParameter(5); //sigma gausovskeho fitu
 276
 277   Double_t LeftSide = myHistGauss->FindBin(GaussMean - NumSigmaSignal*GaussSigma);  // hranice signalu
 278   Double_t RightSide = myHistGauss->FindBin(GaussMean + NumSigmaSignal*GaussSigma);
 279
 280   Double_t NumberBinsSignal = (RightSide - LeftSide) + 1; //pocet binov ktore obsahuju signal
 281
 282   Double_t LeftB_Side_Left = myHistGauss->FindBin(myRange[0]); //bacground on left side of peak
 283   Double_t LeftB_Side_Right = myHistGauss->FindBin(GaussMean - NumSigmaBackground*GaussSigma);
 284
 285   Double_t RightB_Side_Left = myHistGauss->FindBin(GaussMean + NumSigmaBackground*GaussSigma); //bacground on Right side of peak
 286   Double_t RightB_Side_Right = myHistGauss->FindBin(myRange[1]);
 287
 288
 289
 290   Double_t NumberBinsBackground = (LeftB_Side_Right - LeftB_Side_Left) + 1 + (RightB_Side_Right - RightB_Side_Left) + 1; // pocet binov co
 291
 292
 293   //
 294   // Bin Counting Method
 295   //
 296
 297   // Compute the whole content for Signal in 4 sigma
 298   Float_t SignalContentCount = h1MassLambdaPtBin->Integral(LeftSide,RightSide );
 299
 300
 301   // Create a histo with the Noise only.
 302
 303   Double_t  nCountNoise   = 0;
 304   TH1D *myHistNoise = new TH1D(*h1MassLambdaPtBin);
 305   myHistNoise->SetName("myHistNoise");
 306   myHistNoise->SetTitle("myHistNoise");
 307   myHistNoise->Scale(0);
 308   myHistNoise->SetFillColor(myAliceBlue);
 309
 310
 311
 312   for (Int_t iBin = LeftB_Side_Left; iBin <= LeftB_Side_Right; iBin++) {
 313
 314     nCountNoise      +=h1MassLambdaPtBin ->GetBinContent(iBin);
 315     myHistNoise->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin));
 316     myHistNoise->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin));
 317   }
 318
 319   for (Int_t iBin = RightB_Side_Left; iBin <= RightB_Side_Right; iBin++) {
 320
 321     nCountNoise      +=h1MassLambdaPtBin ->GetBinContent(iBin);
 322     myHistNoise->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin));
 323     myHistNoise->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin));
 324   }
 325
 326   // Create a histo with the Signal in 4 sigma.
 327   int totInHistSignal = 0;
 328   TH1D *myHistSignal = new TH1D(*h1MassLambdaPtBin);
 329   myHistSignal->SetName("myHistSignal");
 330   myHistSignal->SetTitle("myHistSignal");
 331   myHistSignal->Scale(0);
 332   myHistSignal->SetFillStyle(3004);
 333
 334
 335   for (Int_t iBin = LeftSide; iBin <= RightSide; iBin++){
 336     myHistSignal->SetBinContent(iBin,h1MassLambdaPtBin->GetBinContent(iBin));
 337     myHistSignal->SetBinError(iBin,h1MassLambdaPtBin->GetBinError(iBin));
 338     totInHistSignal += h1MassLambdaPtBin->GetBinContent(iBin);
 339
 340   }
 341
 342
 343   //
 344   // Background Fit
 345   //
 346
 347
 348   cout<<"zaciatok background fit casti "<<endl;
 349   cout<<endl;
 350   cout<<endl;
 351
 352   //   if  (rOrderPoly == 1 ) {
 353   //     TF1*fitNoise = new TF1("fitNoise",myPol1,LeftB_Side_Left,RightB_Side_Right,4);
 354   //
 355   //     fitNoise->FixParameter(2,(GaussMean - 6*GaussSigma));
 356   //     fitNoise->FixParameter(3,(GaussMean + 6*GaussSigma));
 357   //   }
 358
 359
 360   TF1*fitNoise = new TF1("fitNoise",myPol2,myRange[0],myRange[1],5);
 361   fitNoise->SetParLimits(2,-1E8,0);
 362   fitNoise->SetParameter(0,fPolGauss->GetParameter(0));
 363   fitNoise->SetParameter(1,fPolGauss->GetParameter(1));
 364   fitNoise->SetParameter(2,fPolGauss->GetParameter(2));
 365
 366   fitNoise->FixParameter(3,(GaussMean - NumSigmaBackground*GaussSigma));
 367   fitNoise->FixParameter(4,(GaussMean + NumSigmaBackground*GaussSigma));
 368
 369
 370
 371   Double_t inteGrapmyHistNoise = 0;
 372   Double_t inteGrapmyHistSignal = 0;
 373
 374
 375   myHistNoise->Fit("fitNoise","irem+");
 376   //myHistNoise->Fit("fitNoise","Lirem+");
 377
 378   Double_t hMax=myHistGauss->GetMaximum();
 379   myHistNoise->SetMaximum(hMax);
 380
 381   myHistNoise->GetYaxis()->SetTitle("counts");
 382   myHistNoise->GetYaxis()->SetTitleOffset(1.2);
 383   myHistNoise->SetMaximum(hMax);
 384   myHistNoise->Draw("HIST");
 385   myHistSignal->SetFillColor(myAliceRed);
 386   myHistSignal->Draw("A BAR E SAME");
 387   myHistGauss->SetLineWidth(1);
 388   myHistGauss->Draw("A H E SAME");
 389
 390   //fitNoise->Draw("same");
 391
 392   //****************************************
 393   // Bin Counting
 394   //****************************************
 395
 396   // Noise Under Peak
 397
 398   //     TF1*fNoise = new TF1("fNoise","[0]+x*[1]",LeftB_Side_Left,RightB_Side_Right);
 399   //    fNoise->FixParameter(0,fitNoise->GetParameter(0));
 400   //     fNoise->FixParameter(1,fitNoise->GetParameter(1));
 401
 402
 403
 404   TF1*fNoise = new TF1("fNoise","[0]+x*[1]+x*x*[2]",myRange[0],myRange[1]);
 405   fNoise->FixParameter(0,fitNoise->GetParameter(0));
 406   fNoise->FixParameter(1,fitNoise->GetParameter(1));
 407   fNoise->FixParameter(2,fitNoise->GetParameter(2));
 408
 409
 410   fNoise->Draw("same");
 411
 412   //   Double_t nNoiseUnderPeak    = fNoise->Integral(GaussMean - NumSigmaSignal*GaussSigma,GaussMean + NumSigmaSignal*GaussSigma)/(myHistNoise->GetBinWidth(1));
 413   Double_t lBinWidth = myHistNoise->GetBinWidth(1);
 414   cout<<"sirka binu"<< lBinWidth<<endl;
 415   Double_t nNoiseUnderPeak    = fNoise->Integral(h1MassLambdaPtBin->GetBinCenter(LeftSide)-0.5*lBinWidth,h1MassLambdaPtBin->GetBinCenter(RightSide)+0.5*lBinWidth)/lBinWidth;
 416
 417
 418   Double_t nNoiseUnderPeakErr = TMath::Sqrt(nNoiseUnderPeak);
 419
 420
 421   // Signal
 422   Double_t   Signal[2]; // Signal[0]=signal and Signal[1]=error
 423   Signal[0] = Signal[1]=0.0;
 424   Signal[0] = totInHistSignal-nNoiseUnderPeak;
 425   Signal[1] = TMath::Sqrt(Signal[0]+nNoiseUnderPeak);
 426
 427   //    DrawSpectrumLambda->SetBinContent(rBin,Signal[0]/DrawSpectrumLambda->GetBinWidth(rBin);
 428   //    DrawSpectrumLambda->SetBinError(rBin,Signal[1]/DrawSpectrumLambda->GetBinWidth(rBin));
 429   DrawSpectrumLambda->SetBinContent(rBin,Signal[0]);
 430   DrawSpectrumLambda->SetBinError(rBin,Signal[1]);
 431
 432   printf("BoInfo: Noise intervals [%.3f;%.3f],[%.3f;%.3f]\n",myRange[0],(GaussMean - NumSigmaBackground*GaussSigma),(GaussMean +NumSigmaBackground *GaussSigma),myRange[1]);
 433   printf("BoInfo: nBinNoise(left+right)=%d  nCountNoise(left+right)=%f \n",NumberBinsBackground,nCountNoise);
 434   printf("BoInfo: nBinCentral=%d  totInHistSignal=%f \n",NumberBinsSignal,totInHistSignal);
 435   printf("BoInfo: Noise under peak %.2f +-%.2f \n",nNoiseUnderPeak, nNoiseUnderPeakErr);
 436   printf("BoInfo: Signal %f \n",Signal[0]);
 437
 438
 439   // printf only and keep only on plot the Signal, noise and S/N.
 440   Char_t     cIntervalSignal[50], cSignalOverNoise[50], cNoiseUnderPeak[50], cFitNoiseInfo[50];
 441   Double_t   SignalOverNoise    = 0.0;
 442   Double_t   SignalOverNoiseErr = 0.0;
 443
 444   sprintf(cIntervalSignal,"Signal [%.3f;%.3f] = %.0f #pm %.0f",LeftSide,RightSide,Signal[0], Signal[1]);
 445   sprintf(cNoiseUnderPeak,"Noise under peak = %.0f #pm %.0f",nNoiseUnderPeak, nNoiseUnderPeakErr);
 446   SignalOverNoise    = Signal[0]/nNoiseUnderPeak;
 447   SignalOverNoiseErr = TMath::Sqrt(Signal[0]+2*nNoiseUnderPeak);
 448   sprintf(cSignalOverNoise,"S/N= %.2f #pm %.2f",SignalOverNoise,SignalOverNoiseErr);
 449   cout<<"cSignalOverNoise"<<cSignalOverNoise<<endl;
 450   sprintf(cFitNoiseInfo,"#Chi^{2}/ndf = %.1f/%d",fitNoise->GetChisquare(),fitNoise->GetNDF());
 451
 452   cout<<cIntervalSignal<<endl;
 453   cout<<cNoiseUnderPeak<<endl;
 454   cout<<cFitNoiseInfo<<endl;
 455
 456 }
 457
 458 /////////////////////////////////////koniec funkcie////////////////////////////////////////
 459
 460
 461
 462 /////////////////////////////////////////funkcie///////////////////////////////////////////////////////
 463
 464
 465 void myPadSetUp(TPad *currentPad=0, float rLeft=0, float rTop=0, float rRight=0, float rBottom=0){
 466   currentPad->SetLeftMargin(rLeft);
 467   currentPad->SetTopMargin(rTop);
 468   currentPad->SetRightMargin(rRight);
 469   currentPad->SetBottomMargin(rBottom);
 470   return;
 471 }
 472
 473
 474 double myGauss(double *x, double *par){
 475   double absc  = x[0];
 476   double norm  = par[0];
 477   double mean  = par[1];
 478   double width = par[2];
 479   //  double mygauss = (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
 480   double mygauss = (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
 481   return mygauss;
 482 }
 483
 484 double myPolGauss(double *x, double *pars){
 485   double absc   = x[0];
 486   double level  = pars[0];
 487   double slope  = pars[1];
 488   double slope2 = pars[2];
 489   double norm   = pars[3];
 490   double mean   = pars[4];
 491   double width  = pars[5];
 492   //  double ordo = level + slope * absc +slope2 * absc*absc + (norm/TMath::Sqrt(2*TMath::Pi()))*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
 493   double ordo = level + slope * absc +slope2 * absc*absc + (norm)*TMath::Exp(-0.5*(absc-mean)*(absc-mean)/(width*width));
 494   return ordo;
 495 }
 496
 497 double myPol1(double *x, double *pars) {
 498
 499   double absc    = x[0];
 500   double level   = pars[0];
 501   double slope   = pars[1];
 502   double abscMin = pars[2];
 503   double abscMax = pars[3];
 504
 505   if (absc >= abscMin && absc < abscMax) {
 506     TF1::RejectPoint();
 507     return 0;
 508   }
 509
 510   double ordo = level + slope * absc;
 511   return ordo;
 512
 513 }
 514
 515 double myPol2(double *x, double *pars) {
 516
 517   double absc    = x[0];
 518   double level   = pars[0];
 519   double slope   = pars[1];
 520   double slope2  = pars[2];
 521   double abscMin = pars[3];
 522   double abscMax = pars[4];
 523
 524   if (absc >= abscMin && absc < abscMax) {
 525     TF1::RejectPoint();
 526     return 0;
 527   }
 528
 529   double ordo = level + slope * absc +slope2 * absc*absc;
 530   return ordo;
 531
 532 }
 533
 534
 535
 536 void LoadLibs(){
 537
 538   gSystem->Load("libCore");
 539   gSystem->Load("libTree");
 540   gSystem->Load("libGeom");
 541   gSystem->Load("libVMC");
 542   gSystem->Load("libPhysics");
 543   gSystem->Load("libSTEERBase");
 544   gSystem->Load("libESD");
 545   gSystem->Load("libAOD");
 546   gSystem->Load("libANALYSIS");
 547   gSystem->Load("libANALYSISalice");
 548 }