]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Adding the task of Jitendra and Zaida for efficiency monitoring.
authorakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 14:35:23 +0000 (14:35 +0000)
committerakalweit <akalweit@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 17 Oct 2013 14:35:23 +0000 (14:35 +0000)
PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.cxx [new file with mode: 0644]
PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.h [new file with mode: 0644]
PWGPP/EvTrkSelection/AliSingleTrackEffCuts.cxx [new file with mode: 0644]
PWGPP/EvTrkSelection/AliSingleTrackEffCuts.h [new file with mode: 0644]
PWGPP/EvTrkSelection/macros/AddSingleTrackEfficiencyTask.C [new file with mode: 0644]
PWGPP/EvTrkSelection/macros/RebinCFContainer.C [new file with mode: 0644]
PWGPP/EvTrkSelection/macros/RunCFSingleTrackEfficiencyTask.C [new file with mode: 0644]

diff --git a/PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.cxx b/PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.cxx
new file mode 100644 (file)
index 0000000..79db627
--- /dev/null
@@ -0,0 +1,648 @@
+#include "TCanvas.h"
+#include "TParticle.h"
+#include "TH1I.h"
+#include <TDatabasePDG.h>
+
+#include "AliAnalysisDataSlot.h"
+#include "AliAnalysisDataContainer.h"
+#include "AliAnalysisManager.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+#include "AliAODTrack.h"
+#include "AliAODTracklets.h"
+#include "AliMultiplicity.h"
+#include "AliCFManager.h"
+#include "AliCFCutBase.h"
+#include "AliCFContainer.h"
+#include "AliESDEvent.h"
+#include "AliESDtrack.h"
+#include "AliESDtrackCuts.h"
+#include "AliESDVertex.h"
+#include "AliInputEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliVEvent.h"
+#include "AliAODEvent.h"
+#include "AliVVertex.h"
+#include "AliLog.h"
+#include "AliStack.h"
+#include "AliGenEventHeader.h"
+#include "AliAODMCHeader.h"
+
+#include "AliSingleTrackEffCuts.h"
+#include "AliCFSingleTrackEfficiencyTask.h"
+
+
+ClassImp(AliCFSingleTrackEfficiencyTask)
+
+//__________________________________________________________________________
+AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask() :
+  fReadAODData(0),
+  fCFManager(0x0),
+  fQAHistList(0x0),
+  fTrackCuts(0x0),
+  fMCCuts(0x0),
+  fTriggerMask(AliVEvent::kAny),
+  fSetFilterBit(kFALSE),
+  fbit(0),
+  fHistEventsProcessed(0x0)
+{
+  //
+  //Default constructor
+  //
+}
+
+//___________________________________________________________________________
+AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask(const Char_t* name,AliESDtrackCuts *trackcuts, AliSingleTrackEffCuts * mccuts) :
+  AliAnalysisTaskSE(name),
+  fReadAODData(0),
+  fCFManager(0x0),
+  fQAHistList(0x0),
+  fTrackCuts(trackcuts),
+  fMCCuts(mccuts),
+  fTriggerMask(AliVEvent::kAny),
+  fSetFilterBit(kFALSE),
+  fbit(0),
+  fHistEventsProcessed(0x0)
+{
+  //
+  // Constructor. Initialization of Inputs and Outputs
+  //
+  Info("AliCFSingleTrackEfficiencyTask","Calling Constructor");
+
+  // Output slot #1 writes into a TList container (nevents histogran)
+  DefineOutput(1,TH1I::Class());
+  // Output slot #2 writes into a TList container (distributions)
+  DefineOutput(2,AliCFContainer::Class());
+  // Output slot #3 writes the QA list
+  DefineOutput(3,TList::Class());
+  // Output slot #3 writes the ESD track cuts used
+  DefineOutput(4,AliESDtrackCuts::Class());
+  // Output slot #4 writes the particle and event selection object
+  DefineOutput(5,AliSingleTrackEffCuts::Class());
+}
+
+//_________________________________________________________________________________________________________________
+AliCFSingleTrackEfficiencyTask& AliCFSingleTrackEfficiencyTask::operator=(const AliCFSingleTrackEfficiencyTask& c) 
+{
+  //
+  // Assignment operator
+  //
+  if (this!=&c) {
+    AliAnalysisTaskSE::operator=(c) ;
+
+    fReadAODData = c.fReadAODData ;
+    fCFManager  = c.fCFManager;
+    fQAHistList = c.fQAHistList ;
+
+    if(c.fTrackCuts) { delete fTrackCuts; fTrackCuts = new AliESDtrackCuts(*(c.fTrackCuts)); }
+    if(c.fMCCuts) { delete fMCCuts; fMCCuts = new AliSingleTrackEffCuts(*(c.fMCCuts)); }
+    fTriggerMask = c.fTriggerMask;
+
+    fSetFilterBit  = c.fSetFilterBit;
+    fbit = c.fbit;
+    fHistEventsProcessed = c.fHistEventsProcessed;
+  }
+  return *this;
+}
+
+//________________________________________________________________________________________________________
+AliCFSingleTrackEfficiencyTask::AliCFSingleTrackEfficiencyTask(const AliCFSingleTrackEfficiencyTask& c) :
+  AliAnalysisTaskSE(c),
+  fReadAODData(c.fReadAODData),
+  fCFManager(c.fCFManager),
+  fQAHistList(c.fQAHistList),
+  fTrackCuts(c.fTrackCuts),
+  fMCCuts(c.fMCCuts),
+  fTriggerMask(c.fTriggerMask),
+  fSetFilterBit(c.fSetFilterBit),
+  fbit(c.fbit),
+  fHistEventsProcessed(c.fHistEventsProcessed)
+{
+  //
+  // Copy Constructor
+  //
+}
+
+//___________________________________________________________________________
+AliCFSingleTrackEfficiencyTask::~AliCFSingleTrackEfficiencyTask()
+{
+  //
+  // Destructor
+  //
+  Info("~AliCFSingleTrackEfficiencyTask","Calling Destructor");
+
+  if (fCFManager)           delete fCFManager;
+  if (fHistEventsProcessed) delete fHistEventsProcessed;
+  if (fQAHistList) { fQAHistList->Clear(); delete fQAHistList; }
+  if (fTrackCuts)           delete fTrackCuts;
+  if (fMCCuts)              delete fMCCuts;
+}
+
+//_________________________________________________
+void AliCFSingleTrackEfficiencyTask::Init() 
+{
+  //
+  // Initialization, checks + copy cuts
+  //
+  if(!fMCCuts) {
+    AliFatal(" MC Cuts not defined");
+    return;
+  }
+  if(!fTrackCuts) {
+    AliFatal(" Track Cuts not defined");
+    return;
+  }
+
+  AliESDtrackCuts* copyfTrackCuts = new AliESDtrackCuts(*fTrackCuts);
+  const char* nameoutput = GetOutputSlot(4)->GetContainer()->GetName();
+  copyfTrackCuts->SetName(nameoutput);
+  // Post the data
+  PostData(4,copyfTrackCuts);
+
+  AliSingleTrackEffCuts* copyfMCCuts = new AliSingleTrackEffCuts(*fMCCuts);
+  nameoutput = GetOutputSlot(5)->GetContainer()->GetName();
+  copyfMCCuts->SetName(nameoutput);
+  // Post the data
+  PostData(5,copyfMCCuts);
+
+  return;
+}
+
+//_________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::UserExec(Option_t *)
+{
+  //
+  // User Exec
+  //
+
+  Info("UserExec","") ;
+
+  AliVEvent* event = fInputEvent;
+
+  if(!fInputEvent) {
+    AliFatal("NO EVENT FOUND!");
+    return;
+  }
+  if(!fMCEvent) {
+    AliFatal("NO MC INFO FOUND");
+    return;
+  }
+
+  fHistEventsProcessed->Fill(0.5); // # of Event proceed        
+  Bool_t IsEventMCSelected = kFALSE;
+  Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
+
+  //
+  // Step 0: MC Gen Event Selection
+  //
+  if(isAOD) {
+    if(!event && AODEvent() && IsStandardAOD()) {
+      // In case there is an AOD handler writing a standard AOD, use the AOD 
+      // event in memory rather than the input (ESD) event.    
+      event = dynamic_cast<AliAODEvent*> (AODEvent());
+    }
+    IsEventMCSelected = fMCCuts->IsMCEventSelected(event);//AODs
+  } else {
+    IsEventMCSelected = fMCCuts->IsMCEventSelected(fMCEvent);//ESDs
+  }
+
+  // pass to the manager the event info to the cuts that need it 
+  if(!isAOD) fCFManager->SetMCEventInfo(fMCEvent);
+  else fCFManager->SetMCEventInfo(event);
+  fCFManager->SetRecEventInfo(event);
+
+  if(!IsEventMCSelected) {
+    AliDebug(3,"MC event not passing the quality criteria \n");
+    PostData(1,fHistEventsProcessed);
+    PostData(2,fCFManager->GetParticleContainer());
+    PostData(3,fQAHistList);
+    return;
+  }
+  fHistEventsProcessed->Fill(1.5); // # of Event after passing MC cuts
+       
+
+  //
+  // Step 1-3: Check the MC generated particles
+  // 
+  if(!isAOD) CheckESDMCParticles();
+  else CheckAODMCParticles();
+
+  //
+  // Step 4-7: Reconstructed event and track selection
+  //
+  Bool_t isRecoEventOk = fMCCuts->IsRecoEventSelected(event);
+
+  if(isRecoEventOk) {
+    fHistEventsProcessed->Fill(2.5); // # of Event after passing all cuts
+    CheckReconstructedParticles();
+  }
+  else AliDebug(3,"Event not passing quality criteria\n");
+
+
+  PostData(1,fHistEventsProcessed);
+  PostData(2,fCFManager->GetParticleContainer());
+  PostData(3,fQAHistList);
+
+  return;
+}
+
+
+//___________________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::Terminate(Option_t*)
+{
+  //
+  // Terminate
+  //
+
+  Info("Terminate","");
+  AliAnalysisTaskSE::Terminate();
+
+  /*
+  //draw some example plots....
+  AliCFContainer *cont= dynamic_cast<AliCFContainer*> (GetOutputData(2));
+       
+  TH1D* h00 =   cont->ShowProjection(5,0) ;
+  TH1D* h01 =   cont->ShowProjection(5,1) ;
+  TH1D* h02 =   cont->ShowProjection(5,2) ;
+  TH1D* h03 =   cont->ShowProjection(5,3) ;
+  TH1D* h04 =   cont->ShowProjection(5,4) ;
+  TH1D* h05 =   cont->ShowProjection(5,5) ;
+  TH1D* h06 =   cont->ShowProjection(5,6) ;
+       
+  h00->SetMarkerStyle(23) ;
+  h01->SetMarkerStyle(24) ;
+  h02->SetMarkerStyle(25) ;
+  h03->SetMarkerStyle(26) ;
+  h04->SetMarkerStyle(27) ;
+  h05->SetMarkerStyle(28) ;
+  h06->SetMarkerStyle(28) ;
+       
+       
+       
+  TCanvas * c =new TCanvas("c","",1400,800);
+  c->Divide(4,2);
+       
+  c->cd(1);
+  h00->Draw("p");
+  c->cd(2);
+  h01->Draw("p");
+  c->cd(3);
+  h02->Draw("p");
+  c->cd(4);
+  h03->Draw("p");
+  c->cd(5);
+  h04->Draw("p");
+  c->cd(6);
+  h05->Draw("p");
+  c->cd(7);
+  h06->Draw("p");
+       
+  c->SaveAs("plots.eps");
+  */
+       
+}
+
+
+//___________________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::UserCreateOutputObjects() 
+{
+  //
+  // UserCreateOutputObjects
+  //
+  
+  Info("CreateOutputObjects","CreateOutputObjects of task %s", GetName());
+
+  //slot #1
+  OpenFile(1);
+       
+  const char* nameoutput=GetOutputSlot(1)->GetContainer()->GetName();
+  //       fHistEventsProcessed = new TH1I(nameoutput"fHistEventsProcessed","fHistEventsProcessed",3,0,3) ;
+  fHistEventsProcessed = new TH1I(nameoutput,"fHistEventsProcessed",3,0,3) ;
+  fHistEventsProcessed->GetXaxis()->SetBinLabel(1,"All events");
+  fHistEventsProcessed->GetXaxis()->SetBinLabel(2,"Good MC events");
+  fHistEventsProcessed->GetXaxis()->SetBinLabel(3,"Good Reconstructed events");
+       
+  PostData(1,fHistEventsProcessed) ;
+  PostData(2,fCFManager->GetParticleContainer()) ;
+  PostData(3,fQAHistList) ;
+       
+  return;
+}
+
+
+//_________________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::CheckESDMCParticles()
+{
+  //
+  // Check ESD generated particles
+  //
+  if (!fMCEvent) {
+    AliFatal("NO MC INFO FOUND");
+    return;
+  }
+
+  Double_t containerInput[6] = { 0., 0., 0., 0., 0., 0. };
+
+  TArrayF vtxPos(3);
+  AliGenEventHeader *genHeader = NULL;
+  genHeader = fMCEvent->GenEventHeader();
+  genHeader->PrimaryVertex(vtxPos);
+  containerInput[4] = vtxPos[2]; // z-vtx position
+
+  Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
+  containerInput[5] = multiplicity; //reconstructed number of tracklets
+
+  // loop on the MC Generated Particles
+  for (Int_t ipart=0; ipart<fMCEvent->GetNumberOfTracks(); ipart++) {
+
+    AliMCParticle *mcPart = (AliMCParticle*)fMCEvent->GetTrack(ipart);
+    containerInput[0] = (Float_t)mcPart->Pt();
+    containerInput[1] = mcPart->Eta();
+    containerInput[2] = mcPart->Phi();
+    containerInput[3] = mcPart->Theta();
+
+    // Step 1. Particle passing through Generation criteria and filling
+    if( !fMCCuts->IsMCParticleGenerated(mcPart) ) {
+      AliDebug(3,"MC Particle not passing through genetations criteria\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCGenCut);
+
+    // Step 2. Particle passing through Kinematic criteria and filling
+    if( !fMCCuts->IsMCParticleInKineAcceptance(mcPart) ) {
+      AliDebug(3,"MC Particle not in the kine acceptance\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCKineCut);
+
+    // Step 3. Particle passing through Track ref criteria and filling
+    // did leave signal (enough clusters) on the detector
+    if( !fMCCuts->IsMCParticleInReconstructable(mcPart) ) {
+      AliDebug(3,"MC Particle not in the reconstructible\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCAccpCut);
+
+  }// end of particle loop
+
+  return;
+}
+
+
+//________________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::CheckAODMCParticles()
+{
+  //
+  // Check AOD generated particles
+  //
+  if (!fInputEvent) {
+    AliFatal("NO EVENT FOUND!");
+    return;
+  }
+
+  AliAODEvent* event = dynamic_cast<AliAODEvent*>(fInputEvent);
+  if(!event && AODEvent() && IsStandardAOD()) {
+    // In case there is an AOD handler writing a standard AOD, use the AOD 
+    // event in memory rather than the input (ESD) event.    
+    event = dynamic_cast<AliAODEvent*> (AODEvent());
+  }
+  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(event->FindListObject(AliAODMCParticle::StdBranchName()));
+  if (!mcArray) {
+    AliError("Could not find Monte-Carlo in AOD");
+    return;
+  }
+  AliAODMCHeader *mcHeader=NULL;
+  mcHeader = dynamic_cast<AliAODMCHeader*>(event->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+  if (!mcHeader) {
+    AliError("Could not find MC Header in AOD");
+    return;
+  }
+
+  Double_t containerInput[6] = { 0., 0., 0., 0., 0., 0. };
+  // Set the z-vertex position
+  containerInput[4]  = mcHeader->GetVtxZ();
+  // Multiplicity of the event defined as Ntracklets |eta|<1.0
+  Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
+  containerInput[5] = multiplicity;
+
+  for (Int_t ipart=0; ipart<mcArray->GetEntriesFast(); ipart++) {
+
+    AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle*>(mcArray->At(ipart));
+    if (!mcPart){
+      AliError("Failed casting particle from MC array!, Skipping particle");
+      continue;
+    }
+    containerInput[0] = (Float_t)mcPart->Pt();
+    containerInput[1] = mcPart->Eta();
+    containerInput[2] = mcPart->Phi();
+    containerInput[3] = mcPart->Theta();
+
+    // Step 1. Particle passing through Generation criteria and filling
+    if( !fMCCuts->IsMCParticleGenerated(mcPart) ) {
+      AliDebug(3,"MC Particle not passing quality criteria\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCGenCut);
+
+    // Step 2. Particle passing through Kinematic criteria and filling
+    if( !fMCCuts->IsMCParticleInKineAcceptance(mcPart) ) {
+      AliDebug(3,"MC Particle not in the acceptance\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCKineCut);
+
+    // Step 3. Particle passing through Track Ref criteria and filling
+    // but no info available for Track ref in AOD fillng same as above
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepMCAccpCut);
+  }
+
+  return;
+}
+
+
+//_______________________________________________________________________________
+AliESDtrack * AliCFSingleTrackEfficiencyTask::ConvertTrack(AliAODTrack *track)
+{
+  //
+  // Convert an AOD track to an  ESD track to apply ESDtrackCuts
+  //
+
+  AliAODEvent *aodEvent = (AliAODEvent*)fInputEvent;
+  const AliAODVertex *primary = aodEvent->GetPrimaryVertex();
+  Double_t pos[3],cov[6];
+  primary->GetXYZ(pos);
+  primary->GetCovarianceMatrix(cov);
+  const AliESDVertex vESD(pos,cov,100.,100);
+
+  AliESDtrack *esdTrack =  new AliESDtrack(track);
+  // set the TPC cluster info
+  esdTrack->SetTPCClusterMap(track->GetTPCClusterMap());
+  esdTrack->SetTPCSharedMap(track->GetTPCSharedMap());
+  esdTrack->SetTPCPointsF(track->GetTPCNclsF());
+  // needed to calculate the impact parameters
+  esdTrack->RelateToVertex(&vESD,0.,3.);
+  //  std::cout << " primary vtx "<< primary << std::endl;
+  //  std::cout << " esdvtx "<< vESD.GetName() << std::endl;
+  //  std::cout<< " esdtrack pt "<< esdTrack.Pt() << " and status " << esdTrack.GetStatus() <<endl;
+  //  std::cout << " aod track "<< track<< " and status " << track->GetStatus() << std::endl;
+  //  std::cout << " esd track "<< esdTrack<< " and status " << esdTrack->GetStatus() << std::endl;
+  return esdTrack;
+}
+
+
+//___________________________________________________________________________
+void AliCFSingleTrackEfficiencyTask::CheckReconstructedParticles()
+{
+  //
+  // Check reconstructed particles
+  //
+
+  AliVEvent* event = fInputEvent;
+  Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
+  if(!event && AODEvent() && IsStandardAOD()) {
+   event  = dynamic_cast<AliAODEvent*> (AODEvent());
+  }
+
+  Double_t containerInput[6] = { 0, 0, 0, 0, 0, 0};   // contains reconstructed quantities
+  Double_t containerInputMC[6] = { 0, 0, 0, 0, 0, 0}; // contains generated quantities
+  
+  const AliVVertex *vertex = event->GetPrimaryVertex();
+  // set the z-vertex position
+  containerInput[4] = vertex->GetZ();
+  containerInputMC[4] = containerInput[4];
+  // set the event multiplicity as Ntracklets in |eta|<1.0
+  Double_t multiplicity = (Double_t)GetNumberOfTrackletsInEtaRange(-1.0,1.0);
+  containerInput[5] = multiplicity;
+  containerInputMC[5] = multiplicity;
+
+  // Reco tracks track loop
+  AliVParticle* track = NULL;
+  for (Int_t iTrack = 0; iTrack<event->GetNumberOfTracks(); iTrack++) {
+
+    track = event->GetTrack(iTrack);
+    if(!track) {
+      AliDebug(3,Form("Track %d not found",iTrack));
+      continue;
+    }
+
+    Double_t mom[3];
+    track->PxPyPz(mom);
+    Double_t pt = TMath::Sqrt(mom[0]*mom[0]+mom[1]*mom[1]);
+    containerInput[0] = pt;
+    containerInput[1] = track->Eta();
+    containerInput[2] = track->Phi();
+    containerInput[3] = track->Theta();
+
+    //
+    // Step 4. Track that are recostructed, filling
+    //
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepReconstructed);
+
+    //
+    // Step 5. Track that are recostructed and pass acceptance cuts, filling
+    //
+    if (!fMCCuts->IsRecoParticleKineAcceptance(track)) continue;
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoKineCuts);
+
+    // is track associated to particle ? if yes + implimenting the physical primary..
+    Int_t label = TMath::Abs(track->GetLabel());
+    if (label<=0) {
+      AliDebug(3,"Particle not matching MC label \n");
+      continue;
+    }
+
+    //
+    // Step 6-7. Track that are recostructed + Kine + Quality criteria, filling
+    //
+
+    // check particle selections at MC level
+    AliVParticle *mcPart  = (AliVParticle*)fMCEvent->GetTrack(label);
+    if(!mcPart) continue;
+    containerInputMC[0] = (Float_t)mcPart->Pt();
+    containerInputMC[1] = mcPart->Eta();
+    containerInputMC[2] = mcPart->Phi();
+    containerInputMC[3] = mcPart->Theta();
+
+    if (!fMCCuts->IsMCParticleGenerated(mcPart)) continue;
+
+
+    // for filter bit selection
+    AliAODTrack *aodTrack = dynamic_cast<AliAODTrack*>(track);
+    if(isAOD && fSetFilterBit) if (!aodTrack->TestFilterMask(BIT(fbit))) continue;
+
+    Bool_t isESDtrack = track->IsA()->InheritsFrom("AliESDtrack");
+    AliESDtrack *tmptrack = NULL;
+    if(isESDtrack) {
+      tmptrack = dynamic_cast<AliESDtrack*>(track); // ESDs
+    } else {
+      tmptrack = ConvertTrack(aodTrack); // AODs
+    }
+
+    // exclude global constrained and TPC only tracks (filter bits 128 and 512)
+    Int_t id = tmptrack->GetID();
+    if(isAOD && id<0) {
+      AliDebug(3,"Track removed bc corresponts to either filter bit 128 or 512 (TPC only tracks)\n");
+      continue;
+    }
+
+    // Apply ESD track cuts
+    if( !fTrackCuts->IsSelected(tmptrack) ){
+      AliDebug(3,"Reconstructed track not passing quality criteria\n");
+      if(isAOD) delete tmptrack;
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInputMC,kStepReconstructedMC);
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoQualityCuts);
+    if(isAOD) delete tmptrack;
+
+    //
+    // Step8, PID check
+    //
+    if( !fMCCuts->IsRecoParticlePID(track) ){
+      AliDebug(3,"Reconstructed track not passing PID criteria\n");
+      continue;
+    }
+    fCFManager->GetParticleContainer()->Fill(containerInput,kStepRecoPID);
+
+  }
+  return;
+}
+
+//______________________________________________________________________
+Int_t AliCFSingleTrackEfficiencyTask::GetNumberOfTrackletsInEtaRange(Double_t mineta, Double_t maxeta)
+{
+  //
+  // counts tracklets in given eta range
+  //
+
+  AliAODEvent* event = dynamic_cast<AliAODEvent*>(fInputEvent);
+  Bool_t isAOD = fInputEvent->IsA()->InheritsFrom("AliAODEvent");
+  if(!event && AODEvent() && IsStandardAOD()) {
+   event  = dynamic_cast<AliAODEvent*> (AODEvent());
+  }
+
+  Int_t count=0;
+
+  if(isAOD) {
+    AliAODTracklets* tracklets = event->GetTracklets();
+    Int_t nTr=tracklets->GetNumberOfTracklets();
+    for(Int_t iTr=0; iTr<nTr; iTr++){
+      Double_t theta=tracklets->GetTheta(iTr);
+      Double_t eta=-TMath::Log(TMath::Tan(theta/2.));
+      if(eta>mineta && eta<maxeta) count++;
+    }
+  } else {  
+    AliESDEvent *esdEvent = dynamic_cast<AliESDEvent*>(fInputEvent);
+    const AliMultiplicity *mult = esdEvent->GetMultiplicity();
+    if (mult) {
+      if (mult->GetNumberOfTracklets()>0) {    
+       for (Int_t n=0; n<mult->GetNumberOfTracklets(); n++) {
+         Double_t eta = -TMath::Log( TMath::Tan( mult->GetTheta(n) / 2.0 ) );
+         if(eta>mineta && eta<maxeta) count++;
+       }
+      }
+    }
+  }
+
+  return count;
+}
diff --git a/PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.h b/PWGPP/EvTrkSelection/AliCFSingleTrackEfficiencyTask.h
new file mode 100644 (file)
index 0000000..f8c6802
--- /dev/null
@@ -0,0 +1,101 @@
+#ifndef AliCFSINGLETRACKEFFICIENCYTASK_H
+#define AliCFSINGLETRACKEFFICIENCYTASK_H
+
+#include "AliAnalysisTaskSE.h"
+#include "AliSingleTrackEffCuts.h"
+
+class TH1I;
+class TParticle;
+class AliStack;
+class AliCFManager;
+class AliAODTrack;
+class AliESDtrack;
+class AliESDtrackCuts;
+class AliESDVertex;
+class AliVVertex;
+class AliVParticle;
+
+
+class AliCFSingleTrackEfficiencyTask : public AliAnalysisTaskSE {
+
+ public:
+
+  enum {
+    kStepMCGenCut         = 0, // selected generated particles, event selection
+    kStepMCKineCut        = 1, // generated particles after acceptance cuts
+    kStepMCAccpCut        = 2, // particles with a minimum number of clusters in detector (ESD only)
+    kStepReconstructed    = 3, // reconstructed tracks
+    kStepRecoKineCuts     = 4, // reconstructed tracks after acceptance
+    kStepReconstructedMC  = 5, // tracks passing ESD+MC trackCuts (kine properties)
+    kStepRecoQualityCuts  = 6, // tracks passing ESD+MC trackCuts (reco properties)
+    kStepRecoPID          = 7  // tracks after PID criteria
+  };
+
+  AliCFSingleTrackEfficiencyTask();
+  AliCFSingleTrackEfficiencyTask(const Char_t* name, AliESDtrackCuts *trackcuts, AliSingleTrackEffCuts * mccuts);
+  AliCFSingleTrackEfficiencyTask(const AliCFSingleTrackEfficiencyTask& c);
+  AliCFSingleTrackEfficiencyTask& operator= (const AliCFSingleTrackEfficiencyTask& c);
+  virtual ~AliCFSingleTrackEfficiencyTask();
+
+  // ANALYSIS FRAMEWORK STUFF to loop on data and fill output objects
+  virtual void UserCreateOutputObjects();
+  virtual void UserExec(Option_t *option);
+  virtual void Terminate(Option_t *);
+  virtual void Init();
+  virtual void LocalInit() { Init(); }
+
+  // CORRECTION FRAMEWORK RELATED FUNCTIONS
+  void           SetCFManager(AliCFManager* io) { fCFManager = io; }   // global correction manager
+  AliCFManager * GetCFManager() const { return fCFManager; }           // get corr manager
+  void           SetQAList(TList* list) { fQAHistList = list; }        // set CF QA list (empty, not used, but needed)
+
+  // Setters
+  // analyze AOD:1 or ESD:0 data
+  void SetReadAODData (Bool_t flag=kTRUE) { fReadAODData=flag; }
+  // select track filter bit:1 or not:0 and set the bit
+  void SetFilterBit (Bool_t flag=kTRUE) { fSetFilterBit=flag; }
+  void SetFilterType (Int_t fbittype) { fbit=fbittype; }
+  // select trigger event mask
+  void SetTriggerMask(ULong64_t mask=0) { fTriggerMask=mask; }
+
+  // Getters
+  // analyze AOD:1 or ESD:0 data
+  Bool_t IsReadAODData()   const { return fReadAODData; }
+  // select trigger event mask
+  ULong64_t GetTriggerMask() { return fTriggerMask; }
+  // reconstructed track cuts
+  AliESDtrackCuts *GetTrackCuts(){ return (AliESDtrackCuts*)fTrackCuts; } 
+  // particle and event cuts 
+  AliSingleTrackEffCuts *GetSingleTrackEffCuts() { return (AliSingleTrackEffCuts*)fMCCuts; }
+
+
+ protected:
+
+  // Check ESD generated particles
+  void CheckESDMCParticles();
+  // Check AOD generated particles
+  void CheckAODMCParticles();
+  // Check reconstructed particles
+  void CheckReconstructedParticles();
+  // Convert AOD track to an ESD track
+  AliESDtrack* ConvertTrack(AliAODTrack *track);
+  // Count the number of tracklets in given eta range
+  Int_t GetNumberOfTrackletsInEtaRange(Double_t mineta, Double_t maxeta);
+
+  Bool_t fReadAODData;       // flag for AOD/ESD input files
+  AliCFManager *fCFManager;  // pointer to the CF manager slot 2
+  TList *fQAHistList;        // list of QA histograms slot 3
+
+  AliESDtrackCuts *fTrackCuts;    // track cuts (reconstructed level) slot 4
+  AliSingleTrackEffCuts *fMCCuts; // particle cuts used slot 5
+  ULong64_t fTriggerMask;         // event selection trigger mask
+
+  Bool_t fSetFilterBit; // flag to decide if applying filter-bit selection to tracks
+  Int_t  fbit;          // filter-bit selection to tracks
+
+  TH1I  *fHistEventsProcessed;   //! histo for monitoring the number of events processed slot 1
+
+  ClassDef(AliCFSingleTrackEfficiencyTask,1)
+};
+
+#endif
diff --git a/PWGPP/EvTrkSelection/AliSingleTrackEffCuts.cxx b/PWGPP/EvTrkSelection/AliSingleTrackEffCuts.cxx
new file mode 100644 (file)
index 0000000..14b69cd
--- /dev/null
@@ -0,0 +1,634 @@
+#include "AliVEvent.h"
+#include "AliMCEvent.h"
+#include "TParticle.h"
+#include "AliVParticle.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliVVertex.h"
+#include "AliLog.h"
+#include "AliGenEventHeader.h"
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliAODMCHeader.h"
+#include "AliAODMCParticle.h"
+#include "AliMCParticle.h"
+#include "AliVParticle.h"
+#include "AliStack.h"
+#include "AliMCEventHandler.h"
+#include "AliPIDResponse.h"
+
+#include "AliSingleTrackEffCuts.h"
+
+
+ClassImp(AliSingleTrackEffCuts)
+
+//_____________________________________________________
+AliSingleTrackEffCuts::AliSingleTrackEffCuts():
+AliAnalysisCuts(),
+  fisAOD(kTRUE),
+  fIsPdgCode(kFALSE),
+  fPdgCode(0),
+  fEtaMin(-12),
+  fEtaMax(12),
+  fYMin(-12),
+  fYMax(12),
+  fPtMin(-15),
+  fPtMax(15),
+  fIsCharged(kTRUE),
+  fTriggerMask(AliVEvent::kAny),
+  fMinVtxType(0),
+  fMinVtxContr(1),
+  fMaxVtxZ(10.),
+  fCutOnZVertexSPD(0),
+  fnClusITS(0),
+  fnClusTPC(0),
+  fnClusTOF(0),
+  fnClusMUON(0),
+  fusePid(kFALSE),
+  fParticlePid(AliPID::kPion),
+  fuseTPCPid(kTRUE),
+  fnPTPCBins(0),
+  fnPTPCBinLimits(0),
+  fPTPCBinLimits(0),
+  fnSigmaTPC(0),
+  fPmaxTPC(9999.),
+  fuseTOFPid(kTRUE),
+  fnPTOFBins(0),
+  fnPTOFBinLimits(0),
+  fPTOFBinLimits(0),
+  fnSigmaTOF(0),
+  fPmaxTOF(9999.)
+{
+  //
+  // Default constructor
+  //
+}
+
+//________________________________________________________________________________
+AliSingleTrackEffCuts::AliSingleTrackEffCuts(const char* name, const char* title):
+AliAnalysisCuts(name,title),
+  fisAOD(kTRUE),
+  fIsPdgCode(kFALSE),
+  fPdgCode(0),
+  fEtaMin(-12),
+  fEtaMax(12),
+  fYMin(-12),
+  fYMax(12),
+  fPtMin(-15),
+  fPtMax(15),
+  fIsCharged(kTRUE),
+  fTriggerMask(AliVEvent::kAny),
+  fMinVtxType(0),
+  fMinVtxContr(1),
+  fMaxVtxZ(10.),
+  fCutOnZVertexSPD(0),
+  fnClusITS(0),
+  fnClusTPC(0),
+  fnClusTOF(0),
+  fnClusMUON(0),
+  fusePid(kFALSE),
+  fParticlePid(AliPID::kPion),
+  fuseTPCPid(kTRUE),
+  fnPTPCBins(0),
+  fnPTPCBinLimits(0),
+  fPTPCBinLimits(0),
+  fnSigmaTPC(0),
+  fPmaxTPC(9999.),
+  fuseTOFPid(kTRUE),
+  fnPTOFBins(0),
+  fnPTOFBinLimits(0),
+  fPTOFBinLimits(0),
+  fnSigmaTOF(0),
+  fPmaxTOF(9999.)
+{
+  //
+  // Default constructor
+  //
+}
+
+//_________________________________________________________________________________
+AliSingleTrackEffCuts::AliSingleTrackEffCuts(const AliSingleTrackEffCuts &source):
+  AliAnalysisCuts(source),
+  fisAOD(source.fisAOD),
+  fIsPdgCode(source.fIsPdgCode),
+  fPdgCode(source.fPdgCode),
+  fEtaMin(source.fEtaMin),
+  fEtaMax(source.fEtaMax),
+  fYMin(source.fYMin),
+  fYMax(source.fYMax),
+  fPtMin(source.fPtMin),
+  fPtMax(source.fPtMax),
+  fIsCharged(source.fIsCharged),
+  fTriggerMask(source.fTriggerMask),
+  fMinVtxType(source.fMinVtxType),
+  fMinVtxContr(source.fMinVtxContr),
+  fMaxVtxZ(source.fMaxVtxZ),
+  fCutOnZVertexSPD(source.fCutOnZVertexSPD),
+  fnClusITS(source.fnClusITS),
+  fnClusTPC(source.fnClusTPC),
+  fnClusTOF(source.fnClusTOF),
+  fnClusMUON(source.fnClusMUON),
+  fusePid(source.fusePid),
+  fParticlePid(source.fParticlePid),
+  fuseTPCPid(source.fuseTPCPid),
+  fnPTPCBins(source.fnPTPCBins),
+  fnPTPCBinLimits(source.fnPTPCBinLimits),
+  fPTPCBinLimits(source.fPTPCBinLimits),
+  fnSigmaTPC(source.fnSigmaTPC),
+  fPmaxTPC(source.fPmaxTPC),
+  fuseTOFPid(source.fuseTOFPid),
+  fnPTOFBins(source.fnPTOFBins),
+  fnPTOFBinLimits(source.fnPTOFBinLimits),
+  fPTOFBinLimits(source.fPTOFBinLimits),
+  fnSigmaTOF(source.fnSigmaTOF),
+  fPmaxTOF(source.fPmaxTOF)
+{
+  //
+  // Copy constructor
+  //
+}
+
+
+//_________________________________________________________________________________________
+AliSingleTrackEffCuts &AliSingleTrackEffCuts::operator=(const AliSingleTrackEffCuts &source)
+{
+  //
+  // assignment operator
+  //
+  if(&source == this) return *this;
+
+  fisAOD = source.fisAOD;
+
+  fIsPdgCode = source.fIsPdgCode;
+  fPdgCode = source.fPdgCode;
+  fEtaMin = source.fEtaMin;
+  fEtaMax = source.fEtaMax;
+  fYMin = source.fYMin;
+  fYMax = source.fYMax;
+  fPtMin = source.fPtMin;
+  fPtMax = source.fPtMax;
+  fIsCharged = source.fIsCharged;
+
+  fTriggerMask = source.fTriggerMask;
+  fMinVtxType = source.fMinVtxType;
+  fMinVtxContr = source.fMinVtxContr;
+  fMaxVtxZ = source.fMaxVtxZ;
+  fCutOnZVertexSPD = source.fCutOnZVertexSPD;
+
+  fnClusITS = source.fnClusITS;
+  fnClusTPC = source.fnClusTPC;
+  fnClusTOF = source.fnClusTOF;
+  fnClusMUON = source.fnClusMUON;
+
+  fusePid = source.fusePid;
+  fParticlePid = source.fParticlePid;
+  fuseTPCPid = source.fuseTPCPid;
+  fnPTPCBins = source.fnPTPCBins;
+  fnPTPCBinLimits = source.fnPTPCBinLimits;
+  if(source.fPTPCBinLimits && source.fnSigmaTPC) 
+    SetTPCSigmaPtBins(source.fnPTPCBins,source.fPTPCBinLimits,source.fnSigmaTPC);
+  fPmaxTPC = source.fPmaxTPC;
+  fuseTOFPid = source.fuseTOFPid;
+  fnPTOFBins = source.fnPTOFBins;
+  fnPTOFBinLimits = source.fnPTOFBinLimits;
+  if(source.fPTOFBinLimits && source.fnSigmaTOF) 
+    SetTOFSigmaPtBins(source.fnPTOFBins,source.fPTOFBinLimits,source.fnSigmaTOF);
+  fPmaxTOF = source.fPmaxTOF;
+
+  return *this;
+}
+
+//______________________________________________
+AliSingleTrackEffCuts::~AliSingleTrackEffCuts()
+{
+  //
+  // Destructor
+  //
+
+  if(fPTPCBinLimits){
+    delete [] fPTPCBinLimits;
+    fPTPCBinLimits=NULL;
+  }
+  if(fnSigmaTPC){
+    delete [] fnSigmaTPC;
+    fnSigmaTPC=NULL;
+  }
+  if(fPTOFBinLimits){
+    delete [] fPTOFBinLimits;
+    fPTOFBinLimits=NULL;
+  }
+  if(fnSigmaTOF){
+    delete [] fnSigmaTOF;
+    fnSigmaTOF=NULL;
+  }
+
+}
+
+//______________________________________________
+Bool_t AliSingleTrackEffCuts::IsMCEventSelected(TObject* obj)
+{
+  //
+  // Event selection at MC level 
+  //  |zvtx|<=fMaxVtxZ
+  //
+  if(!obj) return  kFALSE;
+  Bool_t isSelected = kTRUE;
+  
+  AliMCEvent *event=0;
+  AliGenEventHeader *genHeader=0; //ESDs
+  AliAODMCHeader *mcHeader=0; // AODs
+  Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODEvent");
+
+  if(!isAOD) {
+    event = dynamic_cast<AliMCEvent*>(obj);
+    if (!event) return  kFALSE;
+    genHeader = event->GenEventHeader();
+    if (!genHeader) return  kFALSE;
+  } else {
+    AliAODEvent *aodEvent = dynamic_cast<AliAODEvent*> (obj);
+    if (!aodEvent) return  kFALSE;
+    mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
+    if (!mcHeader) {
+      AliError("Could not find MC Header in AOD");
+      return kFALSE;
+    }
+
+    // Check for not null trigger mask to evict pPb MC buggy-events
+    Int_t runnumber = aodEvent->GetRunNumber();
+    if(aodEvent->GetTriggerMask()==0 && 
+       (runnumber>=195344 && runnumber<=195677)){
+      AliDebug(3,"Event rejected because of null trigger mask");
+      return kFALSE;
+    }
+  }
+
+
+  // Cut on Z Vertex ( |z-vtx| <= fMaxVtxZ )
+  TArrayF vtxPos(3);
+  Double_t zMCVertex=-1000;
+  if(!isAOD) {
+    genHeader->PrimaryVertex(vtxPos);
+    zMCVertex = vtxPos[2];
+  } else {
+    zMCVertex = mcHeader->GetVtxZ();
+  }  
+  if( TMath::Abs(zMCVertex)>fMaxVtxZ ) isSelected = kFALSE;
+
+
+  return isSelected;
+}
+
+//__________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsMCParticleGenerated(TObject* obj)
+{
+  //
+  // Check generated particles (primary, charged, pdgcode) 
+  //
+  
+  if(!obj) return kFALSE;
+  if(!obj->InheritsFrom("AliVParticle")) {
+    AliError("object must derived from AliVParticle !");
+    return kFALSE;
+  }
+  AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
+  if(!particle) return kFALSE;
+
+  Bool_t isSelected = kTRUE;
+
+  // Check particle Pdg Code
+  if(fIsPdgCode && TMath::Abs( particle->PdgCode() )!= fPdgCode) isSelected = kFALSE;
+
+  // Charge selection
+  if(fIsCharged && !(particle->Charge()!=0)) isSelected = kFALSE;
+
+  // Selection of Physical Primary particles
+  if(!fisAOD) { // check on ESDs
+    AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());          
+    AliMCEvent* mcevent = mcinfo->MCEvent();
+    AliStack* stack = ((AliMCEvent*)mcevent)->Stack();
+    AliMCParticle* mcPart = dynamic_cast<AliMCParticle *>(obj);             
+    if(!stack->IsPhysicalPrimary(mcPart->GetLabel())) {
+      isSelected = kFALSE;
+    }
+  } else { // Check on AODs
+    AliAODMCParticle* mcPart = dynamic_cast<AliAODMCParticle *>(obj);      
+    if(!mcPart->IsPhysicalPrimary()) isSelected = kFALSE;
+  }
+
+  return isSelected;
+}
+
+
+//_____________________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsMCParticleInKineAcceptance(TObject *obj)
+{
+  //
+  // Check generated particles (eta, y, pt)
+  //
+
+  if(!obj) return  kFALSE;
+  if(!obj->InheritsFrom("AliVParticle")) AliError("object must derived from AliVParticle !");
+  AliVParticle* particle = dynamic_cast<AliVParticle *>(obj);
+  if(!particle) return kFALSE;
+
+  Bool_t isSelected = kTRUE;
+
+  // Cut on eta
+  if(particle->Eta()<fEtaMin || particle->Eta()>fEtaMax) isSelected = kFALSE;
+
+  // Cut on y
+  Double_t energy = particle->E();
+  Double_t pz = particle->Pz();
+  Double_t particleY = (TMath::Abs(energy-pz)>1e-10) ?  0.5*TMath::Log( (energy+pz)/(energy-pz) ) : 1e6;
+  if(particleY<fYMin || particleY>fYMax) isSelected = kFALSE;
+
+  // Cut on pt
+  if(particle->Pt()<fPtMin || particle->Pt()>fPtMax) isSelected = kFALSE;
+
+  return isSelected;
+}
+
+
+//_______________________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsMCParticleInReconstructable(TObject *obj)
+{
+  //
+  // Check if particle has left enough hits in the detectors (only at ESD level)
+  //
+
+  if(!obj) return kFALSE;
+  TString className(obj->ClassName());
+  if (className.CompareTo("AliMCParticle") != 0) {
+    AliError("obj must point to an AliMCParticle !");
+    return kTRUE; // <===================================== FIX ME !!
+  }
+
+  AliMCParticle * part = dynamic_cast<AliMCParticle*>(obj);
+  if(!part) return kFALSE;
+
+  Bool_t isSelected = kTRUE;
+
+  Int_t nHitsITS=0, nHitsTPC=0, nHitsTRD=0, nHitsTOF=0, nHitsMUON=0;
+  for(Int_t iTrackRef=0; iTrackRef<part->GetNumberOfTrackReferences(); iTrackRef++) {
+    AliTrackReference * trackRef = part->GetTrackReference(iTrackRef);
+    if(trackRef){
+      Int_t detectorId = trackRef->DetectorId();
+      switch(detectorId) {
+      case AliTrackReference::kITS  : nHitsITS++  ;
+      case AliTrackReference::kTPC  : nHitsTPC++  ;
+      case AliTrackReference::kTRD  : nHitsTRD++  ; 
+      case AliTrackReference::kTOF  : nHitsTOF++  ;
+      case AliTrackReference::kMUON : nHitsMUON++ ; 
+      }
+    }
+  }
+
+  if(nHitsITS<fnClusITS) isSelected = kFALSE;
+  if(nHitsTPC<fnClusTPC) isSelected = kFALSE;
+  if(nHitsTOF<fnClusTOF) isSelected = kFALSE;
+  if(nHitsMUON<fnClusMUON) isSelected = kFALSE;
+
+  return isSelected;
+}
+
+
+//_____________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsRecoEventSelected(TObject* obj)
+{
+  //
+  // Event selection at reconstructed level (trigger, zvtx)
+  //
+  AliVEvent *event = dynamic_cast<AliVEvent*>(obj);
+  if(!event) return kFALSE;
+
+  Bool_t isSelected = kTRUE;
+
+  // Check if event is accepted by the Physics selection
+  UInt_t trigFired = ((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected();
+  Bool_t isEvtSelected = (trigFired & fTriggerMask);
+  if(!isEvtSelected) isSelected = kFALSE;
+
+  // Vertex selection
+  Bool_t isVtxSelected = IsVertexSelected(event);
+  if(!isVtxSelected) isSelected = kFALSE;
+
+  return isSelected;
+}
+
+
+//______________________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsRecoParticleKineAcceptance(TObject *obj)
+{
+  //
+  // Check if reconstructed particle is in the acceptance (eta, y, pt)
+  //
+
+  Bool_t isSelected = kTRUE;
+
+  AliVParticle *track = dynamic_cast<AliVParticle*>(obj);
+  // Cut on eta
+  if(track->Eta()<fEtaMin || track->Eta()>fEtaMax) isSelected = kFALSE;
+
+  // Cut on y
+  if(track->Y()<fYMin || track->Y()>fYMax) isSelected = kFALSE;
+
+  // Cut on pt
+  if(track->Pt()<fPtMin || track->Pt() >fPtMax) isSelected = kFALSE;
+
+  return isSelected;
+}
+
+//_______________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsVertexSelected(AliVEvent *event)
+{
+  //
+  // Check if the reconstructed vertex is selected
+  //
+
+  Bool_t accept = kTRUE;
+  Bool_t isAOD = event->IsA()->InheritsFrom("AliAODEvent");
+
+  const AliVVertex *vertex = event->GetPrimaryVertex();
+  if(!vertex){
+    accept = kFALSE;
+    AliInfo("no vtx");
+    return accept;
+  }
+
+  // Cut on vertex type  
+  TString title=vertex->GetTitle();
+  if(title.Contains("Z") && fMinVtxType>1){
+    accept=kFALSE;
+  } else if(title.Contains("3D") && fMinVtxType>2){
+    accept=kFALSE;
+  }
+  
+  // cut on minimum number of contributors
+  if(vertex->GetNContributors()<fMinVtxContr){
+    AliInfo(Form("too few contributors %d",vertex->GetNContributors()));
+    accept=kFALSE;
+  }
+
+  // cut on absolute |z| of the vertex
+  if(TMath::Abs(vertex->GetZ())>fMaxVtxZ) {
+    AliInfo("outside the Vtx range");
+    accept=kFALSE;
+  }
+
+  // cut on distance of SPD and TRK vertexes
+  if(fCutOnZVertexSPD==1) {
+    const AliVVertex *vSPD = NULL;
+    if(isAOD) {
+      vSPD = ((AliAODEvent*)event)->GetPrimaryVertexSPD(); 
+    }else {    
+      vSPD = ((AliESDEvent*)event)->GetPrimaryVertexSPD();  
+    }
+
+    if(vSPD && vSPD->GetNContributors()>=fMinVtxContr) {      
+      if(TMath::Abs(vSPD->GetZ()-vertex->GetZ())>0.5) accept = kFALSE;
+    }
+  }
+
+  return accept;
+}
+
+//_________________________________________________________________________________________________
+void AliSingleTrackEffCuts::SetTPCSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
+{
+  //
+  // Set TPC Pid number-of-Sigma in P bins
+  //
+
+  // Set the pt bins
+  if(fPTPCBinLimits) {
+    delete [] fPTPCBinLimits;
+    fPTPCBinLimits = NULL;
+    printf("Changing the TPC cut P bins\n");
+  }
+  if(fnSigmaTPC) {
+    delete [] fnSigmaTPC;
+    fnSigmaTPC = NULL;
+    printf("Changing the TPC Sigma cut per p bin\n");
+  }
+
+  fnPTPCBins = nPtBins;
+  fnPTPCBinLimits = nPtBins+1;
+  fPTPCBinLimits = new Float_t[fnPTPCBinLimits];
+  for(Int_t ib=0; ib<nPtBins+1; ib++) fPTPCBinLimits[ib]=pBinLimits[ib];
+  fnSigmaTPC = new Float_t[fnPTPCBins];
+  for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTPC[ib]=sigmaBin[ib];
+
+  return;
+}
+//_____________________________________________________________________________________________________
+void AliSingleTrackEffCuts::SetTOFSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin)
+{
+  //
+  // Set TOF Pid number-of-Sigma in P bins
+  //
+
+  // Set the pt bins
+  if(fPTOFBinLimits) {
+    delete [] fPTOFBinLimits;
+    fPTOFBinLimits = NULL;
+    printf("Changing the TOF cut P bins\n");
+  }
+  if(fnSigmaTOF) {
+    delete [] fnSigmaTOF;
+    fnSigmaTOF = NULL;
+    printf("Changing the TOF Sigma cut per p bin\n");
+  }
+
+  fnPTOFBins = nPtBins;
+  fnPTOFBinLimits = nPtBins+1;
+  fPTOFBinLimits = new Float_t[fnPTOFBinLimits];
+  for(Int_t ib=0; ib<nPtBins+1; ib++) fPTOFBinLimits[ib]=pBinLimits[ib];
+  fnSigmaTOF = new Float_t[fnPTOFBins];
+  for(Int_t ib=0; ib<nPtBins; ib++) fnSigmaTOF[ib]=sigmaBin[ib];
+
+  return;
+}
+
+//___________________________________________________________________________
+Bool_t AliSingleTrackEffCuts::CheckTPCPIDStatus(AliAODTrack *track) const{
+  //
+  // Check TPC PID status
+  //
+
+  if ((track->GetStatus() & AliESDtrack::kTPCin)==0) return kFALSE;
+  UShort_t nTPCClus=track->GetTPCClusterMap().CountBits();
+  if (nTPCClus<70) return kFALSE;
+  return kTRUE;
+}
+
+//___________________________________________________________________________
+Bool_t AliSingleTrackEffCuts::CheckTOFPIDStatus(AliAODTrack *track) const{
+  //
+  // Check TOC PID status
+  //
+
+  if ((track->GetStatus()&AliESDtrack::kTOFout)==0)   return kFALSE;
+  if ((track->GetStatus()&AliESDtrack::kTIME)==0)     return kFALSE;
+  if ((track->GetStatus()&AliESDtrack::kTOFpid)==0)   return kFALSE;
+  if (!(track->GetStatus()&AliESDtrack::kTOFmismatch)==0) return kFALSE;
+  return kTRUE;
+}
+
+//___________________________________________________________________________
+Bool_t AliSingleTrackEffCuts::IsRecoParticlePID(TObject *obj)
+{
+  //
+  // Check Particle PID (AOD only for now!)
+  //
+
+  Bool_t isSelected = kFALSE;
+  Bool_t isAOD = obj->IsA()->InheritsFrom("AliAODTrack");
+
+  if(!isAOD || !GetUsePid()) return isSelected;
+  if(!fuseTPCPid && !fuseTOFPid) return isSelected;
+
+  AliAODTrack *track = dynamic_cast<AliAODTrack*>(obj);
+  if(!track) { cout<<"No track found"<<endl; return isSelected; }
+
+  // AliAODPid *pid = track->GetDetPid();
+  // if(!pid) { cout<<"No AliAODPid found"<<endl; return isSelected; }
+
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  AliInputEventHandler *inputHandler=(AliInputEventHandler*)mgr->GetInputEventHandler();
+  AliPIDResponse *pidResp=inputHandler->GetPIDResponse();
+  if(!pidResp) { cout<<"No PidResponse found"<<endl; return isSelected;}
+
+  Double_t pPart = track->P();
+  
+  // Check detector status
+  Bool_t okTPC = CheckTPCPIDStatus(track);
+  Bool_t okTOF = CheckTOFPIDStatus(track);
+
+  //  Check Number of Sigmas
+  Double_t nsigmaTPC=pidResp->NumberOfSigmasTPC((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
+  Double_t nsigmaTOF=pidResp->NumberOfSigmasTOF((AliVParticle*)track,(AliPID::EParticleType)fParticlePid);
+  Bool_t isTPCPid=false, isTOFPid=false;
+
+  // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
+  if(fuseTPCPid && okTPC) {
+    for(Int_t j=0; j<fnPTPCBins; j++) {
+      //      cout<<" checking bin: ("<<fPTPCBinLimits[j]<<","<<fPTPCBinLimits[j+1]<<") should be nsigma < "<<fnSigmaTPC[j]<<endl;
+      if ((pPart>fPTPCBinLimits[j]) && (pPart<fPTPCBinLimits[j+1]) && nsigmaTPC<fnSigmaTPC[j]) isTPCPid=true;
+    }
+    if(pPart>fPmaxTPC) isTPCPid=true;
+  }
+
+  // If use TPC and TPC infos are ok, check whether the sigma is ok in the given p range
+  if(fuseTOFPid && okTOF) {
+    for(Int_t j=0; j<fnPTOFBins; j++) {
+      if ((pPart>fPTOFBinLimits[j]) && (pPart<fPTOFBinLimits[j+1]) && nsigmaTOF<fnSigmaTOF[j]) isTPCPid=true;
+    }
+    if(pPart>fPmaxTOF) isTOFPid=true;
+  }
+
+  isSelected = (isTPCPid || isTOFPid) ? true : false;
+
+  return isSelected;
+}
diff --git a/PWGPP/EvTrkSelection/AliSingleTrackEffCuts.h b/PWGPP/EvTrkSelection/AliSingleTrackEffCuts.h
new file mode 100644 (file)
index 0000000..82d4d85
--- /dev/null
@@ -0,0 +1,165 @@
+#ifndef ALISINGLETRACKEFFICUTS_H
+#define ALISINGLETRACKEFFICUTS_H
+
+#include <TString.h>
+#include "TObject.h"
+
+#include "AliVEvent.h"
+#include "AliMCEvent.h"
+#include "AliAnalysisCuts.h"
+#include "AliPID.h"
+#include "AliAODTrack.h"
+#include "AliSingleTrackEffCuts.h"
+
+
+class AliSingleTrackEffCuts : public AliAnalysisCuts 
+{
+ public:
+
+  AliSingleTrackEffCuts();
+  AliSingleTrackEffCuts(const char* name, const char* title);
+
+  AliSingleTrackEffCuts(const AliSingleTrackEffCuts& source);
+  AliSingleTrackEffCuts& operator=(const AliSingleTrackEffCuts& source);
+
+  virtual ~AliSingleTrackEffCuts();
+
+  virtual Bool_t IsSelected(TList* /* list */) { return kFALSE; }
+
+  //
+  // Event and track selections
+  //
+  // Event selection at MC level (z-vtx)
+  Bool_t IsMCEventSelected(TObject *obj);
+  // Event selection at reconstructed level (trigger, zvtx)
+  Bool_t IsRecoEventSelected(TObject *obj);
+  //
+  // Check generated particles (primary, charged, pdgcode)
+  Bool_t IsMCParticleGenerated(TObject *obj);
+  // Check generated particles (eta, pt)
+  Bool_t IsMCParticleInKineAcceptance(TObject *obj);
+  // Check if particle has left enough hits in the detectors (only at ESD level)
+  Bool_t IsMCParticleInReconstructable(TObject *obj);
+  // Check if reconstructed particle is in the acceptance (eta, pt)
+  Bool_t IsRecoParticleKineAcceptance(TObject *obj);
+  // Check if reconstructed particle passes the PID selections
+  //  (only at AOD level for now) if TPC or TOF accept
+  Bool_t IsRecoParticlePID(TObject *obj);
+
+  //
+  // Setters
+  //
+  // Set eta range for acceptance cut (both MC and reco level)
+  void SetEtaRange(Float_t etamin, Float_t etamax){ fEtaMin=etamin; fEtaMax=etamax; }
+  // Set rapidity range for acceptance cut (both MC and reco level)
+  void SetYRange(Float_t ymin, Float_t ymax){ fYMin=ymin; fYMax=ymax; }
+  // Set pt range for acceptance cut (both MC and reco level)
+  void SetPtRange(Float_t ptmin, Float_t ptmax){ fPtMin=ptmin; fPtMax=ptmax; }
+  // Set PDG code to be checked (by default =0, no check)
+  void SetPdgCode(Int_t pdgCode){ fPdgCode = pdgCode; fIsPdgCode=kTRUE; }
+  Int_t GetPdgCode() { return fPdgCode; }
+  // Set if check for charged particle
+  void SetIsCharged(Bool_t charge){ fIsCharged=charge; }
+
+  // Set/get for AOD/ESD analysis
+  void SetIsAOD(Bool_t flag){ fisAOD = flag; }
+  Bool_t IsAOD(){ return fisAOD; }
+
+  // Set for vertex type (0: not cut; 1: SPDZ; 2: SPD3D; 3: Tracks)
+  void SetMinVtxType(Int_t type=3) { fMinVtxType=type; }
+  void SetUseEventsWithOnlySPDVertex(Bool_t flag=kTRUE){ 
+    if(flag) fMinVtxType=1;
+    else fMinVtxType=3;
+  }
+  // Set minimum vertex contributors
+  void SetMinVtxContr(Int_t contr=1) { fMinVtxContr=contr; }
+  // Set maximum z-vtx cut
+  void SetMaxVtxZ(Float_t z=1e6) { fMaxVtxZ=z; }
+  // Appky or not cut on difference SPD-TPC vtx (0: no cut, 1: |zvtx-SPD - zvtx-TPC|<0.5cm)
+  void SetCutOnZVertexSPD(Int_t cut) { fCutOnZVertexSPD=cut; }
+
+  // Select event trigger mask
+  void SetTriggerMask(ULong64_t mask=0) { fTriggerMask=mask; }
+  UInt_t GetTriggerMask(){ return fTriggerMask; }
+
+  // Set minimum number of ITS, TPC, TOF or MUON clusters
+  void SetNumberOfClusters(Int_t nITS, Int_t nTPC, Int_t nTOF, Int_t nMUON){
+    fnClusITS = nITS; fnClusTPC = nTPC; fnClusTOF = nTOF; fnClusMUON = nMUON;
+  }
+
+  // PID setters (flag, particle specie, detector, limits)
+  void SetUsePid(Bool_t flag=kTRUE) { fusePid=flag; }
+  void SetParticleSpecie(AliPID::EParticleType type=AliPID::kPion) { fParticlePid=type; }
+  void SetUseTPCPid(Bool_t flag=kTRUE) { fuseTPCPid=flag; }
+  void SetUseTOFPid(Bool_t flag=kTRUE) { fuseTOFPid=flag; }
+  // set given number of sigma cut per P bin
+  void SetTPCSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin);
+  void SetTOFSigmaPtBins(Int_t nPtBins, Float_t *pBinLimits, Float_t *sigmaBin);
+  // maximum momentum to use PID
+  void SetMaximumPTPC(Float_t p) { fPmaxTPC=p; }
+  void SetMaximumPTOF(Float_t p) { fPmaxTOF=p; }
+  
+
+  //
+  // Getters
+  //
+  Bool_t GetUsePid() const{ return fusePid; }
+  Int_t  GetParticleSpecie() const { return fParticlePid; }
+  Float_t *GetPTPCBinLimits() const { return fPTPCBinLimits; }
+  Int_t  GetNPTPCBins() const {return fnPTPCBins; }
+  Float_t *GetPTOFBinLimits() const { return fPTOFBinLimits; }
+  Int_t  GetNPTOFBins() const { return fnPTOFBins; }
+
+
+ protected:
+
+  Bool_t IsVertexSelected(AliVEvent *event);
+  Bool_t CheckTPCPIDStatus(AliAODTrack *track) const;
+  Bool_t CheckTOFPIDStatus(AliAODTrack *track) const;
+
+  Bool_t fisAOD;  // flag wether it is AOD:1 or ESD:0 analysis
+
+  Bool_t fIsPdgCode; // flag to check pdg code
+  Int_t fPdgCode;    // particle pdg code
+
+  Float_t fEtaMin;   // minimum eta cut
+  Float_t fEtaMax;   // maximum eta cut
+  Float_t fYMin;     // minimum Y cut
+  Float_t fYMax;     // maximum Y cut
+  Float_t fPtMin;    // minimum Pt cut
+  Float_t fPtMax;    // maximum Pt cut
+  Bool_t  fIsCharged; // check if particle is charged (MC level)
+
+  UInt_t  fTriggerMask;   // event trigger mask
+  Int_t   fMinVtxType;    // 0: not cut; 1: SPDZ; 2: SPD3D; 3: Tracks
+  Int_t   fMinVtxContr;   // minimum vertex contributors
+  Float_t fMaxVtxZ;       // maximum |z| of primary vertex
+  Int_t fCutOnZVertexSPD; // 0: no cut, 1: |zvtx-SPD - zvtx-TPC|<0.5cm
+
+  Int_t fnClusITS;   // minimum number of ITS clusters
+  Int_t fnClusTPC;   // minimum number of TPC clusters
+  Int_t fnClusTOF;   // minimum number of TOF clusters
+  Int_t fnClusMUON;  // minimum number of MUON clusters
+
+  Bool_t    fusePid;          // flag to use or not Pid
+  Int_t     fParticlePid;     // integer to define the particle specie to check
+  //
+  Bool_t    fuseTPCPid;       // flag to use TPC Pid
+  Int_t     fnPTPCBins;       // "number of limits", that is fnPBins+1
+  Int_t     fnPTPCBinLimits;  // "number of limits", that is fnPBins+1
+  Float_t*  fPTPCBinLimits;   //[fnPTPCBinLimits]  p bins
+  Float_t*  fnSigmaTPC;       //[fnPTPCBins]
+  Float_t   fPmaxTPC;         // maximum TPC P to use Pid
+  //
+  Bool_t    fuseTOFPid;       // flag to use TOF Pid
+  Int_t     fnPTOFBins;       // "number of limits", that is fnPBins+1
+  Int_t     fnPTOFBinLimits;  // "number of limits", that is fnPBins+1
+  Float_t*  fPTOFBinLimits;   //[fnPTOFBinLimits]  p bins
+  Float_t*  fnSigmaTOF;       //[fnPTOFBins]
+  Float_t   fPmaxTOF;         // maximum TOF P to use Pid
+  
+
+  ClassDef(AliSingleTrackEffCuts,1)  // base class for cuts on AOD reconstructed heavy-flavour decays
+ };
+
+#endif
diff --git a/PWGPP/EvTrkSelection/macros/AddSingleTrackEfficiencyTask.C b/PWGPP/EvTrkSelection/macros/AddSingleTrackEfficiencyTask.C
new file mode 100644 (file)
index 0000000..b6e4a15
--- /dev/null
@@ -0,0 +1,291 @@
+//
+// Particle cuts
+//
+const Double_t etamin = -0.9;
+const Double_t etamax =  0.9;
+const Double_t ptmin = 0.0;
+const Double_t ptmax = 24.0;
+const Double_t phimin = -2*TMath::Pi();
+const Double_t phimax = 2*TMath::Pi();
+const Double_t thetamin = 0;
+const Double_t thetamax = TMath::Pi();
+const Double_t zvtxmin = -10.0;
+const Double_t zvtxmax =  10.0;
+//
+const Int_t mintrackrefsTPC = 5;
+const Int_t mintrackrefsITS = 4;
+const Int_t mintrackrefsTOF = 0;
+const Int_t mintrackrefsMUON = 0;
+const Int_t minclustersTPC = 70;
+const Int_t minclustersITS = 2;
+const Bool_t TPCRefit = kTRUE;
+const Bool_t ITSRefit = kFALSE;
+const Bool_t charge  = kTRUE;
+const Int_t  fBit = 0;
+
+//
+// Container settings
+//
+// Container mutliplicity bins
+const Float_t multmin_0_20 = 0;
+const Float_t multmax_0_20 = 20;
+const Float_t multmin_20_50 = 20;
+const Float_t multmax_20_50 = 50;
+const Float_t multmin_50_102 = 50;
+const Float_t multmax_50_102 = 102;
+//  Container Pt bins
+Double_t ptmin_0_2   = 0.0;
+Double_t ptmax_0_2   = 2.0;
+Double_t ptmin_2_6   = 2.0;
+Double_t ptmax_2_6   = 6.0;
+Double_t ptmin_6_8   = 6.0;
+Double_t ptmax_6_8   = 8.0;
+Double_t ptmin_8_16  = 8.0;
+Double_t ptmax_8_16  = 16.0;
+Double_t ptmin_16_24 = 16.0;
+Double_t ptmax_16_24 = 24.0;
+
+
+AliCFSingleTrackEfficiencyTask *AddSingleTrackEfficiencyTask(const Bool_t readAOD = 0, // Flag to read AOD:1 or ESD:0
+                                                            TString suffix="", // suffix for the output directory
+                                                            AliPID::EParticleType specie=AliPID::kPion, Int_t pdgcode=0 //particle specie
+                                                            )
+{
+
+  Info("AliCFSingleTrackEfficiencyTask","SETUP CONTAINER");
+
+  //
+  // Setting up the container
+  // 
+  // Variables
+  const Int_t nvar = 6; // number of variables on the grid: pt, y, phi, theta, zvtx, multiplicity
+  UInt_t nstep = 8;     // number of container steps
+  const UInt_t ipt = 0;
+  const UInt_t iy  = 1;
+  const UInt_t iphi = 2;
+  const UInt_t itheta = 3;
+  const UInt_t izvtx = 4;
+  const UInt_t imult = 5;
+  //
+  // Containter bining
+  //   A1. Bins variation by hand for pt
+  const Int_t nbinpt_0_2 = 8;  //bins in pt from 0 to 2 GeV
+  const Int_t nbinpt_2_6 = 8;   //bins in pt from 2 to 6 GeV
+  const Int_t nbinpt_6_8 = 2;   //bins in pt from 6 to 8 GeV
+  const Int_t nbinpt_8_16 = 4;  //bins in pt from 8 to 16 GeV
+  const Int_t nbinpt_16_24 = 1; //bins in pt from 16 to 24 GeV
+  //   A2. Bins variation by hand for other variables
+  const Int_t nbin2 = 9; //bins in eta
+  const Int_t nbin3 = 9; //bins in phi
+  const Int_t nbin4 = 9; //bins in theta
+  const Int_t nbin5 = 10; //bins in zvtx
+  //   A3. Bins for multiplicity
+  const Int_t nbinmult = 48;  //bins in multiplicity (total number)    
+  const Int_t nbinmult_0_20 = 20; //bins in multiplicity between 0 and 20
+  const Int_t nbinmult_20_50 = 15; //bins in multiplicity between 20 and 50
+  const Int_t nbinmult_50_102 = 13; //bins in multiplicity between 50 and 102
+
+  //arrays for the number of bins in each dimension
+  Int_t iBin[nvar];
+  iBin[0]=nbinpt_0_2+nbinpt_2_6+nbinpt_6_8+nbinpt_8_16+nbinpt_16_24;
+  iBin[1]=nbin2;
+  iBin[2]=nbin3;
+  iBin[3]=nbin4;
+  iBin[4]=nbin5;
+  iBin[5]=nbinmult;
+
+  //arrays for lower bounds :
+  Double_t *binLimpT = new Double_t[iBin[0]+1];
+  Double_t *binLim2 = new Double_t[iBin[1]+1];
+  Double_t *binLim3 = new Double_t[iBin[2]+1];
+  Double_t *binLim4 = new Double_t[iBin[3]+1];
+  Double_t *binLim5 = new Double_t[iBin[4]+1];
+  Double_t *binLimmult = new Double_t[iBin[5]+1];
+
+  // set the pt bins
+  for(Int_t i=0; i<=nbinpt_0_2; i++) binLimpT[i]=(Double_t)ptmin_0_2 + (ptmax_0_2-ptmin_0_2)/nbinpt_0_2*(Double_t)i ;
+  for(Int_t i=0; i<=nbinpt_2_6; i++) binLimpT[i+nbinpt_0_2]=(Double_t)ptmin_2_6 + (ptmax_2_6-ptmin_2_6)/nbinpt_2_6*(Double_t)i ;
+  for(Int_t i=0; i<=nbinpt_6_8; i++) binLimpT[i+nbinpt_0_2+nbinpt_2_6]=(Double_t)ptmin_6_8 + (ptmax_6_8-ptmin_6_8)/nbinpt_6_8*(Double_t)i ;
+  for(Int_t i=0; i<=nbinpt_8_16; i++) binLimpT[i+nbinpt_0_2+nbinpt_2_6+nbinpt_6_8]=(Double_t)ptmin_8_16 + (ptmax_8_16-ptmin_8_16)/nbinpt_8_16*(Double_t)i ;
+  for(Int_t i=0; i<=nbinpt_16_24; i++) binLimpT[i+nbinpt_0_2+nbinpt_2_6+nbinpt_6_8+nbinpt_8_16]=(Double_t)ptmin_16_24 + (ptmax_16_24-ptmin_16_24)/nbinpt_16_24*(Double_t)i;
+
+  // Other Variables
+  for(Int_t i=0; i<=nbin2; i++) binLim2[i]=(Double_t)etamin + (etamax-etamin)/nbin2*(Double_t)i ;
+  for(Int_t i=0; i<=nbin3; i++) binLim3[i]=(Double_t)phimin + (phimax-phimin)/nbin3*(Double_t)i ;
+  for(Int_t i=0; i<=nbin4; i++) binLim4[i]=(Double_t)thetamin + (thetamax-thetamin)/nbin4*(Double_t)i ;
+  for(Int_t i=0; i<=nbin5; i++) binLim5[i]=(Double_t)zvtxmin + (zvtxmax-zvtxmin)/nbin5*(Double_t)i ;
+
+  // multiplicity bining..
+  for(Int_t i=0; i<=nbinmult_0_20; i++) binLimmult[i]=(Double_t)multmin_0_20 + (multmax_0_20-multmin_0_20)/nbinmult_0_20*(Double_t)i ;
+  for(Int_t i=0; i<=nbinmult_20_50; i++) binLimmult[i+nbinmult_0_20]=(Double_t)multmin_20_50 + (multmax_20_50-multmin_20_50)/nbinmult_20_50*(Double_t)i ;
+  for(Int_t i=0; i<=nbinmult_50_102; i++) binLimmult[i+nbinmult_0_20+nbinmult_20_50]=(Double_t)multmin_50_102 + (multmax_50_102-multmin_50_102)/nbinmult_50_102*(Double_t)i ;
+
+
+  // Container  
+  AliCFContainer* container = new AliCFContainer("container","container for tracks",nstep,nvar,iBin);
+  container -> SetBinLimits(ipt,binLimpT);    // pt
+  container -> SetBinLimits(iy,binLim2);      // eta
+  container -> SetBinLimits(iphi,binLim3);    // phi
+  container -> SetBinLimits(itheta,binLim4);  // theta
+  container -> SetBinLimits(izvtx,binLim5);   // Zvtx
+  container -> SetBinLimits(imult,binLimmult);// multiplicity
+
+  // Variable Titles
+  container -> SetVarTitle(ipt,"pt");
+  container -> SetVarTitle(iy, "y");
+  container -> SetVarTitle(iphi,"phi");
+  container -> SetVarTitle(itheta, "theata");
+  container -> SetVarTitle(izvtx, "Zvtx");
+  container -> SetVarTitle(imult, "Multiplicity");
+
+  // Step Titles
+  container -> SetStepTitle(0, " MC Particle with Generated Cuts");
+  container -> SetStepTitle(1, " MC Particle with Kine Acceptance Cuts");
+  container -> SetStepTitle(2, " MC Particle with Track Ref Acceptance Cuts");
+  container -> SetStepTitle(3, " Total Reconstructed  Particle ");
+  container -> SetStepTitle(4, " Reco Particle With Kine Acceptance Cuts");
+  container -> SetStepTitle(5, " Reco Particle to MC True pt particles ");
+  container -> SetStepTitle(6, " Reco Particle With Quality Cuts");
+  container -> SetStepTitle(7, " Reco PID With Quality Cuts");
+
+
+  // SET TLIST FOR QA HISTOS
+  TList* qaList = new TList();
+  TObjArray* emptyList = new TObjArray(0);
+
+  //CREATE THE INTERFACE TO CORRECTION FRAMEWORK USED IN THE TASK
+  printf("CREATE INTERFACE AND CUTS\n");
+  AliCFManager* man = new AliCFManager();
+
+  man->SetNStepEvent(2);
+  man->SetEventContainer(container);
+  man->SetEventCutsList(0,emptyList);//evtmcList);
+  man->SetEventCutsList(1,emptyList);//evtrecoList);
+  
+  man->SetParticleContainer(container);
+  man->SetParticleCutsList(0,emptyList);//mcGenList);
+  man->SetParticleCutsList(1,emptyList);//mcKineList);
+  man->SetParticleCutsList(2,emptyList);//mcaccList);
+  man->SetParticleCutsList(3,emptyList);//evtrecoPureList);
+  man->SetParticleCutsList(4,emptyList);//recKineList);
+  man->SetParticleCutsList(5,emptyList);
+  man->SetParticleCutsList(6,emptyList);
+  man->SetParticleCutsList(7,emptyList);
+  
+  // Simulated particle & event cuts
+  AliSingleTrackEffCuts* cuts = new AliSingleTrackEffCuts();
+  cuts->SetPtRange(ptmin,ptmax);
+  cuts->SetEtaRange(etamin,etamax);
+  cuts->SetIsCharged(charge);
+  cuts->SetMinVtxContr(1);
+  cuts->SetMaxVtxZ(zvtxmax);
+  cuts->SetNumberOfClusters(mintrackrefsITS,mintrackrefsTPC,mintrackrefsTOF,mintrackrefsMUON);
+  cuts->SetTriggerMask(AliVEvent::kMB);
+  cuts->SetIsAOD(readAOD);
+  //
+  // Pid selection here
+  //
+  if(pdgcode>0){
+    cuts->SetUsePid(true);
+    cuts->SetParticleSpecie(specie);
+    cuts->SetPdgCode(pdgcode);
+    // 
+    const Int_t nlims=1;
+    Float_t plims[nlims+1]={0.,999.}; //TPC limits in momentum [GeV/c]
+    Float_t sigmas[nlims]={3.};
+    cuts->SetUseTPCPid();
+    cuts->SetTPCSigmaPtBins(nlims,plims,sigmas);
+    cuts->SetMaximumPTPC(4.);
+    // 
+    const Int_t nlims2=1;
+    Float_t plims2[nlims2+1]={0.,999.}; //TPC limits in momentum [GeV/c]
+    Float_t sigmas2[nlims2]={3.};
+    cuts->SetUseTOFPid();
+    cuts->SetTOFSigmaPtBins(nlims2,plims2,sigmas2);
+    cuts->SetMaximumPTOF(4.);
+  }
+
+  //
+  //  Track Quality cuts via ESD track cuts
+  //
+  AliESDtrackCuts* QualityCuts = new AliESDtrackCuts();
+  QualityCuts->SetRequireSigmaToVertex(kFALSE);
+  QualityCuts->SetMinNClustersTPC(minclustersTPC);
+  QualityCuts->SetMinNClustersITS(minclustersITS);
+  QualityCuts->SetRequireTPCRefit(TPCRefit);
+  QualityCuts->SetRequireITSRefit(ITSRefit);
+  QualityCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
+  QualityCuts->SetMinDCAToVertexXY(0.);
+  QualityCuts->SetEtaRange(etamin,etamax);
+  QualityCuts->SetPtRange(ptmin,ptmax);
+
+
+  //CREATE THE TASK
+  printf("CREATE CF Single track task\n");
+
+  AliCFSingleTrackEfficiencyTask *task = new AliCFSingleTrackEfficiencyTask("AliCFSingleTrackEfficiencyTask",QualityCuts,cuts);
+  task->SetFilterBit(kTRUE);
+  task->SetFilterType(fBit);
+  task->SelectCollisionCandidates(AliVEvent::kMB);
+  task->SetCFManager(man); //here is set the CF manager
+
+  //
+  // Get the pointer to the existing analysis manager via the static access method.
+  //==============================================================================
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTask", "No analysis manager to connect to.");
+    return NULL;
+  }
+
+  // This task requires an ESD or AOD input handler and an AOD output handler.
+  // Check this using the analysis manager.
+  //===============================================================================
+  TString type = mgr->GetInputEventHandler()->GetDataType();
+  if (!type.Contains("ESD") && !type.Contains("AOD")) {
+    ::Error("AddSingleTrackEfficiencyTask", "AliCFSingleTrackEfficiency task needs the manager to have an ESD or AOD input handler.");
+    return NULL;
+  }
+  
+  mgr->AddTask(task);
+  printf(" Create the output container\n");
+
+  //
+  // Create and connect containers for input/output
+  //
+  // ----- output data -----
+  TString outputfile = AliAnalysisManager::GetCommonFileName();
+  TString input1name="cchain0";
+  TString output2name="HistEventsProcessed", output3name="container",output4name="list",output5name="ESDtrackCuts",output6name="MCtrackCuts";
+  outputfile += ":PWGPP_CFSingleTrack";
+  outputfile += suffix;
+  output2name += suffix;
+  output3name += suffix;
+  output4name += suffix;
+  output5name += suffix;
+  output6name += suffix;
+
+
+  // ------ input data ------
+  AliAnalysisDataContainer *cinput0  = mgr->GetCommonInputContainer();
+  // ----- output data -----
+  // output TH1I for event counting
+  AliAnalysisDataContainer *coutput1 = mgr->CreateContainer(output2name, TH1I::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  // output Correction Framework Container (for acceptance & efficiency calculations)
+  AliAnalysisDataContainer *coutput2 = mgr->CreateContainer(output3name, AliCFContainer::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  // output QA histograms
+  AliAnalysisDataContainer *coutput3 = mgr->CreateContainer(output4name, TList::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  // output ESD track cuts for book keeping
+  AliAnalysisDataContainer *coutput4 = mgr->CreateContainer(output5name, AliESDtrackCuts::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+  // output event and particle selection cuts for book keeping
+  AliAnalysisDataContainer *coutput5 = mgr->CreateContainer(output6name, AliSingleTrackEffCuts::Class(),AliAnalysisManager::kOutputContainer,outputfile.Data());
+
+  mgr->ConnectInput(task,0,mgr->GetCommonInputContainer());
+  mgr->ConnectOutput(task,1,coutput1);
+  mgr->ConnectOutput(task,2,coutput2);
+  mgr->ConnectOutput(task,3,coutput3);
+  mgr->ConnectOutput(task,4,coutput4);
+  mgr->ConnectOutput(task,5,coutput5);
+
+  return task;
+}
diff --git a/PWGPP/EvTrkSelection/macros/RebinCFContainer.C b/PWGPP/EvTrkSelection/macros/RebinCFContainer.C
new file mode 100644 (file)
index 0000000..752e71f
--- /dev/null
@@ -0,0 +1,219 @@
+#include "TDirectoryFile.h"
+#include "THnBase.h"
+#include "THnSparse.h"
+#include "TH1.h"
+#include "TH1D.h"
+#include "TAxis.h"
+#include "TFile.h"
+#include "AliCFGridSparse.h"
+#include "AliCFContainer.h"
+#include "AliSingleTrackEffCuts.h"
+#include "AliCFSingleTrackEfficiencyTask.h"
+#include "AliCFEffGrid.h"
+
+#include "TROOT.h"
+#include "TStyle.h"
+#include "TSystem.h"
+#include "Riostream.h"
+
+//
+// Macro to project (and rebin) the output of CFSingleTrackEfficiencyTask
+//
+// Parameters:
+//   input file name
+//   variable to project and rebin
+//   steps to project (1:GenAcc/LimAcc, 2: RecCuts/GenAcc, 3: RecPID/GenAcc, 4: Rec/GenAcc, 5: RecAcc/GenAcc)
+//   container directory suffix
+//
+
+void RebinCFContainer(const char *infile="AnalysisResults.root",Int_t rebinVar=0, Int_t myEff=3, const char * name="Nch"){
+
+  gSystem->SetIncludePath("-I. -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/RAW -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/PWGHF/vertexingHF -I$ALICE_ROOT/PWGHF/vertexingHF/charmFlow -g");
+
+  gROOT->LoadMacro("AliSingleTrackEffCuts.cxx++g");
+  gROOT->LoadMacro("AliCFSingleTrackEfficiencyTask.cxx++g");
+
+  gROOT->SetStyle("Plain");
+  gStyle->SetPalette(1);
+  gStyle->SetOptStat(0);
+  gStyle->SetPalette(1);
+  gStyle->SetCanvasColor(0);
+  gStyle->SetFrameFillColor(0);
+  gStyle->SetOptTitle(0);
+  
+  
+  TFile* file = TFile::Open(infile,"read");
+  
+  TDirectoryFile *d = 0;
+  AliCFContainer *data = 0;
+  d = (TDirectoryFile*)file->Get(Form("PWGPP_CFSingleTrack%s",name));
+  if(!d) {
+    cout<<" no directory "<<endl;
+    return;
+  }
+  data = (AliCFContainer*) (d->Get("container"));
+  if(!data){
+    cout <<" no container"<<endl;
+  }
+
+  //
+  // *********** NUMERATOR
+  //
+  AliCFGridSparse* gridSparse = 0;             
+  if( myEff == 1 )  gridSparse = (AliCFGridSparse*)data->GetGrid(1); // GenAcc
+  if( myEff == 2 )  gridSparse = (AliCFGridSparse*)data->GetGrid(6); // RecCuts
+  if( myEff == 3 )  gridSparse = (AliCFGridSparse*)data->GetGrid(7); // RecPID
+  if( myEff == 4 )  gridSparse = (AliCFGridSparse*)data->GetGrid(3); // Rec (no cuts)
+  if( myEff == 5 )  gridSparse = (AliCFGridSparse*)data->GetGrid(4); // RecAcc
+
+  THnSparse* numData = (THnSparse*)gridSparse->GetGrid();
+  
+  // method 2: defining a new THnSparse changing only the axis of interest
+  Printf("Method 2 ");
+  THnSparse* newnumData = (THnSparse*)numData->Clone("numNew");
+  newnumData->Reset();
+
+  Int_t nLimits=0;
+  Double_t* newLimits =0;
+  if(rebinVar==0) {
+    nLimits = 18;
+    newLimits = new Double_t[nLimits+1];
+    newLimits[0]=0;
+    newLimits[1]=0.5;
+    newLimits[2]=1;
+    newLimits[3]=1.5;
+    newLimits[4]=2;
+    newLimits[5]=2.5;
+    newLimits[6]=3;
+    newLimits[7]=4;
+    newLimits[8]=5;
+    newLimits[9]=6;
+    newLimits[10]=7;
+    newLimits[11]=8;
+    newLimits[12]=9;
+    newLimits[13]=10;
+    newLimits[14]=11;
+    newLimits[15]=12;
+    newLimits[16]=14;
+    newLimits[17]=16;
+    newLimits[18]=24;
+  } else if (rebinVar==1) {
+    nLimits = 9;
+    newLimits = new Double_t[nLimits+1];
+    newLimits[0]=-0.9;
+    newLimits[1]=-0.7;
+    newLimits[2]=-0.5;
+    newLimits[3]=-0.3;
+    newLimits[4]=-0.1;
+    newLimits[5]=0.1;
+    newLimits[6]=0.3;
+    newLimits[7]=0.5;
+    newLimits[8]=0.7;
+    newLimits[9]=0.9;
+  }
+  const Int_t nnewLimits = nLimits;
+  
+  
+  TAxis* axis = (TAxis*)newnumData->GetAxis(rebinVar);
+  axis->Set(nnewLimits,newLimits);
+  newnumData->SetBinEdges(rebinVar,newLimits);
+  
+  newnumData->RebinnedAdd(numData, 1);
+  
+  // checking the bin contents
+  TH1D* h1 = (TH1D*)numData->Projection(rebinVar);
+  Float_t sum = 0;
+  Float_t sumnew = 0;
+  Printf("Original THnSparse");
+  for (Int_t ibin = 1; ibin<=h1->GetNbinsX(); ibin++){
+    Printf("ibin = %d, low edge = %f, content = %f",ibin,h1->GetBinLowEdge(ibin),h1->GetBinContent(ibin));
+    sum+=h1->GetBinContent(ibin);
+  }
+
+  Printf("THnSparse changing only one axis");
+  TH1D* h2num = (TH1D*)newnumData->Projection(rebinVar);
+  for (Int_t ibin = 1; ibin<=h2num->GetNbinsX(); ibin++){
+    Printf("ibin = %d, low edge = %f, content = %f",ibin,h2num->GetBinLowEdge(ibin),h2num->GetBinContent(ibin));
+    sumnew+=h2num->GetBinContent(ibin);
+  }
+  
+  Printf("sum = %f, sumnew = %f",sum, sumnew);
+
+  //  
+  // *********** DENOMINATOR RECPID
+  //
+  Printf("DENOMINATOR");
+  
+  AliCFGridSparse* gridSparse2 = 0;
+  if( myEff == 1 )  gridSparse2 = (AliCFGridSparse*)data->GetGrid(0); // LimAcc
+  else gridSparse2 = (AliCFGridSparse*)data->GetGrid(1);              // GenAcc
+
+  THnSparse* denData = (THnSparse*)gridSparse2->GetGrid();
+  
+  // method 2: defining a new THnSparse changing only the axis of interest
+  Printf("Method 2 ");
+  THnSparse* newdenData = (THnSparse*)denData->Clone("denNew");
+  newdenData->Reset();
+  
+  TAxis* axis2 = (TAxis*)newdenData->GetAxis(rebinVar);
+  axis2->Set(nnewLimits,newLimits);
+  newdenData->SetBinEdges(rebinVar,newLimits);
+  
+  newdenData->RebinnedAdd(denData, 1);
+  
+  // checking the bin contents
+  TH1D* h1d = (TH1D*)denData->Projection(rebinVar);
+  sum = 0;
+  sumnew = 0;
+  Printf("Original THnSparse");
+  for (Int_t ibin = 1; ibin<=h1d->GetNbinsX(); ibin++){
+    Printf("ibin = %d, low edge = %f, content = %f",ibin,h1d->GetBinLowEdge(ibin),h1d->GetBinContent(ibin));
+    sum+=h1d->GetBinContent(ibin);
+  }
+  
+  Printf("THnSparse changing only one axis");
+  TH1D* h2 = (TH1D*)newdenData->Projection(rebinVar);
+  for (Int_t ibin = 1; ibin<=h2->GetNbinsX(); ibin++){
+    Printf("ibin = %d, low edge = %f, content = %f",ibin,h2->GetBinLowEdge(ibin),h2->GetBinContent(ibin));
+    sumnew+=h2->GetBinContent(ibin);
+  }
+  
+  Printf("sum = %f, sumnew = %f",sum, sumnew);
+  
+  TH1D* heff = (TH1D*)h2num->Clone("heff");
+  heff->Divide(h2num, h2,1,1,"B");
+  //  heff->GetXaxis()->SetRangeUser(0,23.5);
+  TFile* fout ;
+  if( myEff == 1 ) {
+    fout = new TFile(Form("efficiency_GenAcc_over_LimAcc_rebinned_%d_%s.root",rebinVar,name),"RECREATE");
+    heff->Write("hpteffCFLimAccGenAcc");
+    h2num->Write("hptLimAcc");
+    h2->Write("hptGenAcc");
+  }
+  if( myEff == 2 ) {
+    fout = new TFile(Form("efficiency_RecPPR_over_GenAcc_rebinned_%d_%s.root",rebinVar,name),"RECREATE");
+    heff->Write("hpteffCFRecPPRGenAcc");
+    h2num->Write("hptRecPPR");
+    h2->Write("hptGenAcc");
+  }
+  if( myEff == 3 ) {
+    fout = new TFile(Form("efficiency_RecPID_over_GenAcc_rebinned_%d_%s.root",rebinVar,name),"RECREATE");
+    heff->Write("hpteffCFRecPIDGenAcc");
+    h2num->Write("hptRecPID");
+    h2->Write("hptGenAcc");
+  }
+  if( myEff == 4 ) {
+    fout = new TFile(Form("efficiency_Rec_over_GenAcc_rebinned_%d_%s.root",rebinVar,name),"RECREATE");
+    heff->Write("hpteffCFRecGenAcc");
+    h2num->Write("hptRec");
+    h2->Write("hptGenAcc");
+  }
+  if( myEff == 5 ) {
+    fout = new TFile(Form("efficiency_RecAcc_over_GenAcc_rebinned_%d_%s.root",rebinVar,name),"RECREATE");
+    heff->Write("hpteffCFRecAccGenAcc");
+    h2num->Write("hptRecAcc");
+    h2->Write("hptGenAcc");
+  }
+  fout->Close();
+
+}
diff --git a/PWGPP/EvTrkSelection/macros/RunCFSingleTrackEfficiencyTask.C b/PWGPP/EvTrkSelection/macros/RunCFSingleTrackEfficiencyTask.C
new file mode 100644 (file)
index 0000000..f663a8d
--- /dev/null
@@ -0,0 +1,155 @@
+class AliAnalysisGrid;
+class AliAnalysisAlien;
+
+//_______________________________| Loading Libraries |________________________________
+void Load() {
+    
+  gSystem->SetIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS  -I$ALICE_ROOT/OADB -g"); 
+
+  //load the required aliroot libraries
+  gSystem->Load("libTree.so");
+  gSystem->Load("libGeom.so");
+  gSystem->Load("libPhysics.so");
+  gSystem->Load("libVMC.so");
+  gSystem->Load("libMinuit.so");
+  gSystem->Load("libSTEERBase.so");
+  gSystem->Load("libESD.so");
+  gSystem->Load("libAOD.so"); 
+  gSystem->Load("libANALYSIS.so");
+  gSystem->Load("libOADB.so");
+  gSystem->Load("libANALYSISalice.so");
+  gSystem->Load("libCORRFW.so");
+}
+
+//_______________________________| Running on Grid |________________________________
+void RunCFSingleTrackEfficiencyTask()
+{
+
+  const Bool_t readAOD = 1;
+
+  TString        analysisMode    =  "local"; // "local", "grid", or "proof"
+  TString           inputMode    =  "list"; // "list", "xml", or "dataset"
+  Long64_t           nentries    =   123567890,firstentry=0;
+  Bool_t       useAlienPlugin    =   kTRUE;
+  TString          pluginmode    =   "test";
+  TString testfileslistWithPlugin="filesAOD.txt"; // list of local files to test a la "files.txt" to use by the plugin
+
+  TBenchmark benchmark;
+  benchmark.Start("AliCFSingleTrackEfficiencyTask");
+
+  Load();
+
+  if(analysisMode=="grid") {
+    TGrid::Connect("alien://") ;    //  Create an AliRunTagCuts and an AliEventTagCuts Object
+  }
+
+  if(useAlienPlugin) {
+    AliAnalysisGrid *alienHandler = CreateAlienHandler(pluginmode,testfileslistWithPlugin,readAOD);
+    if(!alienHandler) return;
+  }
+
+  printf("CREATE ANALYSIS MANAGER\n");
+  AliAnalysisManager *mgr = new AliAnalysisManager("My Manager","My Manager");
+  mgr->SetDebugLevel(10);
+  if(useAlienPlugin) mgr->SetGridHandler(alienHandler);
+
+  AliMCEventHandler*  mcHandler = new AliMCEventHandler();
+  if (!readAOD) mgr->SetMCtruthEventHandler(mcHandler);
+        
+  AliInputEventHandler* dataHandler;
+  if   (readAOD) dataHandler = new AliAODInputHandler();
+  else           dataHandler = new AliESDInputHandler();
+  mgr->SetInputEventHandler(dataHandler);
+        
+  //CREATE THE TASK
+  printf("Prepare to create the task\n");
+
+  // Run physics selection if not reading AODs
+  if (!readAOD) {
+    gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPhysicsSelection.C");
+    AliPhysicsSelectionTask* physSelTask = AddTaskPhysicsSelection(kTRUE);
+  }
+
+  // Add new tasks
+  //
+  // First add the task for the PID response setting
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDResponse.C");
+  AliAnalysisTaskSE *setupTask = AddTaskPIDResponse(kTRUE,kTRUE);
+        
+  gROOT->LoadMacro("$ALICE_ROOT/ANALYSIS/macros/AddTaskPIDqa.C");
+  AliAnalysisTaskPIDqa *pidQA = AddTaskPIDqa();
+        
+  gROOT->LoadMacro("AliSingleTrackEffCuts.cxx++g");
+  gROOT->LoadMacro("AliCFSingleTrackEfficiencyTask.cxx++g");
+  gROOT->LoadMacro("AddSingleTrackEfficiencyTask.C");
+  AliCFSingleTrackEfficiencyTask *task = AddSingleTrackEfficiencyTask(readAOD,"Nch");
+  AliCFSingleTrackEfficiencyTask *taskPi = AddSingleTrackEfficiencyTask(readAOD,"Pion",AliPID::kPion,211);
+  AliCFSingleTrackEfficiencyTask *taskKa = AddSingleTrackEfficiencyTask(readAOD,"Kaon",AliPID::kKaon,321);
+  // Run the analysis
+  TChain * analysisChain=0;
+        
+  if(analysisChain) printf("CHAIN HAS %d ENTRIES\n",(Int_t)analysisChain->GetEntries());
+  if(!mgr->InitAnalysis()) return;
+  mgr->PrintStatus();
+  if(analysisMode=="grid" && !useAlienPlugin) analysisMode="local";
+  if(analysisMode!="proof") {
+    mgr->StartAnalysis(analysisMode.Data(),analysisChain,nentries,firstentry);
+  }
+        
+  benchmark.Stop("AliCFSingleTrackEfficiencyTask");
+  benchmark.Show("AliCFSingleTrackEfficiencyTask");
+        
+  return;  
+}
+
+
+//_______________________________| CreateAlienHandler |________________________________
+AliAnalysisGrid* CreateAlienHandler(TString pluginmode="test", TString testfileslistWithPlugin="", const Bool_t readAOD=kTRUE)
+{
+
+  AliAnalysisAlien *plugin = new AliAnalysisAlien();
+  plugin->SetRunMode(pluginmode.Data());
+  plugin->SetUser("zconesa");
+  plugin->SetAPIVersion("V1.1x");
+  plugin->SetROOTVersion("v5-34-08");
+  plugin->SetAliROOTVersion("v5-05-03-AN");
+  plugin->SetNtestFiles(1);
+
+  // Set data file list to test on local mode  
+  plugin->SetFileForTestMode(testfileslistWithPlugin.Data());
+
+  // Set data search pattern for DATA grid Mode
+  plugin->SetGridDataDir("/alice/sim/2013/LHC13d3"); // specify MC sample
+  if(readAOD) plugin->SetDataPattern("AOD/*AliAOD.root");// specify AOD set
+  else plugin->SetDataPattern("*/AliESDs.root");
+  
+  Int_t totruns=0;
+  gROOT->LoadMacro("$ALICE_ROOT/PWGHF/vertexingHF/AddGoodRuns.C");
+  totruns += AddGoodRuns(plugin,"LHC13b","LHC13b");
+  totruns += AddGoodRuns(plugin,"LHC13c","LHC13c");
+  plugin->SetNrunsPerMaster(totruns);
+  
+  // plugin->AddDataFile("/alice/cern.ch/user/z/zconesa/sim/LHC13d3/195483_195529.xml");
+
+  plugin->SetGridWorkingDir("sim/LHC13d3/ST290713");
+  plugin->SetGridOutputDir("out"); 
+  
+  plugin->SetExecutable("ST290713.sh");
+  plugin->SetAnalysisSource("AliSingleTrackEffCuts.cxx AliCFSingleTrackEfficiencyTask.cxx");
+  plugin->SetAdditionalLibs("AliSingleTrackEffCuts.h AliSingleTrackEffCuts.cxx AliCFSingleTrackEfficiencyTask.cxx AliCFSingleTrackEfficiencyTask.h");  
+  plugin->AddIncludePath("-I. -I$ROOTSYS/include -I$ALICE_ROOT -I$ALICE_ROOT/include -I$ALICE_ROOT/ITS -I$ALICE_ROOT/TPC -I$ALICE_ROOT/CONTAINERS -I$ALICE_ROOT/STEER/STEER -I$ALICE_ROOT/STEER/STEERBase -I$ALICE_ROOT/STEER/ESD -I$ALICE_ROOT/STEER/AOD -I$ALICE_ROOT/TRD -I$ALICE_ROOT/macros -I$ALICE_ROOT/ANALYSIS -I$ALICE_ROOT/OADB -g"); 
+  
+  plugin->SetDefaultOutputs(kTRUE);
+  // merging via jdl
+  plugin->SetMergeViaJDL(kTRUE);
+  plugin->SetOneStageMerging(kFALSE);
+  plugin->SetMaxMergeStages(2);
+  
+  plugin->SetSplitMaxInputFileNumber(5);
+  
+  plugin->SetAnalysisMacro("ST290713.C");
+  plugin->SetJDLName("TaskHFST290713.jdl");
+  
+  return plugin;
+}