PWGPP/TPC/AliPerformanceEff.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceEff class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms/graphs) are stored in the folder which is
   8 // a data member of AliPerformanceEff.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 //------------------------------------------------------------------------------
  12
  13 /*
  14
  15   // after running comparison task, read the file, and get component
  16   gROOT->LoadMacro("$ALICE_ROOT/PWGPP/Macros/LoadMyLibs.C");
  17   LoadMyLibs();
  18   TFile f("Output.root");
  19   AliPerformanceEff * compObj = (AliPerformanceEff*)coutput->FindObject("AliPerformanceEff");
  20
  21   // Analyse comparison data
  22   compObj->Analyse();
  23
  24   // the output histograms/graphs will be stored in the folder "folderEff"
  25   compObj->GetAnalysisFolder()->ls("*");
  26
  27   // user can save whole comparison object (or only folder with anlysed histograms)
  28   // in the seperate output file (e.g.)
  29   TFile fout("Analysed_Eff.root","recreate");
  30   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  31   fout.Close();
  32
  33 */
  34
  35 #include <TAxis.h>
  36 #include <TH1D.h>
  37 #include "THnSparse.h"
  38
  39 //
  40 #include "AliESDtrack.h"
  41 #include "AliRecInfoCuts.h"
  42 #include "AliMCInfoCuts.h"
  43 #include "AliLog.h"
  44 #include "AliESDVertex.h"
  45 #include "AliExternalTrackParam.h"
  46 #include "AliTracker.h"
  47 #include "AliESDEvent.h"
  48 #include "AliMCEvent.h"
  49 #include "AliMCParticle.h"
  50 #include "AliHeader.h"
  51 #include "AliGenEventHeader.h"
  52 #include "AliStack.h"
  53 #include "AliPerformanceEff.h"
  54
  55 using namespace std;
  56
  57 ClassImp(AliPerformanceEff)
  58
  59 //_____________________________________________________________________________
  60 AliPerformanceEff::AliPerformanceEff():
  61   AliPerformanceObject("AliPerformanceEff"),
  62
  63   // histograms
  64   fEffHisto(0),
  65   fEffSecHisto(0),
  66
  67   // Cuts
  68   fCutsRC(0),
  69   fCutsMC(0),
  70
  71   // histogram folder
  72   fAnalysisFolder(0)
  73 {
  74   // default consttructor
  75   Init();
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliPerformanceEff::AliPerformanceEff(Char_t* name="AliPerformanceEff",Char_t*title="AliPerformanceEff",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  80   AliPerformanceObject(name,title),
  81
  82   // histograms
  83   fEffHisto(0),
  84   fEffSecHisto(0),
  85
  86   // Cuts
  87   fCutsRC(0),
  88   fCutsMC(0),
  89
  90   // histogram folder
  91   fAnalysisFolder(0)
  92 {
  93   // named constructor
  94   //
  95   SetAnalysisMode(analysisMode);
  96   SetHptGenerator(hptGenerator);
  97
  98   Init();
  99 }
 100
 101
 102 //_____________________________________________________________________________
 103 AliPerformanceEff::~AliPerformanceEff()
 104 {
 105 // destructor
 106
 107   if(fEffHisto)  delete  fEffHisto; fEffHisto=0;
 108   if(fEffSecHisto)  delete  fEffSecHisto; fEffSecHisto=0;
 109   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 110 }
 111
 112 //_____________________________________________________________________________
 113 void AliPerformanceEff::Init()
 114 {
 115   // Init histograms
 116   //
 117
 118   // set pt bins
 119   Int_t nPtBins = 50;
 120   Double_t ptMin = 1.e-2, ptMax = 20.;
 121
 122   Double_t *binsPt = 0;
 123
 124   if (IsHptGenerator())  {
 125         ptMax = 100.;
 126   }
 127    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 128
 129   /*
 130   Int_t nPtBins = 31;
 131   Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 132   Double_t ptMin = 0., ptMax = 10.;
 133
 134   if(IsHptGenerator() == kTRUE) {
 135     nPtBins = 100;
 136     ptMin = 0.; ptMax = 100.;
 137   }
 138   */
 139
 140   //mceta:mcphi:mcpt:pid:recStatus:findable:charge
 141   Int_t binsEffHisto[9]={30,144,nPtBins,5,2,2,3,fgkMaxClones+1,fgkMaxFakes+1};
 142   Double_t minEffHisto[9]={-1.5,0.,ptMin,0.,0.,0.,-1.5,0,0};
 143   Double_t maxEffHisto[9]={ 1.5,2.*TMath::Pi(), ptMax,5.,2.,2.,1.5,fgkMaxClones,fgkMaxFakes};
 144
 145   fEffHisto = new THnSparseF("fEffHisto","mceta:mcphi:mcpt:pid:recStatus:findable:charge:nclones:nfakes",9,binsEffHisto,minEffHisto,maxEffHisto);
 146   fEffHisto->SetBinEdges(2,binsPt);
 147
 148   fEffHisto->GetAxis(0)->SetTitle("#eta_{mc}");
 149   fEffHisto->GetAxis(1)->SetTitle("#phi_{mc} (rad)");
 150   fEffHisto->GetAxis(2)->SetTitle("p_{Tmc} (GeV/c)");
 151   fEffHisto->GetAxis(3)->SetTitle("pid");
 152   fEffHisto->GetAxis(4)->SetTitle("recStatus");
 153   fEffHisto->GetAxis(5)->SetTitle("findable");
 154   fEffHisto->GetAxis(6)->SetTitle("charge");
 155   fEffHisto->GetAxis(7)->SetTitle("nClones");
 156   fEffHisto->GetAxis(8)->SetTitle("nFakes");
 157   fEffHisto->Sumw2();
 158
 159   //mceta:mcphi:mcpt:pid:recStatus:findable:mcR:mother_phi:mother_eta:charge
 160   Int_t binsEffSecHisto[12]={30,60,nPtBins,5,2,2,100,60,30,3,fgkMaxClones+1,fgkMaxFakes+1};
 161   Double_t minEffSecHisto[12]={-1.5,0.,ptMin,0.,0.,0.,0.,0.,-1.5,-1.5,0,0};
 162   Double_t maxEffSecHisto[12]={ 1.5,2.*TMath::Pi(), ptMax,5.,2.,2.,200,2.*TMath::Pi(),1.5,1.5,fgkMaxClones,fgkMaxFakes};
 163
 164   fEffSecHisto = new THnSparseF("fEffSecHisto","mceta:mcphi:mcpt:pid:recStatus:findable:mcR:mother_phi:mother_eta:charge:nclones:nfakes",12,binsEffSecHisto,minEffSecHisto,maxEffSecHisto);
 165   fEffSecHisto->SetBinEdges(2,binsPt);
 166
 167   fEffSecHisto->GetAxis(0)->SetTitle("#eta_{mc}");
 168   fEffSecHisto->GetAxis(1)->SetTitle("#phi_{mc} (rad)");
 169   fEffSecHisto->GetAxis(2)->SetTitle("p_{Tmc} (GeV/c)");
 170   fEffSecHisto->GetAxis(3)->SetTitle("pid");
 171   fEffSecHisto->GetAxis(4)->SetTitle("recStatus");
 172   fEffSecHisto->GetAxis(5)->SetTitle("findable");
 173   fEffSecHisto->GetAxis(6)->SetTitle("mcR (cm)");
 174   fEffSecHisto->GetAxis(7)->SetTitle("mother_phi (rad)");
 175   fEffSecHisto->GetAxis(8)->SetTitle("mother_eta");
 176   fEffSecHisto->GetAxis(9)->SetTitle("charge");
 177   fEffSecHisto->GetAxis(10)->SetTitle("nClones");
 178   fEffSecHisto->GetAxis(11)->SetTitle("nFakes");
 179   fEffSecHisto->Sumw2();
 180
 181   // init cuts
 182   if(!fCutsMC)
 183     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 184   if(!fCutsRC)
 185     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 186
 187   // init folder
 188   fAnalysisFolder = CreateFolder("folderEff","Analysis Efficiency Folder");
 189 }
 190
 191 //_____________________________________________________________________________
 192 void AliPerformanceEff::ProcessTPC(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
 193 {
 194   // Fill TPC only efficiency comparison information
 195   if(!esdEvent) return;
 196   if(!mcEvent) return;
 197
 198   AliStack *stack = mcEvent->Stack();
 199   if (!stack) {
 200     AliDebug(AliLog::kError, "Stack not available");
 201     return;
 202   }
 203
 204   Int_t *labelsRec =  NULL;
 205   labelsRec =  new Int_t[esdEvent->GetNumberOfTracks()];
 206   if(!labelsRec)
 207   {
 208      Printf("Cannot create labelsRec");
 209      return;
 210   }
 211   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRec[i] = 0; }
 212
 213   Int_t *labelsAllRec =  NULL;
 214   labelsAllRec =  new Int_t[esdEvent->GetNumberOfTracks()];
 215   if(!labelsAllRec) {
 216      delete  [] labelsRec;
 217      Printf("Cannot create labelsAllRec");
 218      return;
 219   }
 220   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRec[i] = 0; }
 221
 222   // loop over rec. tracks
 223   AliESDtrack *track=0;
 224   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 225   {
 226     track = esdEvent->GetTrack(iTrack);
 227     if(!track) continue;
 228     if(track->Charge()==0) continue;
 229
 230     // if not fUseKinkDaughters don't use tracks with kink index > 0
 231     if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
 232
 233     //Int_t label = TMath::Abs(track->GetLabel());
 234         Int_t label = track->GetTPCLabel(); //Use TPC-only label for TPC-only efficiency analysis
 235     labelsAllRec[iTrack]=label;
 236
 237     // TPC only
 238     if(IsRecTPC(track) != 0) {
 239       labelsRec[iTrack]=label;
 240     }
 241
 242   }
 243
 244   //
 245   // MC histograms for efficiency studies
 246   //
 247   Int_t nPart  = stack->GetNtrack();
 248   //Int_t nPart  = stack->GetNprimary();
 249   for (Int_t iMc = 0; iMc < nPart; ++iMc)
 250   {
 251         if (iMc == 0) continue;         //Cannot distinguish between track or fake track
 252     TParticle* particle = stack->Particle(iMc);
 253     if (!particle) continue;
 254     if (!particle->GetPDG()) continue;
 255     if (particle->GetPDG()->Charge() == 0.0) continue;
 256
 257     // physical primary
 258      Bool_t prim = stack->IsPhysicalPrimary(iMc);
 259      if(!prim) continue;
 260
 261     // --- check for double filling in stack
 262     // use only particles with no daughters in the list of primaries
 263     Int_t nDaughters = 0;// particle->GetNDaughters();
 264
 265     for( Int_t iDaught=0; iDaught < particle->GetNDaughters(); iDaught++ ) {
 266       if( particle->GetDaughter(iDaught) < stack->GetNprimary() )
 267         nDaughters++;
 268     }
 269
 270     if( nDaughters > 0 )
 271       continue;
 272     // --- check for double filling in stack
 273
 274     /*Bool_t findable = kFALSE;
 275     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 276     {
 277       // check findable
 278       if(iMc == labelsAllRec[iRec])
 279       {
 280         findable = IsFindable(mcEvent,iMc);
 281         break;
 282       }
 283     }*/
 284     Bool_t findable = IsFindable(mcEvent,iMc);
 285
 286     Bool_t recStatus = kFALSE;
 287     Int_t nClones = 0, nFakes = 0;
 288     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 289     {
 290       // check reconstructed
 291       if(iMc == labelsRec[iRec])
 292       {
 293         if (recStatus && nClones < fgkMaxClones) nClones++;
 294         recStatus = kTRUE;
 295       }
 296           //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
 297           if (labelsRec[iRec] < 0 && -labelsRec[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
 298     }
 299
 300     // Only 5 charged particle species (e,mu,pi,K,p)
 301     if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
 302
 303     // transform Pdg to Pid
 304     Int_t pid = TransformToPID(particle);
 305
 306     Float_t mceta =  particle->Eta();
 307     Float_t mcphi =  particle->Phi();
 308     if(mcphi<0) mcphi += 2.*TMath::Pi();
 309     Float_t mcpt = particle->Pt();
 310     Float_t charge = 0.;
 311     if (particle->GetPDG()->Charge() < 0)  charge = -1.;
 312     else if (particle->GetPDG()->Charge() > 0)  charge = 1.;
 313
 314     // Fill histograms
 315     Double_t vEffHisto[9] = {mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes};
 316     fEffHisto->Fill(vEffHisto);
 317   }
 318   if(labelsRec) delete [] labelsRec; labelsRec = 0;
 319   if(labelsAllRec) delete [] labelsAllRec; labelsAllRec = 0;
 320 }
 321
 322 //_____________________________________________________________________________
 323 void AliPerformanceEff::ProcessTPCSec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
 324 {
 325   // Fill TPC only efficiency comparison information for secondaries
 326
 327   if(!esdEvent) return;
 328
 329   Int_t *labelsRecSec =  NULL;
 330   labelsRecSec =  new Int_t[esdEvent->GetNumberOfTracks()];
 331   if(!labelsRecSec)
 332   {
 333      Printf("Cannot create labelsRecSec");
 334      return;
 335   }
 336   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecSec[i] = 0; }
 337
 338   Int_t *labelsAllRecSec =  NULL;
 339   labelsAllRecSec =  new Int_t[esdEvent->GetNumberOfTracks()];
 340   if(!labelsAllRecSec) {
 341      delete [] labelsRecSec;
 342      Printf("Cannot create labelsAllRecSec");
 343      return;
 344   }
 345   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecSec[i] = 0; }
 346
 347   // loop over rec. tracks
 348   AliESDtrack *track=0;
 349   Int_t multAll=0, multRec=0;
 350   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 351   {
 352     track = esdEvent->GetTrack(iTrack);
 353     if(!track) continue;
 354     if(track->Charge()==0) continue;
 355
 356     // if not fUseKinkDaughters don't use tracks with kink index > 0
 357     if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
 358
 359     //Int_t label = TMath::Abs(track->GetLabel());
 360         Int_t label = track->GetTPCLabel(); //Use TPC-only label for TPC-only efficiency analysis
 361     labelsAllRecSec[multAll]=label;
 362     multAll++;
 363
 364     // TPC only
 365     if(IsRecTPC(track) != 0) {
 366       labelsRecSec[multRec]=label;
 367       multRec++;
 368     }
 369   }
 370
 371   //
 372   // MC histograms for efficiency studies
 373   //
 374   if(mcEvent)  {
 375
 376   AliStack *stack = mcEvent->Stack();
 377   if (stack) {
 378
 379   Int_t nPart  = stack->GetNtrack();
 380   //Int_t nPart  = stack->GetNprimary();
 381   for (Int_t iMc = 0; iMc < nPart; ++iMc)
 382   {
 383         if (iMc == 0) continue;         //Cannot distinguish between track or fake track
 384     TParticle* particle = stack->Particle(iMc);
 385     if (!particle) continue;
 386     if (!particle->GetPDG()) continue;
 387     if (particle->GetPDG()->Charge() == 0.0) continue;
 388
 389     // physical primary
 390     Bool_t prim = stack->IsPhysicalPrimary(iMc);
 391
 392     // only secondaries which can be reconstructed at TPC
 393     if(prim) continue;
 394
 395     //Float_t radius = TMath::Sqrt(particle->Vx()*particle->Vx()+particle->Vy()*particle->Vy()+particle->Vz()*particle->Vz());
 396     //if(radius > fCutsMC->GetMaxR()) continue;
 397
 398     // only secondary electrons from gamma conversion
 399     //if( TMath::Abs(particle->GetPdgCode())!=fCutsMC->GetEM() ||   particle->GetUniqueID() != 5) continue;
 400
 401     /*Bool_t findable = kFALSE;
 402     for(Int_t iRec=0; iRec<multAll; ++iRec)
 403     {
 404       // check findable
 405       if(iMc == labelsAllRecSec[iRec])
 406       {
 407         findable = IsFindable(mcEvent,iMc);
 408         break;
 409       }
 410     }*/
 411         Bool_t findable = IsFindable(mcEvent,iMc);
 412
 413     Bool_t recStatus = kFALSE;
 414         Int_t nClones = 0, nFakes = 0;
 415     for(Int_t iRec=0; iRec<multRec; ++iRec)
 416     {
 417       // check reconstructed
 418       if(iMc == labelsRecSec[iRec])
 419       {
 420         if (recStatus && nClones < fgkMaxClones) nClones++;
 421                 recStatus = kTRUE;
 422       }
 423           //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
 424           if (labelsRecSec[iRec] < 0 && -labelsRecSec[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
 425     }
 426
 427     // Only 5 charged particle species (e,mu,pi,K,p)
 428     if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
 429
 430     // transform Pdg to Pid
 431     Int_t pid = TransformToPID(particle);
 432
 433     Float_t mceta =  particle->Eta();
 434     Float_t mcphi =  particle->Phi();
 435     if(mcphi<0) mcphi += 2.*TMath::Pi();
 436     Float_t mcpt = particle->Pt();
 437     Float_t mcR = particle->R();
 438
 439     // get info about mother
 440     Int_t motherLabel = particle->GetMother(0);
 441     if(motherLabel < 0) continue;
 442     TParticle *mother = stack->Particle(motherLabel);
 443     if(!mother) continue;
 444
 445     Float_t mother_eta = mother->Eta();
 446     Float_t mother_phi = mother->Phi();
 447     if(mother_phi<0) mother_phi += 2.*TMath::Pi();
 448
 449     Float_t charge = 0.;
 450     if (particle->GetPDG()->Charge() < 0)  charge = -1.;
 451     else if (particle->GetPDG()->Charge() > 0)  charge = 1.;
 452
 453     // Fill histograms
 454     Double_t vEffSecHisto[12] = { mceta, mcphi, mcpt, pid, recStatus, findable, mcR, mother_phi, mother_eta, charge, nClones, nFakes };
 455     fEffSecHisto->Fill(vEffSecHisto);
 456   }
 457   }
 458   }
 459
 460   if(labelsRecSec) delete [] labelsRecSec; labelsRecSec = 0;
 461   if(labelsAllRecSec) delete [] labelsAllRecSec; labelsAllRecSec = 0;
 462 }
 463
 464
 465
 466
 467 //_____________________________________________________________________________
 468 void AliPerformanceEff::ProcessTPCITS(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
 469 {
 470   // Fill efficiency comparison information
 471
 472   if(!esdEvent) return;
 473   if(!mcEvent) return;
 474
 475   AliStack *stack = mcEvent->Stack();
 476   if (!stack) {
 477     AliDebug(AliLog::kError, "Stack not available");
 478     return;
 479   }
 480
 481   Int_t *labelsRecTPCITS =  NULL;
 482   labelsRecTPCITS =  new Int_t[esdEvent->GetNumberOfTracks()];
 483   if(!labelsRecTPCITS)
 484   {
 485      Printf("Cannot create labelsRecTPCITS");
 486      return;
 487   }
 488   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecTPCITS[i] = 0; }
 489
 490   Int_t *labelsAllRecTPCITS =  NULL;
 491   labelsAllRecTPCITS =  new Int_t[esdEvent->GetNumberOfTracks()];
 492   if(!labelsAllRecTPCITS) {
 493      delete [] labelsRecTPCITS;
 494      Printf("Cannot create labelsAllRecTPCITS");
 495      return;
 496   }
 497   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecTPCITS[i] = 0; }
 498
 499   // loop over rec. tracks
 500   AliESDtrack *track=0;
 501   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 502   {
 503     track = esdEvent->GetTrack(iTrack);
 504     if(!track) continue;
 505     if(track->Charge()==0) continue;
 506
 507     // if not fUseKinkDaughters don't use tracks with kink index > 0
 508     if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
 509
 510     //Int_t label = TMath::Abs(track->GetLabel());
 511         Int_t label = track->GetLabel();  //Use global label for combined efficiency analysis
 512     labelsAllRecTPCITS[iTrack]=label;
 513
 514     // iTPC+ITS
 515     if(IsRecTPCITS(track) != 0)
 516       labelsRecTPCITS[iTrack]=label;
 517   }
 518
 519   //
 520   // MC histograms for efficiency studies
 521   //
 522   //Int_t nPart  = stack->GetNtrack();
 523   Int_t nPart  = stack->GetNprimary();
 524   for (Int_t iMc = 0; iMc < nPart; ++iMc)
 525   {
 526         if (iMc == 0) continue;         //Cannot distinguish between track or fake track
 527     TParticle* particle = stack->Particle(iMc);
 528     if (!particle) continue;
 529     if (!particle->GetPDG()) continue;
 530     if (particle->GetPDG()->Charge() == 0.0) continue;
 531
 532     // physical primary
 533     Bool_t prim = stack->IsPhysicalPrimary(iMc);
 534     if(!prim) continue;
 535
 536     /*Bool_t findable = kFALSE;
 537     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 538     {
 539       // check findable
 540       if(iMc == labelsAllRecTPCITS[iRec])
 541       {
 542         findable = IsFindable(mcEvent,iMc);
 543         break;
 544       }
 545     }*/
 546         Bool_t findable = IsFindable(mcEvent,iMc);
 547
 548     Bool_t recStatus = kFALSE;
 549         Int_t nClones = 0, nFakes = 0;
 550     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 551     {
 552       // check reconstructed
 553       if(iMc == labelsRecTPCITS[iRec])
 554       {
 555         if (recStatus && nClones < fgkMaxClones) nClones++;
 556                 recStatus = kTRUE;
 557       }
 558           //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
 559           if (labelsRecTPCITS[iRec] < 0 && -labelsRecTPCITS[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
 560     }
 561
 562     // Only 5 charged particle species (e,mu,pi,K,p)
 563     if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
 564
 565     // transform Pdg to Pid
 566     Int_t pid = TransformToPID(particle);
 567
 568     Float_t mceta =  particle->Eta();
 569     Float_t mcphi =  particle->Phi();
 570     if(mcphi<0) mcphi += 2.*TMath::Pi();
 571     Float_t mcpt = particle->Pt();
 572
 573     Float_t charge = 0.;
 574     if (particle->GetPDG()->Charge() < 0)  charge = -1.;
 575     else if (particle->GetPDG()->Charge() > 0)  charge = 1.;
 576
 577     // Fill histograms
 578     Double_t vEffHisto[9] = { mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes};
 579     fEffHisto->Fill(vEffHisto);
 580   }
 581
 582   if(labelsRecTPCITS) delete [] labelsRecTPCITS; labelsRecTPCITS = 0;
 583   if(labelsAllRecTPCITS) delete [] labelsAllRecTPCITS; labelsAllRecTPCITS = 0;
 584 }
 585
 586 //_____________________________________________________________________________
 587 void AliPerformanceEff::ProcessConstrained(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent)
 588 {
 589   // Process comparison information
 590   if(!esdEvent) return;
 591   if(!mcEvent) return;
 592
 593   AliStack *stack = mcEvent->Stack();
 594   if (!stack) {
 595     AliDebug(AliLog::kError, "Stack not available");
 596     return;
 597   }
 598
 599   Int_t *labelsRecConstrained =  NULL;
 600   labelsRecConstrained =  new Int_t[esdEvent->GetNumberOfTracks()];
 601   if(!labelsRecConstrained)
 602   {
 603      Printf("Cannot create labelsRecConstrained");
 604      return;
 605   }
 606   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsRecConstrained[i] = 0; }
 607
 608   Int_t *labelsAllRecConstrained =  NULL;
 609   labelsAllRecConstrained =  new Int_t[esdEvent->GetNumberOfTracks()];
 610   if(!labelsAllRecConstrained) {
 611      delete [] labelsRecConstrained;
 612      Printf("Cannot create labelsAllRecConstrained");
 613      return;
 614   }
 615   for(Int_t i=0;i<esdEvent->GetNumberOfTracks();i++) { labelsAllRecConstrained[i] = 0; }
 616
 617   // loop over rec. tracks
 618   AliESDtrack *track=0;
 619   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 620   {
 621     track = esdEvent->GetTrack(iTrack);
 622     if(!track) continue;
 623     if(track->Charge()==0) continue;
 624
 625     // if not fUseKinkDaughters don't use tracks with kink index > 0
 626     if(!fUseKinkDaughters && track->GetKinkIndex(0) > 0) continue;
 627
 628     //Int_t label = TMath::Abs(track->GetLabel());
 629         Int_t label = track->GetLabel();
 630     labelsAllRecConstrained[iTrack]=label;
 631
 632     // Constrained
 633     if(IsRecConstrained(track) != 0)
 634       labelsRecConstrained[iTrack]=label;
 635
 636   }
 637
 638   //
 639   // MC histograms for efficiency studies
 640   //
 641
 642
 643   //Int_t nPart  = stack->GetNtrack();
 644   Int_t nPart  = stack->GetNprimary();
 645   for (Int_t iMc = 0; iMc < nPart; ++iMc)
 646   {
 647         if (iMc == 0) continue;         //Cannot distinguish between track or fake track
 648     TParticle* particle = stack->Particle(iMc);
 649     if (!particle) continue;
 650     if (!particle->GetPDG()) continue;
 651     if (particle->GetPDG()->Charge() == 0.0) continue;
 652
 653     // physical primary
 654     Bool_t prim = stack->IsPhysicalPrimary(iMc);
 655     if(!prim) continue;
 656
 657     /*Bool_t findable = kFALSE;
 658     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 659     {
 660       // check findable
 661       if(iMc == labelsAllRecConstrained[iRec])
 662       {
 663         findable = IsFindable(mcEvent,iMc);
 664         break;
 665       }
 666     }*/
 667         Bool_t findable = IsFindable(mcEvent,iMc);
 668
 669     Bool_t recStatus = kFALSE;
 670         Int_t nClones = 0, nFakes = 0;
 671     for(Int_t iRec=0; iRec<esdEvent->GetNumberOfTracks(); ++iRec)
 672     {
 673       // check reconstructed
 674       if(iMc == labelsRecConstrained[iRec])
 675       {
 676         if (recStatus && nClones < fgkMaxClones) nClones++;
 677                 recStatus = kTRUE;
 678       }
 679           //In order to relate the fake track to track parameters, we assign it to the best matching ESD track.
 680           if (labelsRecConstrained[iRec] < 0 && -labelsRecConstrained[iRec] == iMc && nFakes < fgkMaxFakes) nFakes++;
 681     }
 682
 683     // Only 5 charged particle species (e,mu,pi,K,p)
 684     if (fCutsMC->IsPdgParticle(TMath::Abs(particle->GetPdgCode())) == kFALSE) continue;
 685
 686     // transform Pdg to Pid
 687     Int_t pid = TransformToPID(particle);
 688
 689     Float_t mceta =  particle->Eta();
 690     Float_t mcphi =  particle->Phi();
 691     if(mcphi<0) mcphi += 2.*TMath::Pi();
 692     Float_t mcpt = particle->Pt();
 693
 694     Float_t charge = 0.;
 695     if (particle->GetPDG()->Charge() < 0)  charge = -1.;
 696     else if (particle->GetPDG()->Charge() > 0)  charge = 1.;
 697
 698     // Fill histograms
 699     Double_t vEffHisto[9] = { mceta, mcphi, mcpt, pid, recStatus, findable, charge, nClones, nFakes };
 700     fEffHisto->Fill(vEffHisto);
 701   }
 702
 703   if(labelsRecConstrained) delete [] labelsRecConstrained; labelsRecConstrained = 0;
 704   if(labelsAllRecConstrained) delete [] labelsAllRecConstrained; labelsAllRecConstrained = 0;
 705 }
 706
 707 //_____________________________________________________________________________
 708 void AliPerformanceEff::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 709 {
 710   // Process comparison information
 711   //
 712   if(!esdEvent)
 713   {
 714     Error("Exec","esdEvent not available");
 715     return;
 716   }
 717   AliHeader* header = 0;
 718   AliGenEventHeader* genHeader = 0;
 719   AliStack* stack = 0;
 720   TArrayF vtxMC(3);
 721
 722   if(bUseMC)
 723   {
 724     if(!mcEvent) {
 725       Error("Exec","mcEvent not available");
 726       return;
 727     }
 728     // get MC event header
 729     header = mcEvent->Header();
 730     if (!header) {
 731       Error("Exec","Header not available");
 732       return;
 733     }
 734     // MC particle stack
 735     stack = mcEvent->Stack();
 736     if (!stack) {
 737       Error("Exec","Stack not available");
 738       return;
 739     }
 740     // get MC vertex
 741     genHeader = header->GenEventHeader();
 742     if (!genHeader) {
 743       Error("Exec","Could not retrieve genHeader from Header");
 744       return;
 745     }
 746     genHeader->PrimaryVertex(vtxMC);
 747   }
 748   else {
 749     Error("Exec","MC information required!");
 750     return;
 751   }
 752
 753   // use ESD friends
 754   if(bUseESDfriend) {
 755     if(!esdFriend) {
 756       Error("Exec","esdFriend not available");
 757       return;
 758     }
 759   }
 760
 761   //
 762   //  Process events
 763   //
 764   if(GetAnalysisMode() == 0) ProcessTPC(mcEvent,esdEvent);
 765   else if(GetAnalysisMode() == 1) ProcessTPCITS(mcEvent,esdEvent);
 766   else if(GetAnalysisMode() == 2) ProcessConstrained(mcEvent,esdEvent);
 767   else if(GetAnalysisMode() == 5) ProcessTPCSec(mcEvent,esdEvent);
 768   else {
 769     printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 770     return;
 771   }
 772 }
 773
 774 //_____________________________________________________________________________
 775 Int_t AliPerformanceEff::TransformToPID(TParticle *particle)
 776 {
 777 // transform Pdg to Pid
 778 // Pdg convension is different for hadrons and leptons
 779 // (e.g. K+/K- = 321/-321; e+/e- = -11/11 )
 780
 781   Int_t pid = -1;
 782   if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetEM() ) pid = 0;
 783   if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetMuM() ) pid = 1;
 784   if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetPiP() ) pid = 2;
 785   if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetKP() ) pid = 3;
 786   if( TMath::Abs(particle->GetPdgCode())==fCutsMC->GetProt() ) pid = 4;
 787
 788 return pid;
 789 }
 790
 791 //_____________________________________________________________________________
 792 Bool_t AliPerformanceEff::IsFindable(const AliMCEvent *mcEvent, Int_t label)
 793 {
 794 //
 795 // Findfindable tracks
 796 //
 797 if(!mcEvent) return kFALSE;
 798
 799   AliMCParticle *mcParticle = (AliMCParticle*) mcEvent->GetTrack(label);
 800   if(!mcParticle) return kFALSE;
 801
 802   Int_t counter;
 803   Float_t tpcTrackLength = mcParticle->GetTPCTrackLength(AliTracker::GetBz(),0.05,counter,3.0);
 804   //printf("tpcTrackLength %f \n", tpcTrackLength);
 805
 806 return (tpcTrackLength>fCutsMC->GetMinTrackLength());
 807 }
 808
 809 //_____________________________________________________________________________
 810 Bool_t AliPerformanceEff::IsRecTPC(AliESDtrack *esdTrack)
 811 {
 812 //
 813 // Check whether track is reconstructed in TPC
 814 //
 815 if(!esdTrack) return kFALSE;
 816
 817   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 818   if(!track) return kFALSE;
 819
 820   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 821   esdTrack->GetImpactParametersTPC(dca,cov);
 822
 823   Bool_t recStatus = kFALSE;
 824   if(esdTrack->GetTPCNcls()>fCutsRC->GetMinNClustersTPC()) recStatus = kTRUE;
 825
 826   /*
 827   if( TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
 828       TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 829   {
 830     recStatus = kTRUE;
 831   }
 832   */
 833
 834 return recStatus;
 835 }
 836
 837 //_____________________________________________________________________________
 838 Bool_t AliPerformanceEff::IsRecTPCITS(AliESDtrack *esdTrack)
 839 {
 840 //
 841 // Check whether track is reconstructed in TPCITS
 842 //
 843 if(!esdTrack) return kFALSE;
 844
 845   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 846   esdTrack->GetImpactParameters(dca,cov);
 847
 848   Bool_t recStatus = kFALSE;
 849
 850   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return kFALSE; // TPC refit
 851   if ((esdTrack->GetStatus()&AliESDtrack::kITSrefit)==0) return kFALSE; // ITS refit
 852   if (esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return kFALSE;  // min. nb. ITS clusters
 853   //if ((esdTrack->GetStatus()&AliESDtrack::kITSrefit)==0) return kFALSE; // ITS refit
 854   //Int_t clusterITS[200];
 855   //if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return kFALSE;  // min. nb. ITS clusters
 856
 857   recStatus = kTRUE;
 858   /*
 859   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
 860      TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 861   {
 862     recStatus = kTRUE;
 863   }
 864   */
 865
 866 return recStatus;
 867 }
 868
 869 //_____________________________________________________________________________
 870 Bool_t AliPerformanceEff::IsRecConstrained(AliESDtrack *esdTrack)
 871 {
 872 //
 873 // Check whether track is reconstructed in IsRecConstrained
 874 //
 875   if(!esdTrack) return kFALSE;
 876
 877   const AliExternalTrackParam * track = esdTrack->GetConstrainedParam();
 878   if(!track) return kFALSE;
 879
 880   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 881   esdTrack->GetImpactParameters(dca,cov);
 882   //Int_t label = TMath::Abs(esdTrack->GetLabel());
 883
 884   Bool_t recStatus = kFALSE;
 885
 886   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return kFALSE; // TPC refit
 887   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return kFALSE; // min. nb. TPC clusters
 888   Int_t clusterITS[200];
 889   if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return kFALSE;  // min. nb. ITS clusters
 890
 891   recStatus = kTRUE;
 892
 893   /*
 894   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() &&
 895      TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ())
 896   {
 897     recStatus = kTRUE;
 898   }
 899   */
 900
 901 return recStatus;
 902 }
 903
 904 //_____________________________________________________________________________
 905 Long64_t AliPerformanceEff::Merge(TCollection* const list)
 906 {
 907   // Merge list of objects (needed by PROOF)
 908
 909   if (!list)
 910   return 0;
 911
 912   if (list->IsEmpty())
 913   return 1;
 914
 915   TIterator* iter = list->MakeIterator();
 916   TObject* obj = 0;
 917
 918   // collection of generated histograms
 919
 920   Int_t count=0;
 921   while((obj = iter->Next()) != 0)
 922   {
 923     AliPerformanceEff* entry = dynamic_cast<AliPerformanceEff*>(obj);
 924     if (entry == 0) continue;
 925
 926      fEffHisto->Add(entry->fEffHisto);
 927      fEffSecHisto->Add(entry->fEffSecHisto);
 928   count++;
 929   }
 930
 931 return count;
 932 }
 933
 934 //_____________________________________________________________________________
 935 void AliPerformanceEff::Analyse()
 936 {
 937   // Analyse comparison information and store output histograms
 938   // in the folder "folderEff"
 939   //
 940   TH1::AddDirectory(kFALSE);
 941   TObjArray *aFolderObj = new TObjArray;
 942   if(!aFolderObj) return;
 943   //  char title[256];
 944
 945   //
 946   // efficiency vs pt
 947   //
 948
 949   if(GetAnalysisMode() != 5) {
 950
 951   fEffHisto->GetAxis(0)->SetRangeUser(-0.9,0.89); // eta range
 952   fEffHisto->GetAxis(2)->SetRangeUser(0.1,10.);   // pt range
 953
 954   // rec efficiency vs pt
 955   fEffHisto->GetAxis(3)->SetRangeUser(0.,3.99);  // reconstructed
 956
 957   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
 958   aFolderObj->Add(AddHistoEff(2, "ptRecEff", "rec. efficiency", 0));
 959   aFolderObj->Add(AddHistoEff(2, "ptClone", "clone rate", 1));
 960   aFolderObj->Add(AddHistoEff(2, "ptFake", "fake rate", 2));
 961
 962   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
 963   aFolderObj->Add(AddHistoEff(2, "ptRecEffNeg", "rec. efficiency neg.", 0));
 964
 965   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
 966   aFolderObj->Add(AddHistoEff(2, "ptRecEffPos", "rec. efficiency pos.", 0));
 967
 968   // rec efficiency vs pid vs pt
 969   fEffHisto->GetAxis(3)->SetRangeUser(2.,2.);    // pions
 970
 971   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
 972   aFolderObj->Add(AddHistoEff(2, "ptRecEffPi", "rec. efficiency (pions)", 0));
 973
 974   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
 975   aFolderObj->Add(AddHistoEff(2, "ptRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
 976
 977   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
 978   aFolderObj->Add(AddHistoEff(2, "ptRecEffPiPos", "rec. efficiency (pions) pos.", 0));
 979
 980
 981   fEffHisto->GetAxis(3)->SetRangeUser(3.,3.);    // kaons
 982
 983   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
 984   aFolderObj->Add(AddHistoEff(2, "ptRecEffK", "rec. efficiency (kaons)", 0));
 985
 986   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
 987   aFolderObj->Add(AddHistoEff(2, "ptRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
 988
 989   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
 990   aFolderObj->Add(AddHistoEff(2, "ptRecEffKPos", "rec. efficiency (kaons) pos.", 0));
 991
 992
 993   fEffHisto->GetAxis(3)->SetRangeUser(4.,4.);    // protons
 994
 995   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
 996   aFolderObj->Add(AddHistoEff(2, "ptRecEffP", "rec. efficiency (protons)", 0));
 997
 998   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
 999   aFolderObj->Add(AddHistoEff(2, "ptRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1000
1001   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1002   aFolderObj->Add(AddHistoEff(2, "ptRecEffPPos", "rec. efficiency (protons) pos.", 0));
1003
1004   // findable efficiency vs pt
1005   fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1006   fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1007
1008   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1009   aFolderObj->Add(AddHistoEff(2, "ptRecEffF", "rec. efficiency (findable)", 0));
1010   aFolderObj->Add(AddHistoEff(2, "ptCloneF", "clone rate (findable)", 1));
1011   aFolderObj->Add(AddHistoEff(2, "ptFakeF", "fake rate (findable)", 2));
1012
1013
1014   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1015   aFolderObj->Add(AddHistoEff(2, "ptRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1016
1017   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1018   aFolderObj->Add(AddHistoEff(2, "ptRecEffFPos", "rec. efficiency (findable) pos.", 0));
1019
1020   //
1021   // efficiency vs eta
1022   //
1023
1024   fEffHisto->GetAxis(0)->SetRangeUser(-1.5,1.5); // eta range
1025   fEffHisto->GetAxis(2)->SetRangeUser(0.2,10.); // pt range
1026   fEffHisto->GetAxis(5)->SetRangeUser(0.,1.);   // all
1027
1028   // rec efficiency vs eta
1029   fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);  // reconstructed
1030
1031   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1032   aFolderObj->Add(AddHistoEff(0, "etaRecEff", "rec. efficiency", 0));
1033   aFolderObj->Add(AddHistoEff(0, "etaClone", "clone rate", 1));
1034   aFolderObj->Add(AddHistoEff(0, "etaFake", "fake rate", 2));
1035
1036
1037   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1038   aFolderObj->Add(AddHistoEff(0, "etaRecEffNeg", "rec. efficiency neg.", 0));
1039
1040   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1041   aFolderObj->Add(AddHistoEff(0, "etaRecEffPos", "rec. efficiency pos.", 0));
1042
1043   // rec efficiency vs pid vs eta
1044   fEffHisto->GetAxis(3)->SetRangeUser(2.,2.);    // pions
1045
1046   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1047   aFolderObj->Add(AddHistoEff(0, "etaRecEffPi", "rec. efficiency (pions)", 0));
1048
1049   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1050   aFolderObj->Add(AddHistoEff(0, "etaRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
1051
1052   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1053   aFolderObj->Add(AddHistoEff(0, "etaRecEffPiPos", "rec. efficiency (pions) pos.", 0));
1054
1055
1056   fEffHisto->GetAxis(3)->SetRangeUser(3.,3.);    // kaons
1057
1058   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1059   aFolderObj->Add(AddHistoEff(0, "etaRecEffK", "rec. efficiency (kaons)", 0));
1060
1061   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1062   aFolderObj->Add(AddHistoEff(0, "etaRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
1063
1064   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1065   aFolderObj->Add(AddHistoEff(0, "etaRecEffKPos", "rec. efficiency (kaons) pos.", 0));
1066
1067
1068   fEffHisto->GetAxis(3)->SetRangeUser(4.,4.);    // protons
1069
1070   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1071   aFolderObj->Add(AddHistoEff(0, "etaRecEffP", "rec. efficiency (protons)", 0));
1072
1073   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1074   aFolderObj->Add(AddHistoEff(0, "etaRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1075
1076   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1077   aFolderObj->Add(AddHistoEff(0, "etaRecEffPPos", "rec. efficiency (protons) pos.", 0));
1078
1079   // findable efficiency vs eta
1080   fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1081   fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1082
1083   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1084   aFolderObj->Add(AddHistoEff(0, "etaRecEffF", "rec. efficiency (findable)", 0));
1085   aFolderObj->Add(AddHistoEff(0, "etaCloneF", "clone rate (findable)", 1));
1086   aFolderObj->Add(AddHistoEff(0, "etaFakeF", "fake rate (findable)", 2));
1087
1088
1089   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1090   aFolderObj->Add(AddHistoEff(0, "etaRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1091
1092   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1093   aFolderObj->Add(AddHistoEff(0, "etaRecEffFPos", "rec. efficiency (findable) pos.", 0));
1094
1095   //
1096   // efficiency vs phi
1097   //
1098
1099   fEffHisto->GetAxis(0)->SetRangeUser(-0.9,0.9); // eta range
1100   fEffHisto->GetAxis(2)->SetRangeUser(0.2,10.); // pt range
1101   fEffHisto->GetAxis(5)->SetRangeUser(0.,1.);   // all
1102
1103   // rec efficiency vs phi
1104   fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);  // reconstructed
1105
1106   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1107   aFolderObj->Add(AddHistoEff(1, "phiRecEff", "rec. efficiency", 0));
1108   aFolderObj->Add(AddHistoEff(1, "phiClone", "clone rate", 1));
1109   aFolderObj->Add(AddHistoEff(1, "phiFake", "fake rate", 2));
1110
1111   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1112   aFolderObj->Add(AddHistoEff(1, "phiRecEffNeg", "rec. efficiency neg.", 0));
1113
1114   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1115   aFolderObj->Add(AddHistoEff(1, "phiRecEffPos", "rec. efficiency pos.", 0));
1116
1117   // rec efficiency vs pid vs phi
1118   fEffHisto->GetAxis(3)->SetRangeUser(2.,2.);    // pions
1119
1120   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1121   aFolderObj->Add(AddHistoEff(1, "phiRecEffPi", "rec. efficiency (pions)", 0));
1122
1123   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1124   aFolderObj->Add(AddHistoEff(1, "phiRecEffPiNeg", "rec. efficiency (pions) neg.", 0));
1125
1126   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1127   aFolderObj->Add(AddHistoEff(1, "phiRecEffPiPos", "rec. efficiency (pions) pos.", 0));
1128
1129
1130   fEffHisto->GetAxis(3)->SetRangeUser(3.,3.);    // kaons
1131
1132   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1133   aFolderObj->Add(AddHistoEff(1, "phiRecEffK", "rec. efficiency (kaons)", 0));
1134
1135   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1136   aFolderObj->Add(AddHistoEff(1, "phiRecEffKNeg", "rec. efficiency (kaons) neg.", 0));
1137
1138   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1139   aFolderObj->Add(AddHistoEff(1, "phiRecEffKPos", "rec. efficiency (kaons) pos.", 0));
1140
1141
1142   fEffHisto->GetAxis(3)->SetRangeUser(4.,4.);    // protons
1143
1144   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1145   aFolderObj->Add(AddHistoEff(1, "phiRecEffP", "rec. efficiency (protons)", 0));
1146
1147   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1148   aFolderObj->Add(AddHistoEff(1, "phiRecEffPNeg", "rec. efficiency (protons) neg.", 0));
1149
1150   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1151   aFolderObj->Add(AddHistoEff(1, "phiRecEffPPos", "rec. efficiency (protons) pos.", 0));
1152
1153   // findable efficiency vs phi
1154   fEffHisto->GetAxis(3)->SetRangeUser(0.,4.);
1155   fEffHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1156
1157   fEffHisto->GetAxis(6)->SetRangeUser(-2.,2.);   // charge all
1158   aFolderObj->Add(AddHistoEff(1, "phiRecEffF", "rec. efficiency (findable)", 0));
1159   aFolderObj->Add(AddHistoEff(1, "phiCloneF", "clone rate (findable)", 1));
1160   aFolderObj->Add(AddHistoEff(1, "phiFakeF", "fake rate (findable)", 2));
1161
1162
1163   fEffHisto->GetAxis(6)->SetRangeUser(-2.,0.);   // charge negativ
1164   aFolderObj->Add(AddHistoEff(1, "phiRecEffFNeg", "rec. efficiency (findable) neg.", 0));
1165
1166   fEffHisto->GetAxis(6)->SetRangeUser(0.,2.);    // charge positiv
1167   aFolderObj->Add(AddHistoEff(1, "phiRecEffFPos", "rec. efficiency (findable) pos.", 0));
1168   }
1169   else {
1170   //
1171   Float_t minEta=-1.5, maxEta=1.5;
1172   Float_t minR=0.0, maxR=150.0;
1173   Float_t minPt=0.15, maxPt=100.0;
1174
1175   // mother eta range
1176   fEffSecHisto->GetAxis(8)->SetRangeUser(minEta,maxEta);
1177
1178   // particle creation radius range
1179   fEffSecHisto->GetAxis(6)->SetRangeUser(minR,maxR);
1180
1181   //
1182   fEffSecHisto->GetAxis(0)->SetRangeUser(minEta,maxEta);
1183   fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1184
1185   // rec efficiency vs pt
1186
1187   aFolderObj->Add(AddHistoEff(2, "ptRecEff", "rec. efficiency", 0, 1));
1188   aFolderObj->Add(AddHistoEff(2, "ptClone", "clone rate", 1, 1));
1189   aFolderObj->Add(AddHistoEff(2, "ptFake", "fake rate", 2, 1));
1190
1191   // rec efficiency vs pid vs pt
1192
1193   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.); // electrons
1194   aFolderObj->Add(AddHistoEff(2, "ptRecEffEle", "rec. efficiency (electrons)", 0, 1));
1195
1196   fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1197   aFolderObj->Add(AddHistoEff(2, "ptRecEffPi", "rec. efficiency (pions)", 0, 1));
1198
1199   fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1200   aFolderObj->Add(AddHistoEff(2, "ptRecEffK", "rec. efficiency (kaons)", 0, 1));
1201
1202   fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1203   aFolderObj->Add(AddHistoEff(2, "ptRecEffP", "rec. efficiency (protons)", 0, 1));
1204
1205   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1206
1207   // findable efficiency vs pt
1208
1209   fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1210   aFolderObj->Add(AddHistoEff(2, "ptRecEffF", "rec. efficiency (findable)", 0, 1));
1211   aFolderObj->Add(AddHistoEff(2, "ptCloneF", "clone rate (findable)", 1, 1));
1212   aFolderObj->Add(AddHistoEff(2, "ptFakeF", "fake rate (findable)", 2, 1));
1213
1214   //
1215   // efficiency vs eta
1216   //
1217   fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1218   fEffSecHisto->GetAxis(4)->SetRangeUser(0.,1.);   // all
1219   fEffSecHisto->GetAxis(5)->SetRangeUser(0.,1.);   // all
1220
1221   aFolderObj->Add(AddHistoEff(0, "etaRecEff", "rec. efficiency", 0, 1));
1222   aFolderObj->Add(AddHistoEff(0, "etaClone", "clone rate", 1, 1));
1223   aFolderObj->Add(AddHistoEff(0, "etaFake", "fake rate", 2, 1));
1224
1225   // rec efficiency vs pid vs eta
1226   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.); // electrons
1227   aFolderObj->Add(AddHistoEff(0, "etaRecEffEle", "rec. efficiency (electrons)", 0, 1));
1228
1229   fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1230   aFolderObj->Add(AddHistoEff(0, "etaRecEffPi", "rec. efficiency (pions)", 0, 1));
1231
1232   fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1233   aFolderObj->Add(AddHistoEff(0, "etaRecEffK", "rec. efficiency (kaons)", 0, 1));
1234
1235   fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1236   aFolderObj->Add(AddHistoEff(0, "etaRecEffP", "rec. efficiency (protons)", 0, 1));
1237
1238   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1239
1240   // findable efficiency vs eta
1241
1242   fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1243   aFolderObj->Add(AddHistoEff(0, "etaRecEffF", "rec. efficiency (findable)", 0, 1));
1244   aFolderObj->Add(AddHistoEff(0, "etaCloneF", "clone rate (findable)", 1, 1));
1245   aFolderObj->Add(AddHistoEff(0, "etaFakeF", "fake rate (findable)", 2, 1));
1246
1247   //
1248   // efficiency vs phi
1249   //
1250
1251   fEffSecHisto->GetAxis(0)->SetRangeUser(minEta,maxEta);
1252   fEffSecHisto->GetAxis(2)->SetRangeUser(minPt,maxPt);
1253
1254   fEffSecHisto->GetAxis(4)->SetRangeUser(0.,1.);   // all
1255   fEffSecHisto->GetAxis(5)->SetRangeUser(0.,1.);   // all
1256
1257   aFolderObj->Add(AddHistoEff(1, "phiRecEff", "rec. efficiency", 0, 1));
1258   aFolderObj->Add(AddHistoEff(1, "phiClone", "clone rate", 1, 1));
1259   aFolderObj->Add(AddHistoEff(1, "phiFake", "fake rate", 2, 1));
1260
1261   // rec efficiency vs pid vs phi
1262   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,0.);
1263   aFolderObj->Add(AddHistoEff(1, "phiRecEffEle", "rec. efficiency (electrons)", 0, 1));
1264
1265   fEffSecHisto->GetAxis(3)->SetRangeUser(2.,2.); // pions
1266   aFolderObj->Add(AddHistoEff(1, "phiRecEffPi", "rec. efficiency (pions)", 0, 1));
1267
1268   fEffSecHisto->GetAxis(3)->SetRangeUser(3.,3.); // kaons
1269   aFolderObj->Add(AddHistoEff(1, "phiRecEffK", "rec. efficiency (kaons)", 0, 1));
1270
1271   fEffSecHisto->GetAxis(3)->SetRangeUser(4.,4.); // protons
1272   aFolderObj->Add(AddHistoEff(1, "phiRecEffP", "rec. efficiency (protons)", 0, 1));
1273
1274   fEffSecHisto->GetAxis(3)->SetRangeUser(0.,4.);
1275
1276   // findable efficiency vs phi
1277
1278   fEffSecHisto->GetAxis(5)->SetRangeUser(1.,1.); // findable
1279   aFolderObj->Add(AddHistoEff(1, "phiRecEffF", "rec. efficiency (findable)", 0, 1));
1280   aFolderObj->Add(AddHistoEff(1, "phiCloneF", "clone rate (findable)", 1, 1));
1281   aFolderObj->Add(AddHistoEff(1, "phiFakeF", "fake rate (findable)", 2, 1));
1282   }
1283
1284   for (Int_t i = 0;i < fEffHisto->GetNdimensions();i++)
1285   {
1286           fEffHisto->GetAxis(i)->SetRange(1,0);                         //Reset Range
1287   }
1288
1289   // export objects to analysis folder
1290   fAnalysisFolder = ExportToFolder(aFolderObj);
1291
1292   // delete only TObjArray
1293   if(aFolderObj) delete aFolderObj;
1294 }
1295
1296 //_____________________________________________________________________________
1297 TFolder* AliPerformanceEff::ExportToFolder(TObjArray * array)
1298 {
1299   // recreate folder avery time and export objects to new one
1300   //
1301   AliPerformanceEff * comp=this;
1302   TFolder *folder = comp->GetAnalysisFolder();
1303
1304   TString name, title;
1305   TFolder *newFolder = 0;
1306   Int_t i = 0;
1307   Int_t size = array->GetSize();
1308
1309   if(folder) {
1310      // get name and title from old folder
1311      name = folder->GetName();
1312      title = folder->GetTitle();
1313
1314          // delete old one
1315      delete folder;
1316
1317          // create new one
1318      newFolder = CreateFolder(name.Data(),title.Data());
1319      newFolder->SetOwner();
1320
1321          // add objects to folder
1322      while(i < size) {
1323            newFolder->Add(array->At(i));
1324            i++;
1325          }
1326   }
1327
1328 return newFolder;
1329 }
1330
1331
1332 //_____________________________________________________________________________
1333 TFolder* AliPerformanceEff::CreateFolder(TString name,TString title) {
1334 // create folder for analysed histograms
1335 //
1336 TFolder *folder = 0;
1337   folder = new TFolder(name.Data(),title.Data());
1338
1339   return folder;
1340 }
1341
1342 TH1D* WeightedProjection(THnSparseF* src, Int_t axis, Int_t nWeights, Int_t* weightCoords)
1343 {
1344     THnSparseF* tmp = (THnSparseF*) src->Clone();
1345     Int_t i;
1346     for (i = 0;i < tmp->GetNbins();i++)
1347     {
1348         Int_t coords[12];
1349         tmp->GetBinContent(i, coords);
1350         Int_t weight = 0, j;
1351         for (j = 0;j < nWeights;j++)
1352         {
1353                         //The coordinates range from 1 to maxClones / maxFakes + 1, so we have to subtract one
1354             weight += coords[weightCoords[j]] - 1;
1355         }
1356         tmp->SetBinContent(i, weight);
1357     }
1358
1359     TH1D* ret = tmp->Projection(axis);
1360     delete tmp;
1361     return(ret);
1362 }
1363
1364 //_____________________________________________________________________________
1365 TH1D* AliPerformanceEff::AddHistoEff(Int_t axis, const Char_t *name, const Char_t* vsTitle,
1366         const Int_t type, const Int_t secondary) {
1367   // Create and add rec efficiency vs pt, eta, phi
1368
1369   char title[256];
1370
1371   TH1D *recc = NULL;
1372
1373   THnSparseF* EffHisto = secondary ? fEffSecHisto : fEffHisto;
1374
1375   Int_t axis_clone = secondary ? 10 : 7;
1376   Int_t axis_fake = secondary ? 11 : 8;
1377   Int_t axis_all[3] = {4, axis_clone, axis_fake};
1378
1379
1380   if (type == 0) // Efficiency
1381   {
1382       EffHisto->GetAxis(4)->SetRangeUser(0.,1.);    // all
1383       TH1D *all = EffHisto->Projection(axis);
1384
1385       EffHisto->GetAxis(4)->SetRangeUser(1.,1.);    // reconstructed
1386       TH1D *rec = EffHisto->Projection(axis);
1387       recc = (TH1D*)rec->Clone();
1388
1389       if(recc)
1390       {
1391         recc->Divide(rec,all,1,1,"B");
1392         recc->GetYaxis()->SetTitle("efficiency");
1393       }
1394   }
1395   else if (type == 1) // Clone Rate
1396   {
1397       EffHisto->GetAxis(4)->SetRangeUser(1.,1.);    // reconstructed
1398       TH1D *all = WeightedProjection(EffHisto, axis, 3, axis_all);
1399
1400       EffHisto->GetAxis(4)->SetRangeUser(1.,1.);    // reconstructed
1401       TH1D *clone = WeightedProjection(EffHisto, axis, 1, &axis_clone);
1402       recc = (TH1D*) clone->Clone();
1403
1404       if(recc)
1405       {
1406         recc->Divide(clone,all,1,1,"B");
1407         recc->GetYaxis()->SetTitle("clone rate");
1408       }
1409   }
1410   else if (type == 2) // Fake Rate
1411   {
1412       EffHisto->GetAxis(4)->SetRangeUser(1.,1.);    // reconstructed
1413       TH1D *all = WeightedProjection(EffHisto, axis, 3, axis_all);
1414
1415       EffHisto->GetAxis(4)->SetRangeUser(1.,1.);    // reconstructed
1416       TH1D *fake = WeightedProjection(EffHisto, axis, 1, &axis_fake);
1417       recc = (TH1D*) fake->Clone();
1418
1419       if(recc)
1420       {
1421         recc->Divide(fake,all,1,1,"B");
1422         recc->GetYaxis()->SetTitle("fake rate");
1423       }
1424   }
1425
1426   EffHisto->GetAxis(4)->SetRange(1,0);                          //Reset Range
1427
1428   if (recc) { // Coverity fix
1429     recc->SetName(name);
1430
1431     recc->GetXaxis()->SetTitle(fEffHisto->GetAxis(axis)->GetTitle());
1432
1433     snprintf(title,256,"%s vs %s",vsTitle, fEffHisto->GetAxis(axis)->GetTitle());
1434     recc->SetTitle(title);
1435
1436     if (axis == 2 ) recc->SetBit(TH1::kLogX);
1437   }
1438   return recc;
1439 }