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