PWGLF/SPECTRA/Nuclei/B2/AliLnSecondaries.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 // removal of secondaries using DCA distributions
  17 // author: Eulogio Serradilla <eulogio.serradilla@cern.ch>
  18
  19 #include <Riostream.h>
  20 #include <TFile.h>
  21 #include <TString.h>
  22 #include <TF1.h>
  23 #include <TMath.h>
  24 #include <TFractionFitter.h>
  25 #include <TObjArray.h>
  26 #include <TH1D.h>
  27 #include <TH2D.h>
  28 #include <TBackCompFitter.h>
  29
  30 #include "AliLnSecondaries.h"
  31 #include "B2.h"
  32
  33 ClassImp(AliLnSecondaries)
  34
  35 AliLnSecondaries::AliLnSecondaries(const TString& particle, const TString& dataFilename, const TString& simuFilename, const TString& outputFilename, const TString& otag)
  36 : TObject()
  37 , fParticle(particle)
  38 , fDataFilename(dataFilename)
  39 , fSimuFilename(simuFilename)
  40 , fOutputFilename(outputFilename)
  41 , fOutputTag(otag)
  42 , fPtMin(0.5)
  43 , fPtMax(3.0)
  44 , fNbin(1)
  45 , fMinDCAxy(-1.5)
  46 , fMaxDCAxy(1.5)
  47 , fFracProc(0)
  48 , fMatDCAxyMod(AliLnSecondaries::kFlatDCAxy)
  49 , fANucTemplate(0)
  50 , fScMat(1)
  51 , fScFd(1)
  52 , fAddFakeTracks(1)
  53 {
  54 //
  55 // constructor
  56 //
  57         TH1::SetDefaultSumw2(); // switch on histogram errors
  58 }
  59
  60 AliLnSecondaries::~AliLnSecondaries()
  61 {
  62 //
  63 // destructor
  64 //
  65 }
  66
  67 Int_t AliLnSecondaries::Exec()
  68 {
  69 //
  70 // Fractions of primaries and secondaries using templates from simulations
  71 // Secondaries are from feed-down and/or from materials
  72 //
  73         using namespace std;
  74
  75         TFile* fdata = new TFile(fDataFilename.Data(),"read");
  76         if(fdata->IsZombie()) exit(1);
  77
  78         TFile* fsimu = new TFile(fSimuFilename.Data(),"read");
  79         if(fsimu->IsZombie()) exit(1);
  80
  81         TFile* fdebug =  new TFile(Form("debug-%s",fOutputFilename.Data()),"recreate");
  82
  83         if(fOutputTag != "") fdebug->mkdir(fOutputTag.Data());
  84         if(fOutputTag != "") fdebug->cd(fOutputTag.Data());
  85
  86         // --------- ideal values: primaries + no secondaries --------
  87
  88         TH1D* hPidPt = FindObj<TH1D>(fdata, fParticle + "_PID_Pt");
  89         hPidPt->Write();
  90
  91         const char* contrib[] = { "prim", "mat", "fdwn" };
  92
  93         TH1D* hFracPt[3];
  94         for(Int_t i=0; i<3; ++i)
  95         {
  96                 hFracPt[i] = this->ZeroClone(hPidPt, Form("%s_P%s_Pt",fParticle.Data(),contrib[i]));
  97                 hFracPt[i]->SetYTitle(Form("P_{%s}",contrib[i]));
  98         }
  99
 100         // fill
 101
 102         if(fFracProc == kMonteCarlo)
 103         {
 104                 // compute the fractions from the montecarlo simulation
 105                 TH1D* hAll  = FindObj<TH1D>(fsimu, fParticle + "_Sim_Pt");
 106                 TH1D* hPrim = FindObj<TH1D>(fsimu, fParticle + "_Sim_Prim_Pt");
 107                 TH1D* hMat  = FindObj<TH1D>(fsimu, fParticle + "_Sim_Mat_Pt");
 108                 TH1D* hFdwn = FindObj<TH1D>(fsimu, fParticle + "_Sim_Fdwn_Pt");
 109
 110                 delete hFracPt[0];
 111                 delete hFracPt[1];
 112                 delete hFracPt[2];
 113
 114                 hFracPt[0] = Divide(hPrim, hAll, fParticle + "_Pprim_Pt");
 115                 hFracPt[1] = Divide(hMat,  hAll, fParticle + "_Pmat_Pt");
 116                 hFracPt[2] = Divide(hFdwn, hAll, fParticle + "_Pfdwn_Pt");
 117
 118                 hFracPt[1]->Scale(fScMat);
 119                 hFracPt[2]->Scale(fScFd);
 120         }
 121         else if(fParticle == "AntiProton")
 122         {
 123                 // assume only secondaries from feed-down
 124                 TH2D* hDataDCAxyPt     = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
 125                 TH2D* hDataMCDCAxyPt   = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
 126                 TH2D* hPrimDCAxyPt     = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
 127                 TH2D* hFdwnDCAxyPt     = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
 128                 TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
 129                 TH2D* hFakeFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
 130
 131                 if(fAddFakeTracks)
 132                 {
 133                         hPrimDCAxyPt->Add(hFakePrimDCAxyPt);
 134                         hFdwnDCAxyPt->Add(hFakeFdwnDCAxyPt);
 135                 }
 136
 137                 hDataDCAxyPt->Rebin2D(1,fNbin);
 138                 hDataMCDCAxyPt->Rebin2D(1,fNbin);
 139                 hPrimDCAxyPt->Rebin2D(1,fNbin);
 140                 hFdwnDCAxyPt->Rebin2D(1,fNbin);
 141
 142                 this->GetFraction(hFracPt[0], hFracPt[2], hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hFdwnDCAxyPt, "Fdwn");
 143         }
 144         else if(fParticle.Contains("Anti"))
 145         {
 146                 // assume there are no secondaries for antinuclei
 147                 this->GetFraction(hFracPt[0]);
 148         }
 149         else if(fParticle == "Proton")
 150         {
 151                 // secondaries from materials and feed-down
 152                 TH2D* hDataDCAxyPt     = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
 153                 TH2D* hDataMCDCAxyPt   = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
 154                 TH2D* hPrimDCAxyPt     = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
 155                 TH2D* hMatDCAxyPt      = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
 156                 TH2D* hFdwnDCAxyPt     = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fdwn_DCAxy_Pt");
 157                 TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
 158                 TH2D* hFakeMatDCAxyPt  = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Mat_DCAxy_Pt");
 159                 TH2D* hFakeFdwnDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Fdwn_DCAxy_Pt");
 160
 161                 if(fAddFakeTracks)
 162                 {
 163                         hPrimDCAxyPt->Add(hFakePrimDCAxyPt);
 164                         hMatDCAxyPt->Add(hFakeMatDCAxyPt);
 165                         hFdwnDCAxyPt->Add(hFakeFdwnDCAxyPt);
 166                 }
 167
 168                 if(fMatDCAxyMod == kFlatDCAxy)
 169                 {
 170                         hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
 171                 }
 172
 173                 hDataDCAxyPt->Rebin2D(1,fNbin);
 174                 hDataMCDCAxyPt->Rebin2D(1,fNbin);
 175                 hPrimDCAxyPt->Rebin2D(1,fNbin);
 176                 hMatDCAxyPt->Rebin2D(1,fNbin);
 177                 hFdwnDCAxyPt->Rebin2D(1,fNbin);
 178
 179                 this->GetFraction(hFracPt, hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hMatDCAxyPt, hFdwnDCAxyPt);
 180         }
 181         else
 182         {
 183                 // only secondaries from materials for nuclei
 184                 TH2D* hDataDCAxyPt     = FindObj<TH2D>(fdata, fParticle + "_PID_DCAxy_Pt");
 185                 TH2D* hDataMCDCAxyPt   = FindObj<TH2D>(fsimu, fParticle + "_PID_DCAxy_Pt");
 186                 TH2D* hPrimDCAxyPt     = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Prim_DCAxy_Pt");
 187                 TH2D* hFakePrimDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Fake_Prim_DCAxy_Pt");
 188
 189                 if(fAddFakeTracks) hPrimDCAxyPt->Add(hFakePrimDCAxyPt);
 190
 191                 if(fANucTemplate)
 192                 {
 193                         hPrimDCAxyPt   = FindObj<TH2D>(fdata, Form("Anti%s_PID_DCAxy_Pt",fParticle.Data()));
 194                 }
 195
 196                 TH2D* hMatDCAxyPt = FindObj<TH2D>(fsimu, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
 197                 if(fMatDCAxyMod == kFlatDCAxy)
 198                 {
 199                         hMatDCAxyPt = this->GetFlatDCAxyPt(10, hPrimDCAxyPt, fParticle + "_Sim_PID_Mat_DCAxy_Pt");
 200                 }
 201
 202                 hDataDCAxyPt->Rebin2D(1,fNbin);
 203                 hDataMCDCAxyPt->Rebin2D(1,fNbin);
 204                 hPrimDCAxyPt->Rebin2D(1,fNbin);
 205                 hMatDCAxyPt->Rebin2D(1,fNbin);
 206
 207                 this->GetFraction(hFracPt[0], hFracPt[1], hDataDCAxyPt, hDataMCDCAxyPt, hPrimDCAxyPt, hMatDCAxyPt, "Mat");
 208         }
 209
 210         // --------- save --------------
 211
 212         TFile* foutput = new TFile(fOutputFilename.Data(),"recreate");
 213         if(fOutputTag != "")
 214         {
 215                 foutput->mkdir(fOutputTag.Data());
 216                 foutput->cd(fOutputTag.Data());
 217         }
 218
 219         // ----- fractions w.r.t. primaries ------------
 220
 221         TH1D* hMatFracPt = (TH1D*)hFracPt[1]->Clone(Form("%s_Frac_Mat_Pt", fParticle.Data()));
 222         hMatFracPt->SetYTitle("P_{1} / P_{0}");
 223         hMatFracPt->Divide(hFracPt[0]);
 224         hMatFracPt->Write();
 225
 226         TH1D* hFdwnFracPt = (TH1D*)hFracPt[2]->Clone( Form("%s_Frac_Fdwn_Pt",fParticle.Data()));
 227         hFdwnFracPt->SetYTitle("P_{2} / P_{0}");
 228         hFdwnFracPt->Divide(hFracPt[0]);
 229         hFdwnFracPt->Write();
 230
 231         // --------- optionally fit fractions to extrapolate at high pt -----------
 232
 233         TF1* fncMat = this->GetMatFraction(fParticle + "_Frac_Mat_Fit_Pt");
 234         if(fParticle.Contains("Anti"))
 235         {
 236                 fncMat->SetParameters(0,0,0);
 237         }
 238         else if(fParticle=="Proton")
 239         {
 240                 hMatFracPt->Fit(fncMat, "NQ", "", 0.,2);
 241         }
 242         else
 243         {
 244                 hMatFracPt->Fit(fncMat, "NQ", "", 0.,1.6);
 245         }
 246
 247         TF1* fncFdwn = this->GetFdwnFraction(fParticle + "_Frac_Fdwn_Fit_Pt");
 248         if(!fParticle.Contains("Proton"))
 249         {
 250                 fncFdwn->SetParameters(0,0,0);
 251         }
 252         else
 253         {
 254                 hFdwnFracPt->Fit(fncFdwn, "NQ", "", 0.5,3.);
 255         }
 256
 257         fncFdwn->Write();
 258         fncMat->Write();
 259
 260         // ---------- clean ---------
 261
 262         delete fncFdwn;
 263         delete fncMat;
 264
 265         delete hMatFracPt;
 266         delete hFdwnFracPt;
 267
 268         for(Int_t i=0; i<3; ++i)
 269         {
 270                 hFracPt[i]->Write();
 271                 delete hFracPt[i];
 272         }
 273
 274         delete foutput;
 275         delete fdebug;
 276         delete fsimu;
 277         delete fdata;
 278
 279         delete dynamic_cast<TBackCompFitter*>(TVirtualFitter::GetFitter());
 280
 281         return 0;
 282 }
 283
 284 void AliLnSecondaries::GetFraction(TH1D* hPrimPt) const
 285 {
 286 //
 287 // No secondaries only primaries
 288 //
 289         TF1* one = new TF1("one", "1", hPrimPt->GetXaxis()->GetXmin(), hPrimPt->GetXaxis()->GetXmax());
 290         hPrimPt->Reset();
 291         hPrimPt->Sumw2();
 292         hPrimPt->Add(one);
 293         delete one;
 294 }
 295
 296 void AliLnSecondaries::GetFraction(TH1D* hPrimPt, TH1D* hSecPt, const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hSecDCAxyPt, const TString& sec) const
 297 {
 298 //
 299 // slice the DCA distributions and compute the fractions
 300 // 2 contributions (primaries and secondaries)
 301 //
 302         TString nosec = (sec == "Fdwn") ? "Mat" : "Fdwn";
 303
 304         Int_t lowbin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMin);
 305         Int_t hibin  = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax);
 306
 307         for(Int_t i=lowbin; i<hibin; ++i)
 308         {
 309                 TH1D* hDCAxy      = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
 310                 TH1D* hMCDCAxy    = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i);
 311                 TH1D* hPrimDCAxy  = hPrimDCAxyPt->ProjectionY(Form("%s_Prim_DCAxy_%02d",fParticle.Data(),i),i,i);
 312                 TH1D* hSecDCAxy   = hSecDCAxyPt->ProjectionY(Form("%s_%s_DCAxy_%02d",fParticle.Data(),sec.Data(),i),i,i);
 313
 314                 // fractions
 315                 Double_t frac[2] = {0};
 316                 Double_t err[2]  = {0};
 317
 318                 this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hSecDCAxy, i, sec);
 319
 320                 // fill
 321                 TH1D* hFracPt[] = { hPrimPt, hSecPt};
 322                 for(Int_t j=0; j<2; ++j)
 323                 {
 324                         hFracPt[j]->SetBinContent(i, frac[j]);
 325                         hFracPt[j]->SetBinError(i, err[j]);
 326                 }
 327
 328                 // ------------- debug histograms ------------
 329
 330                 hDCAxy->Write();
 331                 hMCDCAxy->Write();
 332                 hPrimDCAxy->Write();
 333                 hSecDCAxy->Write();
 334
 335                 TH1D* hNoSecDCAxy = this->ZeroClone(hMCDCAxy, Form("%s_%s_DCAxy_%02d",fParticle.Data(),nosec.Data(),i));
 336                 TH1D* hNoSecFit   = this->ZeroClone(hMCDCAxy, Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),nosec.Data(),i));
 337
 338                 hNoSecDCAxy->Write();
 339                 hNoSecFit->Write();
 340
 341                 delete hNoSecDCAxy;
 342                 delete hNoSecFit;
 343
 344                 // ----------------------- end debug ----------------------
 345
 346                 delete hDCAxy;
 347                 delete hMCDCAxy;
 348                 delete hPrimDCAxy;
 349                 delete hSecDCAxy;
 350         }
 351 }
 352
 353 void AliLnSecondaries::GetFraction(TH1D* hFracPt[3], const TH2D* hDCAxyPt, const TH2D* hMCDCAxyPt, const TH2D* hPrimDCAxyPt, const TH2D* hMatDCAxyPt, const TH2D* hFdwnDCAxyPt) const
 354 {
 355 //
 356 // slice the DCA distribution and get the fractions for each pt bin
 357 // (3 contributions)
 358 //
 359         Int_t lowbin = hDCAxyPt->GetXaxis()->FindFixBin(fPtMin);
 360         Int_t hibin  = hDCAxyPt->GetXaxis()->FindFixBin(fPtMax);
 361
 362         for(Int_t i=lowbin; i<hibin; ++i)
 363         {
 364                 // slices
 365                 TH1D* hDCAxy     = hDCAxyPt->ProjectionY(Form("%s_Data_DCAxy_%02d",fParticle.Data(),i),i,i);
 366                 TH1D* hMCDCAxy   = hMCDCAxyPt->ProjectionY(Form("%s_SimData_DCAxy_%02d",fParticle.Data(),i),i,i);
 367                 TH1D* hPrimDCAxy = hPrimDCAxyPt->ProjectionY(Form("%s_Prim_DCAxy_%02d",fParticle.Data(),i),i,i);
 368                 TH1D* hMatDCAxy  = hMatDCAxyPt->ProjectionY(Form("%s_Mat_DCAxy_%02d",fParticle.Data(),i),i,i);
 369                 TH1D* hFdwnDCAxy = hFdwnDCAxyPt->ProjectionY(Form("%s_Fdwn_DCAxy_%02d",fParticle.Data(),i),i,i);
 370
 371                 // fractions
 372                 Double_t frac[3] = {0};
 373                 Double_t err[3]  = {0};
 374
 375                 this->GetTFFfractions(frac, err, hDCAxy, hPrimDCAxy, hMatDCAxy, hFdwnDCAxy, i);
 376
 377                 // fill
 378                 for(Int_t k=0; k<3; ++k)
 379                 {
 380                         hFracPt[k]->SetBinContent(i, frac[k]);
 381                         hFracPt[k]->SetBinError(i, err[k]);
 382                 }
 383
 384                 // -------------- debug ---------------------------
 385
 386                 hDCAxy->Write();
 387                 hMCDCAxy->Write();
 388                 hPrimDCAxy->Write();
 389                 hMatDCAxy->Write();
 390                 hFdwnDCAxy->Write();
 391
 392                 // --------------------- end debug ----------------
 393
 394                 delete hDCAxy;
 395                 delete hMCDCAxy;
 396                 delete hPrimDCAxy;
 397                 delete hMatDCAxy;
 398                 delete hFdwnDCAxy;
 399         }
 400 }
 401
 402 Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hSec, Int_t ibin, const TString& sec) const
 403 {
 404 //
 405 // find the fractions of 2 contributions
 406 //
 407         TObjArray* mc = new TObjArray(2);
 408
 409         mc->Add(hPrim);
 410         mc->Add(hSec);
 411
 412         TFractionFitter* fit = new TFractionFitter(hData, mc, "Q"); // initialise
 413         fit->Constrain(1,0.0001,1.);
 414         fit->Constrain(2,0.0001,1.);
 415
 416         Int_t status = fit->Fit();
 417         status = fit->Fit();
 418         status = fit->Fit();
 419
 420         if (status == 0) // get fractions
 421         {
 422                 fit->GetResult(0, frac[0], err[0]);
 423                 fit->GetResult(1, frac[1], err[1]);
 424         }
 425
 426         const char* contrib[] = {"Prim", sec.Data() };
 427         this->WriteTFFdebug(hData, fit, status, ibin, contrib, frac, 2);
 428
 429         delete mc;
 430         delete fit;
 431
 432         return status;
 433 }
 434
 435 Int_t AliLnSecondaries::GetTFFfractions(Double_t* frac, Double_t* err, TH1D* hData, TH1D* hPrim, TH1D* hMat, TH1D* hFdwn, Int_t ibin) const
 436 {
 437 //
 438 // find the fractions of 3 contributions
 439 //
 440         TObjArray* mc = new TObjArray(3);
 441
 442         mc->Add(hPrim);
 443         mc->Add(hMat);
 444         mc->Add(hFdwn);
 445
 446         TFractionFitter* fit = new TFractionFitter(hData, mc, "Q");
 447         fit->Constrain(1,0.0001,1.);
 448         fit->Constrain(2,0.0001,1.);
 449         fit->Constrain(3,0.0001,1.);
 450         //fit->SetRangeX(1,50);
 451
 452         Int_t status = fit->Fit();
 453         status = fit->Fit();
 454         status = fit->Fit();
 455
 456         if (status == 0) // get fractions
 457         {
 458                 fit->GetResult(0, frac[0], err[0]);
 459                 fit->GetResult(1, frac[1], err[1]);
 460                 fit->GetResult(2, frac[2], err[2]);
 461         }
 462
 463         const char* contrib[] = {"Prim", "Mat", "Fdwn"};
 464         this->WriteTFFdebug(hData, fit, status, ibin, contrib, frac, 3);
 465
 466         delete mc;
 467         delete fit;
 468
 469         return status;
 470 }
 471
 472 void AliLnSecondaries::WriteTFFdebug(const TH1D* hData, TFractionFitter* fit, Int_t status, Int_t ibin, const char* contrib[], Double_t* frac, Int_t kmax) const
 473 {
 474 //
 475 // Write TFractionFitter debug histograms
 476 //
 477         if (status == 0)
 478         {
 479                 TH1D* hResult = (TH1D*) fit->GetPlot();
 480
 481                 hResult->SetName(Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
 482                 hResult->Write();
 483
 484                 for(Int_t k=0; k<kmax; ++k)
 485                 {
 486                         TH1D* hMCpred = (TH1D*) fit->GetMCPrediction(k);
 487                         hMCpred->SetName(Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),contrib[k],ibin));
 488
 489                         hMCpred->Sumw2();
 490                         hMCpred->Scale(frac[k]*hResult->Integral()/hMCpred->Integral());
 491
 492                         hMCpred->Write();
 493                 }
 494         }
 495         else // write void histograms
 496         {
 497                 TH1D* hZero = this->ZeroClone(hData,Form("%s_Fit_Data_DCAxy_%02d",fParticle.Data(),ibin));
 498                 hZero->Write();
 499
 500                 for(Int_t k=0; k<kmax; ++k)
 501                 {
 502                         hZero->SetName(Form("%s_Fit_%s_DCAxy_%02d",fParticle.Data(),contrib[k],ibin));
 503                         hZero->Write();
 504                 }
 505                 delete hZero;
 506         }
 507 }
 508
 509 TH1D* AliLnSecondaries::ZeroClone(const TH1D* h, const TString& name) const
 510 {
 511 //
 512 // clone histogram and reset to zero
 513 //
 514         TH1D* clone = dynamic_cast<TH1D*>(h->Clone(name.Data()));
 515         if(clone == 0)
 516         {
 517                 this->Error("ZeroClone", "could not clone %s", h->GetName());
 518                 exit(1);
 519         }
 520         clone->Reset();
 521         clone->Sumw2();
 522
 523         return clone;
 524 }
 525
 526 TH2D* AliLnSecondaries::GetFlatDCAxyPt(Double_t max, const TH2D* hDCAxyPt, const TString& name) const
 527 {
 528 //
 529 // generate a flat dcaxy-pt distribution
 530 //
 531         TF1* fnc = new TF1("gen","1", fMinDCAxy, fMaxDCAxy);
 532         TH2D* h = (TH2D*)hDCAxyPt->Clone(name.Data());
 533         h->Reset();
 534         h->Sumw2();
 535
 536         TAxis* xAxis = h->GetXaxis();
 537         TAxis* yAxis = h->GetYaxis();
 538
 539         for(Int_t i=1; i<xAxis->GetNbins()+1; ++i)
 540         {
 541                 Double_t pt = xAxis->GetBinCenter(i);
 542
 543                 for(Int_t j=1; j<yAxis->GetNbins()+1; ++j)
 544                 {
 545                         for(Int_t k=0; k<max; ++k) h->Fill(pt,fnc->GetRandom());
 546                 }
 547         }
 548
 549         delete fnc;
 550         return h;
 551 }
 552
 553 TF1* AliLnSecondaries::GetMatFraction(const TString& name) const
 554 {
 555 //
 556 // model the material fraction as an exponential + a constant
 557 //
 558         TF1* fnc = new TF1(name.Data(), "[0]*exp(-[1]*x)+[2]",0,10);
 559         fnc->SetParameters(1.7,5.5,0.01);
 560
 561         fnc->SetParLimits(2,0,0.02); // asymptotic value
 562
 563         return fnc;
 564 }
 565
 566 TF1* AliLnSecondaries::GetFdwnFraction(const TString& name) const
 567 {
 568 //
 569 // model the feed-down fraction as an exponential + constant
 570 //
 571         TF1* fnc = new TF1(name, "[0]*exp(-[1]*x)+[2]",0,10);
 572         fnc->SetParameters(0.5,1,0.02);
 573
 574         fnc->SetParLimits(2,0,0.3); // asymptotic value
 575
 576         return fnc;
 577 }