PWGLF/SPECTRA/Nuclei/B2/AliLnPt.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 // reconstructed pt
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #include <Riostream.h>
  20 #include <TFile.h>
  21 #include <TSystem.h>
  22 #include <TH1D.h>
  23 #include <TH2D.h>
  24 #include <TF1.h>
  25 #include <TString.h>
  26 #include <TROOT.h>
  27 #include <TString.h>
  28 #include <RooWorkspace.h>
  29 #include <RooRealVar.h>
  30 #include <RooDataSet.h>
  31 #include <RooDataHist.h>
  32 #include <RooHistPdf.h>
  33 #include <RooMsgService.h>
  34 #include <TMath.h>
  35
  36 #include "AliLnPt.h"
  37 #include "B2.h"
  38
  39 ClassImp(AliLnPt)
  40
  41 AliLnPt::AliLnPt(const TString& particle, Double_t trigEff, const TString& inputFilename, const TString& outputFilename, const TString& otag, const TString& corrFilename, const TString& corrtag)
  42 : TObject()
  43 , fParticle(particle)
  44 , fTrigEff(trigEff)
  45 , fInputFilename(inputFilename)
  46 , fOutputFilename(outputFilename)
  47 , fOutputTag(otag)
  48 , fCorrFilename(corrFilename)
  49 , fCorrTag(corrtag)
  50 , fPtMin(0.5)
  51 , fPtMax(3.0)
  52 , fPid(0)
  53 , fSecondaries(1)
  54 , fEfficiency(1)
  55 , fIsOnlyGen(0)
  56 , fPtPid(1.)
  57 , fBkgMin(2.2)
  58 , fBkgMax(5.)
  59 , fIntMin(2.)
  60 , fIntMax(6.)
  61 , fMakeStats(1)
  62 , fMCtoINEL(0)
  63 , fVtxFactor(1)
  64 , fFitFrac(1)
  65 , fFdwnCorr(1)
  66 , fSameFdwn(0)
  67 , fPidProc(kMassSquared)
  68 , fPidEff(1)
  69 , fDebugLevel(0)
  70 {
  71 //
  72 // constructor
  73 //
  74         TH1::SetDefaultSumw2(); // switch on histogram errors
  75
  76         if(fDebugLevel < 3)
  77         {
  78                 // disable verbose in RooFit
  79                 RooMsgService::instance().setSilentMode(1);
  80                 RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);
  81         }
  82 }
  83
  84 AliLnPt::~AliLnPt()
  85 {
  86 //
  87 // destructor
  88 //
  89 }
  90
  91 Int_t AliLnPt::Exec()
  92 {
  93 //
  94 // Get the true pt distribution
  95 //
  96         TFile* finput = new TFile(fInputFilename.Data(), "read");
  97         if (finput->IsZombie()) exit(1);
  98
  99         TFile* foutput = new TFile(fOutputFilename.Data(), "recreate");
 100         if (foutput->IsZombie()) exit(1);
 101
 102         TFile* fcorr = 0;
 103         TFile* fdebug = 0;
 104
 105         TString species = fParticle;
 106         species.ReplaceAll("Anti","");
 107
 108         TH1D* hStats = FindObj<TH1D>(finput, species + "_Stats");
 109
 110         TH1D* hCorrPt = 0;
 111
 112         if(fIsOnlyGen)
 113         {
 114                 hCorrPt = dynamic_cast<TH1D*>( (FindObj<TH1D>(finput, fParticle + "_Gen_PhS_Prim_Pt"))->Clone( Form("%s_Cor_Pt", fParticle.Data())) );
 115         }
 116         else
 117         {
 118                 fcorr = new TFile(fCorrFilename.Data(),"read");
 119                 if (fcorr->IsZombie()) exit(1);
 120
 121                 fdebug =  new TFile(Form("debug-%s",fOutputFilename.Data()),"recreate");
 122                 if (fdebug->IsZombie()) exit(1);
 123
 124                 fdebug->mkdir(fOutputTag.Data());
 125
 126                 // pointers to the required histograms
 127
 128                 TH1D* hOrigPt  = FindObj<TH1D>(finput, fParticle + "_PID_Pt");
 129
 130                 TH2D* hPidPt = 0;
 131                 if(fPidProc == kTimeOfFlight)        hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DTime_Pt");
 132                 else if(fPidProc == kMassSquared)    hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_M2_Pt");
 133                 else if(fPidProc == kMassDifference) hPidPt = FindObj<TH2D>(finput, fParticle + "_PID_DM2_Pt");
 134                 else fPid = 0;
 135
 136                 // corrections
 137                 TH1D* hFracMatPt   = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Pt");
 138                 TH1D* hFracFdwnPt  = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Pt");
 139                 TF1* fncFracMatPt  = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Mat_Fit_Pt");
 140                 TF1* fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, fParticle + "_Frac_Fdwn_Fit_Pt");
 141
 142                 if(fSameFdwn)
 143                 {
 144                         hFracFdwnPt = FindObj<TH1D>(fcorr, fCorrTag, species + "_Frac_Fdwn_Pt");
 145                         fncFracFdwnPt = FindObj<TF1>(fcorr, fCorrTag, species + "_Frac_Fdwn_Fit_Pt");
 146                 }
 147
 148                 TH1D* hEffVtxPt    = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_Vtx_Pt");
 149                 TH1D* hEffAccTrkPt = FindObj<TH1D>(fcorr, fCorrTag, fParticle + "_Eff_AccTrk_Pt");
 150
 151                 // apply corrections
 152
 153                 fdebug->cd(fOutputTag.Data());
 154
 155                 TH1D* hPidCorrPt = 0;
 156                 if(fPid)
 157                 {
 158                         if(fDebugLevel>1) std::cout << "\t\tpid correction" << std::endl;
 159
 160                         hPidCorrPt = this->PID(hOrigPt, hPidPt, fPtPid, fPtMax, fParticle + "_PidCorr_Pt");
 161                         hPidCorrPt->Scale(1./fPidEff);
 162                 }
 163                 else
 164                 {
 165                         hPidCorrPt = dynamic_cast<TH1D*>(hOrigPt->Clone(Form("%s_PidCorr_Pt", fParticle.Data())));
 166                 }
 167
 168                 TH1D* hSecCorrPt = 0;
 169                 if(fSecondaries)
 170                 {
 171                         if(fDebugLevel>1) std::cout << "\t\tremoval of secondaries" << std::endl;
 172
 173                         if(fFitFrac)
 174                         {
 175                                 hSecCorrPt = this->Secondaries(hPidCorrPt, fncFracMatPt, fncFracFdwnPt, fParticle + "_SecCorr_Pt");
 176                         }
 177                         else
 178                         {
 179                                 hSecCorrPt = this->Secondaries(hPidCorrPt, hFracMatPt, hFracFdwnPt, fParticle + "_SecCorr_Pt");
 180                         }
 181                 }
 182                 else
 183                 {
 184                         hSecCorrPt = dynamic_cast<TH1D*>(hPidCorrPt->Clone(Form("%s_SecCorr_Pt", fParticle.Data())));
 185                 }
 186
 187                 TH1D* hEffCorrPt = 0;
 188                 if(fEfficiency)
 189                 {
 190                         if(fDebugLevel>1) std::cout << "\t\tefficiency correction" << std::endl;
 191
 192                         if(fMCtoINEL)
 193                         {
 194                                 TH1D* hStatsMC = FindObj<TH1D>(fcorr, fCorrTag, species + "_Stats");
 195                                 fVtxFactor = this->GetVertexCorrection(hStats, hStatsMC);
 196                         }
 197
 198                         hEffCorrPt = this->Efficiency(hSecCorrPt, hEffVtxPt, hEffAccTrkPt, fParticle + "_EffCorr_Pt");
 199                 }
 200                 else
 201                 {
 202                         hEffCorrPt = dynamic_cast<TH1D*>(hSecCorrPt->Clone(Form("%s_EffCorr_Pt", fParticle.Data())));
 203                 }
 204
 205                 // corrected pt
 206
 207                 hCorrPt = dynamic_cast<TH1D*>(hEffCorrPt->Clone(Form("%s_Corr_Pt", fParticle.Data())));
 208
 209                 // debug histograms
 210
 211                 fdebug->cd(fOutputTag.Data());
 212
 213                 hOrigPt->Write(); // original distribution
 214                 hPidCorrPt->Write();
 215                 hSecCorrPt->Write();
 216                 hEffCorrPt->Write();
 217                 hCorrPt->Write();
 218
 219                 // free memory
 220                 delete hPidCorrPt;
 221                 delete hSecCorrPt;
 222                 delete hEffCorrPt;
 223         }
 224
 225         // save final pt
 226
 227         TH1D* hPt = dynamic_cast<TH1D*>(hCorrPt->Clone(Form("%s_Pt", fParticle.Data())));
 228         hPt->SetTitle(fParticle.Data());
 229         hPt->Reset();
 230         hPt->Sumw2();
 231
 232         Int_t lowbin = hPt->GetXaxis()->FindFixBin(fPtMin);
 233         Int_t hibin  = hPt->GetXaxis()->FindFixBin(fPtMax);
 234
 235         for(Int_t i=lowbin; i<hibin; ++i)
 236         {
 237                 hPt->SetBinContent(i, hCorrPt->GetBinContent(i));
 238                 hPt->SetBinError(i, hCorrPt->GetBinError(i));
 239         }
 240
 241         foutput->cd();
 242
 243         hPt->Write();
 244
 245         delete hPt;
 246         delete hCorrPt;
 247
 248         delete fdebug;
 249         delete fcorr;
 250
 251         // stats
 252
 253         if(fMakeStats)
 254         {
 255                 TH1D* hStatsFix = new TH1D(Form("%s_Stats", fParticle.Data()), "Stats", 6, 0, 6);
 256                 hStatsFix->GetXaxis()->SetBinLabel(1,"Events");
 257                 hStatsFix->GetXaxis()->SetBinLabel(2,"Trig");
 258                 hStatsFix->GetXaxis()->SetBinLabel(3,"AnaT");
 259                 hStatsFix->GetXaxis()->SetBinLabel(4,"Inel");
 260                 hStatsFix->GetXaxis()->SetBinLabel(5,"Vtx");
 261                 hStatsFix->GetXaxis()->SetBinLabel(6,"Vz");
 262                 hStatsFix->SetStats(0);
 263
 264                 // Fill
 265                 hStatsFix->Fill("Events",hStats->Integral(1,1));
 266                 hStatsFix->Fill("Trig",hStats->Integral(2,2));
 267                 if(fIsOnlyGen)
 268                 {
 269                         hStatsFix->Fill("Ana",hStats->Integral(2,2));
 270                         hStatsFix->Fill("Inel",hStats->Integral(2,2));
 271                 }
 272                 else
 273                 {
 274                         Double_t nTrig = hStats->Integral(2,2);
 275                         Double_t nAna  = hStats->Integral(3,3);
 276                         Double_t nVtx  = hStats->Integral(4,4);
 277                         Double_t nVz   = hStats->Integral(5,5);
 278                         Double_t nAnaT = nTrig; // all triggering events since we have the MC correction
 279                         Double_t nInel = nTrig/fTrigEff;
 280
 281                         if(!fMCtoINEL)
 282                         {
 283                                 // Add events without vertex to get the proportional trig. events
 284                                 Double_t nNoVtx = (nAna/nVtx)*(nTrig-nVtx);
 285                                 nAnaT = nAna + nNoVtx;
 286                                 nInel = nAnaT/fTrigEff;
 287                         }
 288
 289                         hStatsFix->Fill("AnaT",nAnaT);
 290                         hStatsFix->Fill("Inel",nInel);
 291                         hStatsFix->Fill("Vtx",nVtx);
 292                         hStatsFix->Fill("Vz",nVz);
 293                 }
 294
 295                 foutput->cd();
 296                 hStatsFix->Write();
 297                 delete hStatsFix;
 298         }
 299
 300         delete foutput;
 301         delete finput;
 302
 303         return 0;
 304 }
 305
 306 Double_t AliLnPt::GetVertexCorrection(const TH1D* hData, const TH1D* hMC) const
 307 {
 308 //
 309 // vertex correction factor
 310 //
 311         Double_t nAna   = hData->Integral(3,3);
 312         Double_t nVtx   = hData->Integral(4,4);
 313         Double_t nAnaMC = hMC->Integral(3,3);
 314         Double_t nVtxMC = hMC->Integral(4,4);
 315
 316         return (nVtx/nAna)/(nVtxMC/nAnaMC);
 317 }
 318
 319 Double_t AliLnPt::GetM2Width(Double_t pt, Double_t m) const
 320 {
 321 //
 322 // expected M2 width parameterization from LHC12a5b
 323 //
 324         Double_t c0 = 1.90212e-03;
 325         Double_t c1 = 1.29814e-02;
 326
 327         TString species = fParticle;
 328         species.ReplaceAll("Anti","");
 329
 330         if(species == "Deuteron")
 331         {
 332                 c0 = 1.37785e-02;
 333                 c1 = 6.78951e-02;
 334         }
 335         else if(species == "Triton")
 336         {
 337                 c0 = 3.79471e-02;
 338                 c1 = 1.87478e-01;
 339         }
 340         else if(species == "He3")
 341         {
 342                 c0 = 2.05887e-01;
 343                 c1 = 1.17813e-01;
 344         }
 345         else if(species == "Alpha")
 346         {
 347                 c0 = 4.13329e-01;
 348                 c1 = 2.27076e-01;
 349         }
 350
 351         return c0/(pt*pt) + c1*(pt*pt/(m*m) + 1.);
 352 }
 353
 354 Double_t AliLnPt::GetExpectedDT(Double_t pt, Double_t m) const
 355 {
 356 //
 357 // estimate of the expected time of flight (ns) for mass hypothesis m
 358 //
 359         Double_t l = 370.; // cm
 360         Double_t c = 2.99792458e+1; // cm/ns
 361
 362         return (l/c)*TMath::Sqrt(1.+(m/pt)*(m/pt));
 363 }
 364
 365 RooWorkspace* AliLnPt::GetToFModel(Double_t pt, const TString& name) const
 366 {
 367 //
 368 // model for the time of flight distribution
 369 //
 370         using namespace RooFit;
 371
 372         RooWorkspace* w = new RooWorkspace(name.Data());
 373
 374         // signal
 375         w->factory(Form("AliLnGaussianExpTail::Sd(x[%f,%f], mu[0,-0.06,0.06], sigma[0.120,0.05,0.3], tau[0.1,0.08,0.4])", fBkgMin, fBkgMax));
 376
 377         // background
 378         Double_t expTd   = this->GetExpectedDT(pt, GetMass("deuteron"));
 379         Double_t expDTp  = this->GetExpectedDT(pt, GetMass("proton")) - expTd;
 380         Double_t expDTk  = this->GetExpectedDT(pt, GetMass("kaon")) - expTd;
 381         Double_t expDTpi = this->GetExpectedDT(pt, GetMass("pion")) - expTd;
 382
 383         w->factory(Form("AliLnGaussianExpTail::ProtonBkg(x, muP[%f,%f-0.2,%f+0.2], widthP[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTp, expDTp, expDTp));
 384         w->factory(Form("AliLnGaussianExpTail::KaonBkg(x, muK[%f,%f-0.2,%f+0.2], widthK[0.120,0.05,0.3], tauP[0.08,0.01,0.2])",expDTk, expDTk, expDTk));
 385         w->factory(Form("AliLnGaussianExpTail::PionBkg(x, muPi[%f,%f-0.2,%f+0.2], widthPi[0.120,0.05,0.3], tauPi[0.08,0.01,0.2])",expDTpi, expDTpi, expDTpi));
 386
 387         w->factory("SUM::Bkg( Npi[0.3,0,1.e+9]*PionBkg, Nk[0.24,0,1.e+9]*KaonBkg, Np[0.47,0,1.e+9]*ProtonBkg, Np[0.47,0,1.e+9]*ProtonBkg)");
 388
 389         w->factory("SUM::model( Nd[0,1.e+9]*Sd,  Nbkg[0,1.e+9]*Bkg )");
 390
 391         w->var("x")->SetTitle("t - <t> (ns)");
 392
 393         return w;
 394 }
 395
 396 RooWorkspace* AliLnPt::GetM2Model(Double_t pt, Double_t m, const TString& name, Double_t max) const
 397 {
 398 //
 399 // mass squared model for deuterons
 400 //
 401         using namespace RooFit;
 402
 403         RooWorkspace* w = new RooWorkspace(name.Data());
 404
 405         // background
 406
 407         w->factory(Form("Exponential::PionBkg(x[%f,%f],lambdaPi[-0.3,-10.,0.])", fBkgMin, fBkgMax));
 408         w->factory("Exponential::KaonBkg(x,lambdaK[-0.4,-10.,0.])");
 409         w->factory("Exponential::ProtonBkg(x,lambdaP[-0.5,-10.,0.])");
 410
 411         w->factory("SUM::Bkg( Npi[0.3,0.,1.]*PionBkg, Nk[0.24,0.,1.]*KaonBkg, Np[0.47,0.,1.]*ProtonBkg )");
 412
 413         // signal
 414
 415         Double_t sigma = GetM2Width(pt, m);
 416         Double_t mm = (fPidProc == kMassSquared) ? m*m : 0;
 417
 418         if(pt<2.0)
 419         {
 420                 w->factory(Form("AliLnM2::Sd(x,mu[%f,%f,%f], sigma[%f,%f,%f],tau[0.1,0.05,0.4],p[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1, 1.16*pt,pt,1.33*pt));
 421         }
 422         else // neglect non-gaussian tails
 423         {
 424                 w->factory(Form("Gaussian::Sd(x, mu[%f,%f,%f], sigma[%f,%f,%f])",mm, mm-0.15, mm+0.1, sigma, 0.1, sigma+0.1));
 425         }
 426
 427         // model
 428         w->factory(Form("SUM::model( Nd[0.,%f]*Sd, Nbkg[0.,%f]*Bkg )", max, max));
 429
 430         w->var("x")->SetTitle("#it{m}^{2} (GeV^{2}/c^{4})");
 431
 432         return w;
 433 }
 434
 435 void AliLnPt::GetPtFromPid(TH1D* hPt, Double_t ptmin, Double_t ptmax, const TH2D* hPidPt, Double_t pidMin, Double_t pidMax) const
 436 {
 437 //
 438 // Integrate pid distribution
 439 //
 440         Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
 441         Int_t hibin  = hPt->GetXaxis()->FindFixBin(ptmax);
 442
 443         for(Int_t i=lowbin; i<hibin; ++i)
 444         {
 445                 TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
 446
 447                 Int_t pidLowBin = hPid->GetXaxis()->FindFixBin(pidMin);
 448                 Int_t pidHiBin = hPid->GetXaxis()->FindFixBin(pidMax);
 449
 450                 Double_t err = 0;
 451                 Double_t n = hPid->IntegralAndError(pidLowBin, pidHiBin, err);
 452
 453                 hPt->SetBinContent(i, n);
 454                 hPt->SetBinError(i, err);
 455
 456                 delete hPid;
 457         }
 458 }
 459
 460 TH1D* AliLnPt::PID(const TH1D* hPt, const TH2D* hPidPt2, Double_t ptmin, Double_t ptmax, const TString& name)
 461 {
 462 //
 463 // Remove contamination in the pt distribution
 464 //
 465         TH1D* hPidCorrPt = dynamic_cast<TH1D*>(hPt->Clone(name.Data()));
 466         if(hPidCorrPt == 0) return 0;
 467
 468         hPidCorrPt->Reset();
 469         hPidCorrPt->Sumw2();
 470
 471         // set binning
 472         const Int_t kNBin = 2;
 473
 474         TH2D* hPidPt = dynamic_cast<TH2D*>(hPidPt2->Clone("_pid_pt_"));
 475         hPidPt->Rebin2D(1,kNBin);
 476
 477         // integrate around the mean value
 478
 479         this->GetPtFromPid(hPidCorrPt, 0, ptmin, hPidPt, fIntMin, fIntMax);
 480
 481         // slice the m2-pt
 482
 483         using namespace RooFit;
 484
 485         Int_t lowbin = hPt->GetXaxis()->FindFixBin(ptmin);
 486         Int_t hibin  = hPt->GetXaxis()->FindFixBin(ptmax);
 487
 488         Double_t m = GetMass(fParticle);
 489
 490         for(Int_t i=lowbin; i<hibin; ++i)
 491         {
 492                 TH1D* hPid = hPidPt->ProjectionY(Form("%s_PID_%02d",fParticle.Data(),i),i,i);
 493
 494                 RooWorkspace* w = (fPidProc==kTimeOfFlight) ? GetToFModel( hPt->GetBinCenter(i), Form("%s_M2_%02d",fParticle.Data(),i)) : GetM2Model(hPt->GetBinCenter(i), m, Form("%s_M2_%02d",fParticle.Data(),i),hPid->Integral());
 495
 496                 RooRealVar* x = w->var("x");
 497                 RooDataHist data("data", "data to fit", *x, hPid);
 498
 499                 w->import(data);
 500
 501                 w->pdf("model")->fitTo(data, Minos(kTRUE), PrintLevel(-1), Verbose(kFALSE), Warnings(kFALSE), PrintEvalErrors(-1) );
 502
 503                 hPidCorrPt->SetBinContent(i, w->var("Nd")->getVal());
 504                 hPidCorrPt->SetBinError(i, w->var("Nd")->getError());
 505
 506                 w->Write();
 507
 508                 delete w;
 509                 delete hPid;
 510         }
 511
 512         delete hPidPt;
 513         return hPidCorrPt;
 514 }
 515
 516 TH1D* AliLnPt::Secondaries(const TH1D* hPt, const TH1D* hFracMatPt, const TH1D* hFracFdwnPt, const TString& name) const
 517 {
 518 //
 519 // remove contamination of secondaries using fractions
 520 // N_prim = N/(1+fmat+fdwn)
 521 //
 522         TF1* one = new TF1("one","1",0,10);
 523
 524         TH1D* xfactor = (TH1D*)hFracMatPt->Clone(Form("%s_1_fmat_fdwn_",fParticle.Data()));
 525         if(fFdwnCorr) xfactor->Add(hFracFdwnPt);
 526         xfactor->Add(one);
 527
 528         TH1D* hSecCorPt = (TH1D*)hPt->Clone(name.Data());
 529
 530         hSecCorPt->Divide(xfactor);
 531
 532         delete xfactor;
 533         delete one;
 534
 535         return hSecCorPt;
 536 }
 537
 538 TH1D* AliLnPt::Secondaries(const TH1D* hPt, TF1* fncMatPt, TF1* fncFdwnPt, const TString& name) const
 539 {
 540 //
 541 // remove contamination of secondaries using a function fitted to the fractions
 542 // N_prim = N/(1+fmat+fdwn)
 543 //
 544         TF1* fnc = new TF1("_secondaries_","1 + [0]*exp(-[1]*x)+[2] + [3]*exp(-[4]*x)+[5]", 0, 10);
 545
 546         Double_t param[6] = {0};
 547         fncMatPt->GetParameters(param);
 548         if(fFdwnCorr) fncFdwnPt->GetParameters(&param[3]);
 549
 550         Double_t err[6] = {0};
 551         for(Int_t i=0; i<3; ++i) err[i] = fncMatPt->GetParError(i);
 552         if(fFdwnCorr) for(Int_t i=3; i<6; ++i) err[i] = fncFdwnPt->GetParError(i);
 553
 554         fnc->SetParameters(param);
 555         fnc->SetParErrors(err);
 556
 557         TH1D* hSecCorPt = (TH1D*)hPt->Clone(name.Data());
 558         hSecCorPt->Divide(fnc);
 559
 560         delete fnc;
 561
 562         return hSecCorPt;
 563 }
 564
 565 TH1D* AliLnPt::Efficiency(const TH1D* hPt, const TH1D* hEffVtxPt, const TH1D* hEffAccTrkPt, const TString& name) const
 566 {
 567 //
 568 // correct by efficiency
 569 //
 570         TH1D* hEffCorPt = (TH1D*)hPt->Clone(name.Data());
 571         hEffCorPt->Divide(hEffAccTrkPt);
 572
 573         if(fMCtoINEL)
 574         {
 575                 hEffCorPt->Divide(hEffVtxPt);
 576                 hEffCorPt->Scale(fVtxFactor);
 577         }
 578
 579         return hEffCorPt;
 580 }