Update of the HFE package
authorsma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Sep 2010 17:33:46 +0000 (17:33 +0000)
committersma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 16 Sep 2010 17:33:46 +0000 (17:33 +0000)
20 files changed:
PWG3/hfe/AliAnalysisTaskDCA.cxx
PWG3/hfe/AliAnalysisTaskDisplacedElectrons.cxx
PWG3/hfe/AliAnalysisTaskHFE.cxx
PWG3/hfe/AliAnalysisTaskHFE.h
PWG3/hfe/AliHFEcontainer.cxx
PWG3/hfe/AliHFEcuts.cxx
PWG3/hfe/AliHFEcuts.h
PWG3/hfe/AliHFEextraCuts.cxx
PWG3/hfe/AliHFEmcQA.cxx
PWG3/hfe/AliHFEmcQA.h
PWG3/hfe/AliHFEpid.cxx
PWG3/hfe/AliHFEpid.h
PWG3/hfe/AliHFEpidQA.cxx
PWG3/hfe/AliHFEpidTOF.cxx
PWG3/hfe/AliHFEpidTPC.cxx
PWG3/hfe/AliHFEpidTRD.cxx
PWG3/hfe/AliHFEpidTRD.h
PWG3/hfe/AliHFEsecVtx.cxx
PWG3/hfe/AliHFEsecVtx.h
PWG3/hfe/AliHFEtools.cxx

index a61a340..afaf690 100644 (file)
@@ -91,7 +91,7 @@ AliAnalysisTaskDCA::AliAnalysisTaskDCA():
   DefineOutput(2, TList::Class());
 
   fDefaultPID = new AliESDpid;
-  fHFEpid = new AliHFEpid;
+  fHFEpid = new AliHFEpid("PIDforDCAanalysis");
 
 }
 
@@ -131,7 +131,7 @@ AliAnalysisTaskDCA::AliAnalysisTaskDCA(const char * name):
   DefineOutput(2, TList::Class());
   
   fDefaultPID = new AliESDpid;
-  fHFEpid = new AliHFEpid;
+  fHFEpid = new AliHFEpid("PIDforDCAanalysis");
 
 }
 
index 3132add..7b86e45 100644 (file)
@@ -89,7 +89,7 @@ AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons():
   // Initialize pid
   
   fDeDefaultPID = new AliESDpid;
-  fDePID = new AliHFEpid;
+  fDePID = new AliHFEpid("PIDforDisplacedElectronAnalysis");
   
 
 }
@@ -126,7 +126,7 @@ AliAnalysisTaskDisplacedElectrons::AliAnalysisTaskDisplacedElectrons(const char
 
   // Initialize pid
   fDeDefaultPID = new AliESDpid;
-  fDePID = new AliHFEpid;
+  fDePID = new AliHFEpid("PIDforDisplacedElectronAnalysis");
 
 }
 
index b59edb9..0fa0c70 100644 (file)
@@ -87,13 +87,19 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fWeighting(kFALSE)
   , fWeightFactors(NULL)
   , fWeightFactorsFunction(NULL)
+  , fBackGroundFactorsFunction(NULL)
   , fCFM(NULL)
+  , fV0CF(NULL)
   , fHadronicBackground(NULL)
   , fCorrelation(NULL)
   , fPIDperformance(NULL)
   , fSignalToBackgroundMC(NULL)
   , fPID(NULL)
+  , fPIDtagged(NULL)
+  , fPIDpreselect(NULL)
   , fCuts(NULL)
+  , fCutsTagged(NULL)
+  , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
@@ -121,13 +127,19 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fWeighting(kFALSE)
   , fWeightFactors(NULL)
   , fWeightFactorsFunction(NULL)
+  , fBackGroundFactorsFunction(NULL)
   , fCFM(NULL)
+  , fV0CF(NULL)
   , fHadronicBackground(NULL)
   , fCorrelation(NULL)
   , fPIDperformance(NULL)
   , fSignalToBackgroundMC(NULL)
   , fPID(NULL)
+  , fPIDtagged(NULL)
+  , fPIDpreselect(NULL)
   , fCuts(NULL)
+  , fCutsTagged(NULL)
+  , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
@@ -161,13 +173,19 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fWeighting(kFALSE)
   , fWeightFactors(NULL)
   , fWeightFactorsFunction(NULL)
+  , fBackGroundFactorsFunction(NULL)
   , fCFM(NULL)
+  , fV0CF(NULL)
   , fHadronicBackground(NULL)
   , fCorrelation(NULL)
   , fPIDperformance(NULL)
   , fSignalToBackgroundMC(NULL)
   , fPID(NULL)
+  , fPIDtagged(NULL)
+  , fPIDpreselect(NULL)
   , fCuts(NULL)
+  , fCutsTagged(NULL)
+  , fCutspreselect(NULL)
   , fSecVtx(NULL)
   , fElecBackGround(NULL)
   , fMCQA(NULL)
@@ -209,13 +227,18 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fWeighting = fWeighting;
   target.fWeightFactors = fWeightFactors;
   target.fWeightFactorsFunction = fWeightFactorsFunction;
+  target.fBackGroundFactorsFunction = fBackGroundFactorsFunction;
   target.fCFM = fCFM;
+  target.fV0CF = fV0CF;
   target.fHadronicBackground = fHadronicBackground;
   target.fCorrelation = fCorrelation;
   target.fPIDperformance = fPIDperformance;
   target.fSignalToBackgroundMC = fSignalToBackgroundMC;
   target.fPID = fPID;
+  target.fPIDtagged = fPIDtagged;
+  target.fPIDpreselect = fPIDpreselect;
   target.fCuts = fCuts;
+  target.fCutspreselect = fCutspreselect;
   target.fSecVtx = fSecVtx;
   target.fElecBackGround = fElecBackGround;
   target.fMCQA = fMCQA;
@@ -235,6 +258,7 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   //
   return;
   if(fPID) delete fPID;
+  if(fPIDtagged) delete fPIDtagged;
   if(fQA){
     fQA->Clear();
     delete fQA;
@@ -245,6 +269,7 @@ AliAnalysisTaskHFE::~AliAnalysisTaskHFE(){
   }
   if(fWeightFactors) delete fWeightFactors;
   if(fWeightFactorsFunction) delete fWeightFactorsFunction;
+  if(fBackGroundFactorsFunction) delete fBackGroundFactorsFunction;
   if(fHistMCQA){
     fHistMCQA->Clear();
     delete fHistMCQA;
@@ -283,7 +308,7 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   // QA histograms are created if requested
   // Called once per worker
   //
-  fPID = new AliHFEpid;
+  fPID = new AliHFEpid("standardPID"); fPIDtagged = new AliHFEpid("taggedPID");
   AliDebug(3, "Creating Output Objects");
   // Automatic determination of the analysis mode
   AliVEventHandler *inputHandler = dynamic_cast<AliVEventHandler *>(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
@@ -321,6 +346,7 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   if(!fOutput) fOutput = new TList;
   // Initialize correction Framework and Cuts
   fCFM = new AliCFManager;
+  fV0CF = new AliCFManager;
   MakeParticleContainer();
   MakeEventContainer();
   // Temporary fix: Initialize particle cuts with NULL
@@ -332,8 +358,14 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
     fCuts->CreateStandardCuts();
   }
   if(IsAODanalysis()) fCuts->SetAOD();
+  // Make clone for V0 tagging step
+  fCutsTagged = new AliHFEcuts(*fCuts);
+  fCutsTagged->SetName("hfeV0Cuts");
+  fCutsTagged->SetTitle("Cuts for tagged Particles");
   fCuts->Initialize(fCFM);
-  if(fCuts->IsInDebugMode()) fQA->Add(fCuts->GetQAhistograms());
+  fCutsTagged->Initialize(fV0CF);
+  if(fCuts->IsQAOn()) fQA->Add(fCuts->GetQAhistograms());
+  if(fCutsTagged->IsQAOn()) fQA->Add(fCutsTagged->GetQAhistograms());
  
   // add output objects to the List
   fOutput->AddAt(fCFM->GetParticleContainer(), 0);
@@ -343,7 +375,7 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   fOutput->AddAt(fSignalToBackgroundMC, 4);
   fOutput->AddAt(fNElectronTracksEvent, 5);
   fOutput->AddAt(fHadronicBackground, 6);
-
+  fOutput->AddAt(fV0CF, 7);
   // Initialize PID
   if(IsQAOn(kPIDqa)){
     AliInfo("PID QA switched on");
@@ -353,54 +385,28 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   }
   fPID->SetHasMCData(HasMCData());
   if(!fPIDdetectors.Length() && ! fPIDstrategy) AddPIDdetector("TPC");
-  if(fPIDstrategy)
+  if(fPIDstrategy){
     fPID->InitializePID(Form("Strategy%d", fPIDstrategy));
-  else
+    fPIDtagged->InitializePID(Form("Strategy%d", fPIDstrategy));
+  }
+  else{
     fPID->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed 
+    fPIDtagged->InitializePID(fPIDdetectors.Data());     // Only restrictions to TPC allowed 
+  }
 
   // mcQA----------------------------------
   if (HasMCData() && IsQAOn(kMCqa)) {
     AliInfo("MC QA on");
     if(!fMCQA) fMCQA = new AliHFEmcQA;
     if(!fHistMCQA) fHistMCQA = new TList();
-    fHistMCQA->SetName("MCqa");
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm 
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_");       // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm 
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_");    // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_");     // create histograms for charm 
-    fMCQA->CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_");    // create histograms for beauty
-    fMCQA->CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_");    // create histograms for beauty
-    TIter next(gDirectory->GetList());
-    TObject *obj;
-    int counter = 0;
-    TString objname;
-    while ((obj = next.Next())) {
-      objname = obj->GetName();
-      TObjArray *toks = objname.Tokenize("_");
-      if (toks->GetEntriesFast()){
-        TObjString *fpart = (TObjString *)(toks->UncheckedAt(0));
-        if ((fpart->String()).CompareTo("mcqa") == 0) fHistMCQA->AddAt(obj, counter++);
-      }
-    }
+    fMCQA->CreatDefaultHistograms(fHistMCQA);
     fQA->Add(fHistMCQA);
   } 
 
   // secvtx----------------------------------
   if (GetPlugin(kSecVtx)) {
     AliInfo("Secondary Vertex Analysis on");
-    fSecVtx = new AliHFEsecVtx;
+    if(!fSecVtx) fSecVtx = new AliHFEsecVtx;
     fSecVtx->SetHasMCData(HasMCData());
 
     if(!fHistSECVTX) fHistSECVTX = new TList();
@@ -463,9 +469,13 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
     if(workingPID){
       AliDebug(1, "Using ESD PID from the input handler");
       fPID->SetESDpid(workingPID);
+      fPIDtagged->SetESDpid(workingPID);
+      if(fPIDpreselect) fPIDpreselect->SetESDpid(workingPID);
     } else { 
       AliDebug(1, "Using default ESD PID");
       fPID->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
+      fPIDtagged->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData()));
+      if(fPIDpreselect) fPIDpreselect->SetESDpid(AliHFEtools::GetDefaultPID(HasMCData())); 
     }
     ProcessESD();
   }
@@ -538,18 +548,16 @@ void AliAnalysisTaskHFE::ProcessMC(){
   // In case MC QA is on also MC QA loop is done
   //
   AliDebug(3, "Processing MC Information");
-  Double_t nContrib = 0;
-  const AliVVertex *pVertex = fMCEvent->GetPrimaryVertex();
-  if(pVertex) nContrib = pVertex->GetNContributors();
+  Double_t eventstatus = 1.;
   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepGenerated, fMCEvent)) return;
-  fCFM->GetEventContainer()->Fill(&nContrib,AliHFEcuts::kEventStepGenerated);
+  fCFM->GetEventContainer()->Fill(&eventstatus,AliHFEcuts::kEventStepGenerated);
   Int_t nElectrons = 0;
   if(IsESDanalysis()){
     if (HasMCData() && IsQAOn(kMCqa)) {
       AliDebug(2, "Running MC QA");
 
       if(fMCEvent->Stack()){
-       fMCQA->SetMCEvent(fMCEvent);
+             fMCQA->SetMCEvent(fMCEvent);
         fMCQA->SetGenEventHeader(fMCEvent->GenEventHeader());
         fMCQA->Init();
 
@@ -607,9 +615,10 @@ void AliAnalysisTaskHFE::ProcessESD(){
   // Loop over Tracks, filter according cut steps defined in AliHFEcuts
   //
   AliDebug(3, "Processing ESD Event");
-  Double_t nContrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
+  Double_t eventstatus = 1.;
+  fCFM->GetEventContainer()->Fill(&eventstatus, AliHFEcuts::kEventStepRecNoCut);
   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fInputEvent)) return;
-  fCFM->GetEventContainer()->Fill(&nContrib, AliHFEcuts::kEventStepReconstructed);
+  fCFM->GetEventContainer()->Fill(&eventstatus, AliHFEcuts::kEventStepReconstructed);
   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
   if(!fESD){
     AliError("ESD Event required for ESD Analysis")
@@ -632,7 +641,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
     }
   }
 
-
+  Double_t nContrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
   Double_t container[10];
   memset(container, 0, sizeof(Double_t) * 10);
   // container for the output THnSparse
@@ -662,11 +671,25 @@ void AliAnalysisTaskHFE::ProcessESD(){
   //
   AliDebug(3, Form("Number of Tracks: %d", fESD->GetNumberOfTracks()));
   for(Int_t itrack = 0; itrack < fESD->GetNumberOfTracks(); itrack++){
-    
     track = fESD->GetTrack(itrack);
 
+    // fill counts of v0-identified particles
+    Int_t v0pid = -1;
+    if(track->TestBit(BIT(14))) v0pid = AliPID::kElectron;
+    else if(track->TestBit(BIT(15))) v0pid = AliPID::kPion;
+    else if(track->TestBit(BIT(16))) v0pid = AliPID::kProton;
+    if(v0pid > -1)
+      FilterTaggedTrack(track, v0pid);
     AliDebug(3, Form("Doing track %d, %p", itrack, track));
-          
+     
+    //////////////////////////////////////
+    // preselect
+    /////////////////////////////////////
+    if(fPIDpreselect && fCutspreselect) {
+      if(!PreSelectTrack(track)) continue;
+    }
+     
     container[0] = track->Pt();
     container[1] = track->Eta();
     container[2] = track->Phi();
@@ -711,7 +734,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
         };
       } else if(IsGammaElectron(track)) container[4] = container[9] = kGammaConv;
       AliDebug(3, Form("Signal Decision(%f/%f)", container[4], container[9]));
-    }
+    } 
     AliDebug(3, Form("Weight? %f", weight));
     if(signal) {
       alreadyseen = cont.Find(TMath::Abs(track->GetLabel()));
@@ -786,10 +809,17 @@ void AliAnalysisTaskHFE::ProcessESD(){
 
     // Fill Containers
     if(signal) {
-      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack,weight);
-      fCFM->GetParticleContainer()->Fill(&container[5], AliHFEcuts::kStepPID,weight);
+      // Apply weight for background contamination
+      Double_t weightBackGround = 1.0;
+      if(fBackGroundFactorsFunction) {
+       weightBackGround = weightBackGround - fBackGroundFactorsFunction->Eval(TMath::Abs(track->P()));
+       if(weightBackGround < 0.0) weightBackGround = 1.0;
+      }
+      //      
+      fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack,weight*weightBackGround);
+      fCFM->GetParticleContainer()->Fill(&container[5], AliHFEcuts::kStepPID,weight*weightBackGround);
       if(alreadyseen) {
-        fCFM->GetParticleContainer()->Fill(&container[5], (AliHFEcuts::kStepPID + (AliHFEcuts::kNcutStepsESDtrack)),weight);
+        fCFM->GetParticleContainer()->Fill(&container[5], (AliHFEcuts::kStepPID + (AliHFEcuts::kNcutStepsESDtrack)),weight*weightBackGround);
       }
       // dimensions 3&4&5 : pt,eta,phi (MC)
       ((THnSparseF *)fCorrelation->At(1))->Fill(container);
@@ -894,9 +924,10 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   // Function is still in development
   //
   AliDebug(3, "Processing AOD Event");
-  Double_t nContrib = fInputEvent->GetPrimaryVertex()->GetNContributors();
+  Double_t eventstatus = 1.;
+  fCFM->GetEventContainer()->Fill(&eventstatus,AliHFEcuts::kEventStepRecNoCut);
   if(!fCFM->CheckEventCuts(AliHFEcuts::kEventStepReconstructed, fInputEvent)) return;
-  fCFM->GetEventContainer()->Fill(&nContrib,AliHFEcuts::kEventStepReconstructed);
+  fCFM->GetEventContainer()->Fill(&eventstatus,AliHFEcuts::kEventStepReconstructed);
   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
   if(!fAOD){
     AliError("AOD Event required for AOD Analysis")
@@ -953,7 +984,13 @@ void AliAnalysisTaskHFE::ProcessAOD(){
     if(HasMCData()) hfetrack.fMCtrack = mctrack;
     //if(!fPID->IsSelected(&hfetrack)) continue;    // we will do PID here as soon as possible
     // Particle identified - Fill CF Container
-    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack);
+    // Apply weight for background contamination
+    Double_t weightBackGround = 1.0;
+    if(fBackGroundFactorsFunction) {
+      weightBackGround = weightBackGround - fBackGroundFactorsFunction->Eval(TMath::Abs(track->P()));
+      if(weightBackGround < 0.0) weightBackGround = 1.0;
+    }
+    fCFM->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID + 2*AliHFEcuts::kNcutStepsESDtrack, weightBackGround);
     nElectronCandidates++;    
     if(HasMCData()){
       // Track selected: distinguish between true and fake
@@ -1060,19 +1097,90 @@ Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
 }
 
 //____________________________________________________________
+void AliAnalysisTaskHFE::FilterTaggedTrack(AliESDtrack *track, Int_t species){
+  //
+  // Filter tracks tagged by V0 PID class
+  //
+  Int_t offset = AliHFEcuts::kStepRecKineITSTPC;
+  Double_t container[5] ={track->Pt(), track->Eta(), track->Phi(), track->Charge(), species};
+  fV0CF->GetParticleContainer()->Fill(container, 0); // Fill Container without filtering
+  Bool_t survived = kTRUE;
+  for(Int_t icut = AliHFEcuts::kStepRecKineITSTPC; icut < AliHFEcuts::kStepPID; icut++){
+    AliDebug(2, Form("Checking cut %d for species %s", icut, AliPID::ParticleName(species)));
+    if(!fV0CF->CheckParticleCuts(icut, track)){
+      survived = kFALSE;
+      break;
+    }
+    AliDebug(2, Form("Cut passed, filling container %d", icut - offset + 1));
+    fV0CF->GetParticleContainer()->Fill(container, icut - offset + 1);
+  }
+  if(survived){
+    // Apply PID
+    AliHFEpidObject hfetrack;
+    hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
+    hfetrack.fRecTrack = track;
+    if(fPIDtagged->IsSelected(&hfetrack)) fV0CF->GetParticleContainer()->Fill(container, AliHFEcuts::kStepPID - offset + 1); 
+  } 
+}
+//____________________________________________________________
+Bool_t AliAnalysisTaskHFE::PreSelectTrack(AliESDtrack *track) const {
+  //
+  // Preselect tracks
+  //
+  
+
+  Bool_t survived = kTRUE;
+  
+  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecKineITSTPC, track)) {
+    survived = kFALSE;
+    //printf("Did not pass AliHFEcuts::kStepRecKineITSTPC\n");
+  }
+  //else printf("Pass AliHFEcuts::kStepRecKineITSTPC\n");
+  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepRecPrim, track)) {
+    survived = kFALSE;
+    //printf("Did not pass AliHFEcuts::kStepRecPrim\n");
+  }
+  //else printf("Pass AliHFEcuts::kStepRecPrim\n");
+  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsITS, track)) {
+    survived = kFALSE;
+    //printf("Did not pass AliHFEcuts::kStepHFEcutsITS\n");
+  }
+  //else printf("Pass AliHFEcuts::kStepHFEcutsITS\n");
+  if(!fCutspreselect->CheckParticleCuts(AliHFEcuts::kStepHFEcutsTRD, track)) {
+    survived = kFALSE;
+    //printf("Did not pass AliHFEcuts::kStepHFEcutsTRD\n");
+  }
+  //else printf("Pass AliHFEcuts::kStepHFEcutsTRD\n");
+  
+  if(survived){
+    // Apply PID
+    AliHFEpidObject hfetrack;
+    hfetrack.fAnalysisType = AliHFEpidObject::kESDanalysis;
+    hfetrack.fRecTrack = track;
+    if(!fPIDpreselect->IsSelected(&hfetrack)) {
+      //printf("Did not pass AliHFEcuts::kPID\n");
+      survived = kFALSE;
+    }
+    //else printf("Pass AliHFEcuts::kPID\n");
+  }
+
+  return survived; 
+      
+}
+//____________________________________________________________
 void AliAnalysisTaskHFE::MakeEventContainer(){
   //
   // Create the event container for the correction framework and link it
   //
   const Int_t kNvar = 1;  // number of variables on the grid: number of tracks per event
-  const Double_t kNTrackBound[2] = {-0.5, 200.5};
-  const Int_t kNBins = 201;
+  const Double_t kStatBound[2] = {0.5, 1.5};
+  const Int_t kNBins = 1;
 
   AliCFContainer *evCont = new AliCFContainer("eventContainer", "Container for events", AliHFEcuts::kNcutStepsEvent, kNvar, &kNBins);
 
-  Double_t *trackBins = AliHFEtools::MakeLinearBinning(kNBins, kNTrackBound[0], kNTrackBound[1]);
-  evCont->SetBinLimits(0,trackBins);
-  delete[] trackBins;
+  Double_t *statBins = AliHFEtools::MakeLinearBinning(kNBins, kStatBound[0], kStatBound[1]);
+  evCont->SetBinLimits(0, statBins);
+  delete[] statBins;
 
   fCFM->SetEventContainer(evCont);
 }
@@ -1085,13 +1193,13 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
   //
   const Int_t kNvar   = 5;
   //number of variables on the grid:pt,eta, phi, charge
-  const Double_t kPtbound[2] = {0.1, 10.};
+  const Double_t kPtbound[2] = {0.1, 20.};
   const Double_t kEtabound[2] = {-0.8, 0.8};
   const Double_t kPhibound[2] = {0., 2. * TMath::Pi()};
 
   //arrays for the number of bins in each dimension
   Int_t iBin[kNvar];
-  iBin[0] = 40; // bins in pt
+  iBin[0] = 44; // bins in pt
   iBin[1] =  8; // bins in eta 
   iBin[2] = 18; // bins in phi
   iBin[3] =  2; // bins in charge
@@ -1163,6 +1271,16 @@ void AliAnalysisTaskHFE::MakeParticleContainer(){
     fPIDperformance->SetBinEdges(idim, binEdges2[idim]);
     fSignalToBackgroundMC->SetBinEdges(idim, binEdges2[idim]); 
   }
+
+  // create correction framework container for V0-tagged particles, new bin limits in 4th bin
+  iBin[4] = 5;
+  delete binEdges[4];
+  binEdges[4] = AliHFEtools::MakeLinearBinning(iBin[4], 0, iBin[4]); // Numeric precision
+  AliCFContainer *tagged = new AliCFContainer("taggedTrackContainer", "Correction Framework Container for tagged tracks", AliHFEcuts::kNcutStepsESDtrack, kNvar, iBin);
+  for(Int_t ivar = 0; ivar < kNvar; ivar++)
+    tagged->SetBinLimits(ivar, binEdges[ivar]);
+  fV0CF->SetParticleContainer(tagged);
+
   for(Int_t ivar = 0; ivar < kNvar; ivar++)
     delete binEdges[ivar];
   for(Int_t ivar = kNvar; ivar < nDim; ivar++)
@@ -1196,7 +1314,7 @@ void AliAnalysisTaskHFE::PrintStatus() const {
   printf("\n");
   printf("\tQA: \n");
   printf("\t\tPID: %s\n", IsQAOn(kPIDqa) ? "YES" :  "NO");
-  printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsInDebugMode()) ? "YES" : "NO");
+  printf("\t\tCUTS: %s\n", (fCuts != NULL && fCuts->IsQAOn()) ? "YES" : "NO");
   printf("\t\tMC: %s\n", IsQAOn(kMCqa) ? "YES" : "NO");
   printf("\n");
 }
index 84ef254..88b713e 100644 (file)
@@ -80,17 +80,20 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     Bool_t HasMCData() const { return TestBit(kHasMCdata); }
     Bool_t GetPlugin(Int_t plug) const { return TESTBIT(fPlugins, plug); };
     void SetHFECuts(AliHFEcuts * const cuts) { fCuts = cuts; };
+    void SetHFECutsPreselect(AliHFEcuts * const cuts) { fCutspreselect = cuts; };
     void SetHFEElecBackGround(AliHFEelecbackground * const elecBackGround) { fElecBackGround = elecBackGround; };
     void SetQAOn(Int_t qaLevel) { SETBIT(fQAlevel, qaLevel); };
     void SwitchOnPlugin(Int_t plug);
     void SetHasMCData(Bool_t hasMC = kTRUE) { SetBit(kHasMCdata, hasMC); };
     void SetPIDdetectors(Char_t * const detectors){ fPIDdetectors = detectors; }
     void SetPIDStrategy(UInt_t strategy) { fPIDstrategy = strategy; }
+    void SetPIDPreselect(AliHFEpid * const cuts) { fPIDpreselect = cuts; };
     void AddPIDdetector(TString detector);
     void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
     void SetESDAnalysis() { SetBit(kAODanalysis, kFALSE); };
     void SetWeightFactors(TH3D * const weightFactors);
     void SetWeightFactorsFunction(TF1 * const weightFactorsFunction);
+    void SetBackGroundFactorsFunction(TF1 * const backGroundFactorsFunction) { fBackGroundFactorsFunction = backGroundFactorsFunction; };
     void PrintStatus() const;
  
   private:
@@ -127,6 +130,8 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     void ProcessMC();
     void ProcessESD();
     void ProcessAOD();
+    void FilterTaggedTrack(AliESDtrack *track, Int_t species);
+    Bool_t PreSelectTrack(AliESDtrack *track) const;
     Bool_t ProcessMCtrack(AliVParticle *track);
     Bool_t ProcessCutStep(Int_t cutStep, AliVParticle *track, Double_t *container, Bool_t signal, Bool_t alreadyseen,Double_t weight);
     
@@ -138,13 +143,19 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     Bool_t  fWeighting;                   // Weighting or not for the efficiency maps
     TH3D *fWeightFactors;                 // Weight factors
     TF1  *fWeightFactorsFunction;         // Weight factors
+    TF1  *fBackGroundFactorsFunction;     // BackGround factors
     AliCFManager *fCFM;                   //! Correction Framework Manager
+    AliCFManager *fV0CF;                  //! Correction Framework Manager for V0-tagged tracks
     AliCFContainer * fHadronicBackground; //! Container for hadronic Background
     TList *fCorrelation;                  //! response matrix for unfolding  
     THnSparseF *fPIDperformance;          //! info on contamination and yield of electron spectra
     THnSparseF *fSignalToBackgroundMC;    //! Signal To Background Studies on pure MC information
     AliHFEpid *fPID;                      //! PID
+    AliHFEpid *fPIDtagged;                // PID oject for tagged tracks (identical clone without QA)
+    AliHFEpid *fPIDpreselect;             // PID oject for pre-selected tracks (without QA)
     AliHFEcuts *fCuts;                    // Cut Collection
+    AliHFEcuts *fCutsTagged;              // Cut Collection for tagged tracks
+    AliHFEcuts *fCutspreselect;           // Cut Collection for pre-selected tracks
     AliHFEsecVtx *fSecVtx;                //! Secondary Vertex Analysis
     AliHFEelecbackground *fElecBackGround;//! Background analysis
     AliHFEmcQA *fMCQA;                    //! MC QA
index febc5c5..4afd1d9 100644 (file)
@@ -249,7 +249,7 @@ AliCFContainer *AliHFEcontainer::MakeMergedCFContainer(const Char_t *name, const
   Int_t *dummyBinning = new Int_t[fNVars];
   for(UInt_t ibin = 0; ibin < fNVars; ibin++) dummyBinning[ibin] = 1;
   AliCFContainer *cmerged = new AliCFContainer(name, title, nStepMerged, fNVars, dummyBinning);
-  delete [] dummyBinning;
+  delete dummyBinning;
   // Fill container with content
   AliInfo("Filling new container");
   Int_t cstep = 0;
index ed1307b..96cc250 100644 (file)
@@ -72,6 +72,27 @@ ClassImp(AliHFEcuts)
 
 //__________________________________________________________________
 AliHFEcuts::AliHFEcuts():
+  TNamed(),
+  fRequirements(0),
+  fMinClustersTPC(0),
+  fMinTrackletsTRD(0),
+  fCutITSPixel(0),
+  fCheckITSLayerStatus(kTRUE),
+  fMaxChi2clusterTPC(0.),
+  fMinClusterRatioTPC(0.),
+  fSigmaToVtx(0.),
+  fHistQA(0x0),
+  fCutList(0x0),
+  fDebugLevel(0)
+{
+  //
+  // Dummy Constructor
+  //
+}
+
+//__________________________________________________________________
+AliHFEcuts::AliHFEcuts(const Char_t *name, const Char_t *title):
+  TNamed(name, title),
   fRequirements(0),
   fMinClustersTPC(0),
   fMinTrackletsTRD(0),
@@ -90,21 +111,19 @@ AliHFEcuts::AliHFEcuts():
   memset(fProdVtx, 0, sizeof(Double_t) * 4);
   memset(fDCAtoVtx, 0, sizeof(Double_t) * 2);
   memset(fPtRange, 0, sizeof(Double_t) * 2);
-  fCutList = new TObjArray();
-  fCutList->SetOwner();
 }
 
 //__________________________________________________________________
 AliHFEcuts::AliHFEcuts(const AliHFEcuts &c):
-  TObject(c),
+  TNamed(c),
   fRequirements(c.fRequirements),
-  fMinClustersTPC(c.fMinClustersTPC),
-  fMinTrackletsTRD(c.fMinTrackletsTRD),
-  fCutITSPixel(c.fCutITSPixel),
-  fCheckITSLayerStatus(c.fCheckITSLayerStatus),
-  fMaxChi2clusterTPC(c.fMaxChi2clusterTPC),
-  fMinClusterRatioTPC(c.fMinClusterRatioTPC),
-  fSigmaToVtx(c.fSigmaToVtx),
+  fMinClustersTPC(0),
+  fMinTrackletsTRD(0),
+  fCutITSPixel(0),
+  fCheckITSLayerStatus(0),
+  fMaxChi2clusterTPC(0),
+  fMinClusterRatioTPC(0),
+  fSigmaToVtx(0),
   fHistQA(0x0),
   fCutList(0x0),
   fDebugLevel(0)
@@ -112,11 +131,66 @@ AliHFEcuts::AliHFEcuts(const AliHFEcuts &c):
   //
   // Copy Constructor
   //
-  memcpy(fProdVtx, c.fProdVtx, sizeof(Double_t) * 4);
-  memcpy(fDCAtoVtx, c.fDCAtoVtx, sizeof(Double_t) * 2);
-  memcpy(fPtRange, c.fPtRange, sizeof(Double_t) *2);
-  fCutList = dynamic_cast<TObjArray *>(c.fCutList->Clone());
-  fCutList->SetOwner();
+  c.Copy(*this);
+}
+
+//__________________________________________________________________
+AliHFEcuts &AliHFEcuts::operator=(const AliHFEcuts &c){
+  //
+  // Make assignment
+  //
+  if(&c != this) c.Copy(*this);
+  return *this;
+}
+
+//__________________________________________________________________
+void AliHFEcuts::Copy(TObject &c) const {
+  //
+  // Performing copy
+  //
+  AliHFEcuts &target = dynamic_cast<AliHFEcuts &>(c);
+
+  target.fRequirements = fRequirements;
+  target.fMinClustersTPC = fMinClustersTPC;
+  target.fMinTrackletsTRD = fMinTrackletsTRD;
+  target.fCutITSPixel = fCutITSPixel;
+  target.fCheckITSLayerStatus = fCheckITSLayerStatus;
+  target.fMaxChi2clusterTPC = fMaxChi2clusterTPC;
+  target.fMinClusterRatioTPC = fMinClusterRatioTPC;
+  target.fSigmaToVtx = fSigmaToVtx;
+  target.fDebugLevel = 0;
+
+  memcpy(target.fProdVtx, fProdVtx, sizeof(Double_t) * 4);
+  memcpy(target.fDCAtoVtx, fDCAtoVtx, sizeof(Double_t) * 2);
+  memcpy(target.fPtRange, fPtRange, sizeof(Double_t) *2);
+
+  // Copy cut List
+  if(target.fCutList){
+    target.fCutList->Clear();
+    delete target.fCutList;
+  }
+  if(fCutList){
+    target.fCutList = dynamic_cast<TObjArray *>(fCutList->Clone());
+    target.fCutList->SetOwner();
+  }
+  if(target.fHistQA){
+    target.fHistQA->Clear();
+    delete target.fHistQA;
+    target.fHistQA = NULL;
+  }
+  if(fHistQA){
+    // If the QA list was already produced, then we create it new, loop over the cuts and connect all the histos with this list
+    target.fHistQA = new TList;
+    target.fHistQA->SetName(Form("%s_CutQAhistograms", GetName()));
+    fHistQA->SetOwner(kFALSE);
+    TIter cit(target.fCutList);
+    TObjArray *clist = NULL;
+    AliCFCutBase *co = NULL;
+    while((clist = dynamic_cast<TObjArray *>(cit()))){
+      TIter cit1(clist);
+      while((co = dynamic_cast<AliCFCutBase *>(cit1()))) co->SetQAOn(target.fHistQA);
+    }
+  }
 }
 
 //__________________________________________________________________
@@ -140,9 +214,15 @@ void AliHFEcuts::Initialize(AliCFManager *cfm){
   // Publishes the cuts to the correction framework manager
   //
   AliDebug(2, "Called");
-  if(IsInDebugMode()){
-     fHistQA = new TList;
-    fHistQA->SetName("CutQAhistograms");
+  if(fCutList)
+    fCutList->Delete();
+  else{
+    fCutList = new TObjArray;
+    fCutList->SetOwner();
+  }
+  if(IsQAOn()){
+    fHistQA = new TList;
+    fHistQA->SetName(Form("%s_CutQAhistograms", GetName()));
     fHistQA->SetOwner(kFALSE);
   }
  
@@ -184,6 +264,17 @@ void AliHFEcuts::Initialize(AliCFManager *cfm){
 void AliHFEcuts::Initialize(){
   // Call all the setters for the cuts
   AliDebug(2, "Called\n");
+   if(fCutList)
+    fCutList->Delete();
+  else{
+    fCutList = new TObjArray;
+    fCutList->SetOwner();
+  }
+  if(IsQAOn()){
+    fHistQA = new TList;
+    fHistQA->SetName(Form("%s_CutQAhistograms", GetName()));
+    fHistQA->SetOwner(kFALSE);
+  }
   SetParticleGenCutList();
   SetAcceptanceCutList();
   SetRecKineITSTPCCutList();
@@ -251,7 +342,7 @@ void AliHFEcuts::SetParticleGenCutList(){
     genCuts->SetProdVtxRange2D(kTRUE);  // Use ellipse
   }
   genCuts->SetRequirePdgCode(11, kTRUE);
-  if(IsInDebugMode()) genCuts->SetQAOn(fHistQA);
+  if(IsQAOn()) genCuts->SetQAOn(fHistQA);
 
   // Add
   mcCuts->AddLast(genCuts);
@@ -261,7 +352,7 @@ void AliHFEcuts::SetParticleGenCutList(){
     AliCFTrackKineCuts *kineMCcuts = new AliCFTrackKineCuts((Char_t *)"fCutsKineMC", (Char_t *)"MC Kine Cuts");
     kineMCcuts->SetPtRange(fPtRange[0], fPtRange[1]);
     kineMCcuts->SetEtaRange(-0.8, 0.8);
-    if(IsInDebugMode()) kineMCcuts->SetQAOn(fHistQA);
+    if(IsQAOn()) kineMCcuts->SetQAOn(fHistQA);
     mcCuts->AddLast(kineMCcuts);
   }
    
@@ -283,7 +374,7 @@ void AliHFEcuts::SetAcceptanceCutList(){
   accCuts->SetMinNHitITS(3);
   accCuts->SetMinNHitTPC(2);
   accCuts->SetMinNHitTRD(2*fMinTrackletsTRD);
-  if(IsInDebugMode()) accCuts->SetQAOn(fHistQA);
+  if(IsQAOn()) accCuts->SetQAOn(fHistQA);
   
   TObjArray *partAccCuts = new TObjArray();
   partAccCuts->SetName("fPartAccCuts");
@@ -329,7 +420,7 @@ void AliHFEcuts::SetRecKineITSTPCCutList(){
   kineCuts->SetPtRange(fPtRange[0], fPtRange[1]);
   kineCuts->SetEtaRange(-0.8, 0.8);
   
-  if(IsInDebugMode()){
+  if(IsQAOn()){
     trackQuality->SetQAOn(fHistQA);
     hfecuts->SetQAOn(fHistQA);
     kineCuts->SetQAOn(fHistQA);
@@ -364,7 +455,7 @@ void AliHFEcuts::SetRecPrimaryCutList(){
     primaryCut->SetMaxNSigmaToVertex(fSigmaToVtx);
   }
   primaryCut->SetAcceptKinkDaughters(kFALSE);
-  if(IsInDebugMode()) primaryCut->SetQAOn(fHistQA);
+  if(IsQAOn()) primaryCut->SetQAOn(fHistQA);
   
   TObjArray *primCuts = new TObjArray;
   primCuts->SetName("fPartPrimCuts");
@@ -384,7 +475,7 @@ void AliHFEcuts::SetHFElectronITSCuts(){
     hfecuts->SetCheckITSstatus(fCheckITSLayerStatus);
   }
   
-  if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA);
+  if(IsQAOn()) hfecuts->SetQAOn(fHistQA);
   hfecuts->SetDebugLevel(fDebugLevel);
 
   TObjArray *hfeCuts = new TObjArray;
@@ -401,7 +492,7 @@ void AliHFEcuts::SetHFElectronTRDCuts(){
   AliDebug(2, "Called\n");
   AliHFEextraCuts *hfecuts = new AliHFEextraCuts("fCutsHFElectronGroupTRD","Extra cuts from the HFE group");
   if(fMinTrackletsTRD > 0.) hfecuts->SetMinTrackletsTRD(fMinTrackletsTRD);
-  if(IsInDebugMode()) hfecuts->SetQAOn(fHistQA);
+  if(IsQAOn()) hfecuts->SetQAOn(fHistQA);
   hfecuts->SetDebugLevel(fDebugLevel);
   
   TObjArray *hfeCuts = new TObjArray;
@@ -411,21 +502,13 @@ void AliHFEcuts::SetHFElectronTRDCuts(){
 }
 
 //__________________________________________________________________
-void AliHFEcuts::SetDebugMode(){ 
-  //
-  // Switch on QA
-  //
-  SetBit(kDebugMode, kTRUE); 
-}
-
-//__________________________________________________________________
 Bool_t AliHFEcuts::CheckParticleCuts(CutStep_t step, TObject *o){
   //
   // Checks the cuts without using the correction framework manager
   // 
   AliDebug(2, "Called\n");
-  TString stepnames[kNcutStepsTrack] = {"fPartGenCuts", "fPartAccCuts", "fPartRecCuts", "fPartPrimCuts", "fPartHFECuts"};
-  TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(stepnames[step].Data()));
+  TString stepnames[kNcutStepsTrack] = {"fPartGenCuts","fPartSignal","fPartAccCuts","fPartRecNoCuts","fPartRecKineITSTPCCuts", "fPartPrimCuts", "fPartHFECutsITS","fPartHFECutsTRD","fPartHFEPid"};
+ TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(stepnames[step].Data()));
   if(!cuts) return kTRUE;
   TIterator *it = cuts->MakeIterator();
   AliCFCutBase *mycut;
index 1981424..b3ef435 100644 (file)
@@ -20,8 +20,8 @@
 #ifndef ALIHFECUTS_H
 #define ALIHFECUTS_H
 
-#ifndef ROOT_TObject
-#include <TObject.h>
+#ifndef ROOT_TNamed
+#include <TNamed.h>
 #endif
 
 #ifndef ALIHFEEXTRACUTS_H
@@ -35,7 +35,7 @@ class AliMCParticle;
 class TObjArray;
 class TList;
 
-class AliHFEcuts : public TObject{
+class AliHFEcuts : public TNamed{
   public:
     typedef enum{
       kStepMCGenerated = 0,
@@ -50,17 +50,20 @@ class AliHFEcuts : public TObject{
     } CutStep_t;
     typedef enum{
       kEventStepGenerated = 0,
-      kEventStepReconstructed = 1
+      kEventStepRecNoCut = 1,
+      kEventStepReconstructed = 2
     } EventCutStep_t;
     enum{
-      kNcutStepsEvent = 2,
+      kNcutStepsEvent = 3,
       kNcutStepsTrack = 9,
       kNcutStepsESDtrack = 6 
     };    // Additional constants
 
     AliHFEcuts();
+    AliHFEcuts(const Char_t *name, const Char_t *title);
     AliHFEcuts(const AliHFEcuts &c);
     AliHFEcuts &operator=(const AliHFEcuts &c);
+    void Copy(TObject &o) const;
     ~AliHFEcuts();
     
     void Initialize(AliCFManager *cfm);
@@ -70,11 +73,11 @@ class AliHFEcuts : public TObject{
   
     TList *GetQAhistograms() const { return fHistQA; }
     
-    void SetDebugMode();
+    void SetQAOn() {SetBit(kDebugMode, kTRUE); };
+    void UnsetQA() {SetBit(kDebugMode, kFALSE); };
+    Bool_t IsQAOn() const { return TestBit(kDebugMode); };
     void SetAOD() { SetBit(kAOD, kTRUE); }
     void SetESD() { SetBit(kAOD, kFALSE); }
-    void UnsetDebugMode() { SetBit(kDebugMode, kFALSE); }
-    Bool_t IsInDebugMode() const { return TestBit(kDebugMode); };
     Bool_t IsAOD() const { return TestBit(kAOD); }
     Bool_t IsESD() const { return !TestBit(kAOD); }
     
index aa3a6ad..41cd385 100644 (file)
@@ -256,7 +256,8 @@ void AliHFEextraCuts::FillQAhistosESD(AliESDtrack *track, UInt_t when){
   (dynamic_cast<TH1F *>(container->At(3)))->Fill(track->GetTRDntrackletsPID());
   UChar_t itsPixel = track->GetITSClusterMap();
   TH1 *pixelHist = dynamic_cast<TH1F *>(container->At(4));
-  Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
+  //Int_t firstEntry = pixelHist->GetXaxis()->GetFirst();
+  Double_t firstEntry = 0.5;
   if(!((itsPixel & BIT(0)) || (itsPixel & BIT(1))))
     pixelHist->Fill(firstEntry + 3);
   else{
@@ -310,26 +311,26 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){
   TH1 *histo1D = 0x0;
   TH2 *histo2D = 0x0;
   histos[0] = new TList();
-  histos[0]->SetName("BeforeCut");
+  histos[0]->SetName(Form("%s_BeforeCut",GetName()));
   histos[0]->SetOwner();
   histos[1] = new TList();
-  histos[1]->SetName("AfterCut");
+  histos[1]->SetName(Form("%s_AfterCut",GetName()));
   histos[1]->SetOwner();
   TString cutstr[2] = {"before", "after"};
   for(Int_t icond = 0; icond < 2; icond++){
-    histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamR%s", cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
+    histos[icond]->AddAt((histo1D = new TH1F(Form("%s_impactParamR%s",GetName(),cutstr[icond].Data()), "Radial Impact Parameter", 100, 0, 10)), 0);
     histo1D->GetXaxis()->SetTitle("Impact Parameter");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
-    histos[icond]->AddAt((histo1D = new TH1F(Form("impactParamZ%s", cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
+    histos[icond]->AddAt((histo1D = new TH1F(Form("%s_impactParamZ%s",GetName(),cutstr[icond].Data()), "Z Impact Parameter", 200, 0, 20)), 1);
     histo1D->GetXaxis()->SetTitle("Impact Parameter");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
-    histos[icond]->AddAt((histo1D = new TH1F(Form("tpcClr%s", cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
+    histos[icond]->AddAt((histo1D = new TH1F(Form("%s_tpcClr%s",GetName(),cutstr[icond].Data()), "Cluster Ratio TPC", 10, 0, 1)), 2);
     histo1D->GetXaxis()->SetTitle("Cluster Ratio TPC");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
-    histos[icond]->AddAt((histo1D = new TH1F(Form("trdTracklets%s", cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
+    histos[icond]->AddAt((histo1D = new TH1F(Form("%s_trdTracklets%s",GetName(),cutstr[icond].Data()), "Number of TRD tracklets", 7, 0, 7)), 3);
     histo1D->GetXaxis()->SetTitle("Number of TRD Tracklets");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
-    histos[icond]->AddAt((histo1D = new TH1F(Form("itsPixel%s", cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
+    histos[icond]->AddAt((histo1D = new TH1F(Form("%s_itsPixel%s",GetName(),cutstr[icond].Data()), "ITS Pixel Hits", 6, 0, 6)), 4);
     histo1D->GetXaxis()->SetTitle("ITS Pixel");
     histo1D->GetYaxis()->SetTitle("Number of Tracks");
     Int_t first = histo1D->GetXaxis()->GetFirst();
@@ -339,11 +340,11 @@ void AliHFEextraCuts::AddQAHistograms(TList *qaList){
   }
   fQAlist = new TList();
   fQAlist->SetOwner();
-  fQAlist->SetName("HFelectronExtraCuts");
+  fQAlist->SetName(Form("%s_HFelectronExtraCuts",GetName()));
   fQAlist->AddAt(histos[0], 0);
   fQAlist->AddAt(histos[1], 1);
   // Add cut correlation
-  fQAlist->AddAt((histo2D = new TH2F("cutcorrellation", "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
+  fQAlist->AddAt((histo2D = new TH2F(Form("%s_cutcorrelation",GetName()), "Cut Correlation", kNcuts, 0, kNcuts - 1, kNcuts, 0, kNcuts -1)), 2);
   TString labels[kNcuts] = {"MinImpactParamR", "MaxImpactParamR", "MinImpactParamZ", "MaxImpactParamZ", "ClusterRatioTPC", "MinTrackletsTRD", "ITSpixel"};
   Int_t firstx = histo2D->GetXaxis()->GetFirst(), firsty = histo2D->GetYaxis()->GetFirst();
   for(Int_t icut = 0; icut < kNcuts; icut++){
index a9ed852..c513037 100644 (file)
@@ -93,9 +93,45 @@ AliHFEmcQA::~AliHFEmcQA()
 //_______________________________________________________________________________________________
 void AliHFEmcQA::PostAnalyze() const
 {
+        //
+        // Post analysis
+        //
 }
 
 //__________________________________________
+void AliHFEmcQA::CreatDefaultHistograms(TList * const qaList)
+{      
+  //
+  // make default histograms
+  //
+  
+  if(!qaList) return;
+
+  fQAhistos = qaList;
+  fQAhistos->SetName("MCqa");
+
+  CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
+  CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,0,"mcqa_");              // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
+  CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,1,"mcqa_barrel_");       // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
+  CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,2,"mcqa_unitY_");        // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm 
+  CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,3,"mcqa_reccut_");       // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm 
+  CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,4,"mcqa_recpidcut_");    // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kCharm,5,"mcqa_secvtxcut_");     // create histograms for charm 
+  CreateHistograms(AliHFEmcQA::kBeauty,5,"mcqa_secvtxcut_");    // create histograms for beauty
+  CreateHistograms(AliHFEmcQA::kOthers,5,"mcqa_secvtxcut_");    // create histograms for beauty
+}
+  
+//__________________________________________
 void AliHFEmcQA::CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt) 
 {
   // create histograms
@@ -999,23 +1035,108 @@ Int_t AliHFEmcQA::GetSource(TParticle * const mcpart)
 }
 
 //__________________________________________
-void AliHFEmcQA::MakeHistograms(){
-  //
-  // Create the QA histograms
-  //
-  fQAhistos = new TList;
-  fQAhistos->SetName("MCqa");
+Int_t AliHFEmcQA::GetElecSource(TParticle * const mcpart)
+{
+  // decay particle's origin 
 
-  CreateHistograms(AliHFEmcQA::kCharm,0,"mcqa_");               // create histograms for charm
-  CreateHistograms(AliHFEmcQA::kBeauty,0,"mcqa_");              // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,1,"mcqa_barrel_");        // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,1,"mcqa_barrel_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,2,"mcqa_unitY_");         // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,2,"mcqa_unitY_");        // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,3,"mcqa_reccut_");        // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,3,"mcqa_reccut_");       // create histograms for beauty
-  CreateHistograms(AliHFEmcQA::kCharm,4,"mcqa_recpidcut_");     // create histograms for charm 
-  CreateHistograms(AliHFEmcQA::kBeauty,4,"mcqa_recpidcut_");    // create histograms for beauty
+  if ( abs(mcpart->GetPdgCode()) != AliHFEmcQA::kElectronPDG ) return -1;
+
+  Int_t origin = -1;
+  Bool_t isFinalOpenCharm = kFALSE;
+
+  if(!mcpart){
+    AliDebug(1, "no mcparticle, return\n");
+    return -1;
+  }
+
+  Int_t iLabel = mcpart->GetFirstMother();
+  if (iLabel<0){
+    AliDebug(1, "Stack label is negative, return\n");
+    return -1;
+  }
+
+  AliMCParticle *mctrack = NULL;
+  if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(iLabel))))) return -1; 
+  TParticle *partMother = mctrack->Particle();
+  Int_t maPdgcode = partMother->GetPdgCode();
+
+   // if the mother is charmed hadron  
+   if ( int(abs(maPdgcode)/100.) == kCharm || int(abs(maPdgcode)/1000.) == kCharm ) {
+
+     for (Int_t i=0; i<fNparents; i++){
+        if (abs(maPdgcode)==fParentSelect[0][i]){
+          isFinalOpenCharm = kTRUE;
+        }
+     }
+     if (!isFinalOpenCharm) return -1;
+
+     // iterate until you find B hadron as a mother or become top ancester 
+     for (Int_t i=1; i<fgkMaxIter; i++){
+
+        Int_t jLabel = partMother->GetFirstMother();
+        if (jLabel == -1){
+          origin = kDirectCharm;
+          return origin;
+        }
+        if (jLabel < 0){ // safety protection
+          AliDebug(1, "Stack label is negative, return\n");
+          return -1;
+        }
+
+        // if there is an ancester
+        if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(jLabel))))) return -1; 
+        TParticle* grandMa = mctrack->Particle();
+        Int_t grandMaPDG = grandMa->GetPdgCode();
+
+        for (Int_t j=0; j<fNparents; j++){
+           if (abs(grandMaPDG)==fParentSelect[1][j]){
+             origin = kBeautyCharm;
+             return origin;
+           }
+        }
+
+        partMother = grandMa;
+     } // end of iteration 
+   } // end of if
+   else if ( int(abs(maPdgcode)/100.) == kBeauty || int(abs(maPdgcode)/1000.) == kBeauty ) {
+     for (Int_t i=0; i<fNparents; i++){
+        if (abs(maPdgcode)==fParentSelect[1][i]){
+          origin = kDirectBeauty;
+          return origin;
+        }
+     }
+   } // end of if
+   else if ( abs(maPdgcode) == 22 ) {
+     origin = kGamma;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 111 ) {
+     origin = kPi0;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 221 ) {
+     origin = kEta;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 223 ) {
+     origin = kOmega;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 333 ) {
+     origin = kPhi;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 331 ) {
+     origin = kEtaPrime;
+     return origin;
+   } // end of if
+   else if ( abs(maPdgcode) == 113 ) {
+     origin = kRho0;
+     return origin;
+   } // end of if
+   else origin = kElse;
+
+   return origin;
 }
 
 //__________________________________________
index dc3fa5b..fe4e19c 100644 (file)
@@ -45,7 +45,7 @@ class AliHFEmcQA: public TObject {
   public: 
     enum heavyType {kCharm=4, kBeauty=5, kOthers=6, kElectronPDG=11};
     enum qType {kQuark, kantiQuark, kHadron, keHadron, kDeHadron, kElectron, kElectron2nd};
-    enum SourceType {kDirectCharm=1, kDirectBeauty=2, kBeautyCharm=3, kGamma=4, kPi0=5, kElse=6, kMisID=7};
+    enum SourceType {kDirectCharm=1, kDirectBeauty=2, kBeautyCharm=3, kGamma=4, kPi0=5, kElse=6, kMisID=7, kEta=8, kOmega=9, kPhi=10, kEtaPrime=11, kRho0=12};
     enum ProcessType {
       kPairCreationFromq,  kPairCreationFromg,  kFlavourExitation,  kGluonSplitting, kInitialPartonShower, kLightQuarkShower
     };
@@ -56,9 +56,9 @@ class AliHFEmcQA: public TObject {
 
     virtual ~AliHFEmcQA();
 
-    void MakeHistograms();
     TList *GetList() const { return fQAhistos; };
     void PostAnalyze() const;
+    void CreatDefaultHistograms(TList * const qaList); // create default histograms  
     void CreateHistograms(const Int_t kquark, Int_t icut, TString hnopt=""); // create histograms for mc qa analysis
     void SetMCEvent(AliMCEvent* const mcEvent){fMCEvent = mcEvent;} 
     void SetGenEventHeader(AliGenEventHeader* const mcHeader){fMCHeader=mcHeader;} // set stack pointer
@@ -70,7 +70,8 @@ class AliHFEmcQA: public TObject {
     void GetDecayedKine(TParticle *part, const Int_t kquark, const Int_t kdecayed, Int_t icut); // get decay electron kinematics distribution
                void GetDecayedKine(AliAODMCParticle *mcpart, const Int_t kquark, Int_t kdecayed, Int_t icut); // get decay electron kinematics for AOD 
     void EndOfEventAna(const Int_t kquark); // run analysis which should be done at the end of the event loop
-               Int_t GetSource(TParticle * const mcpart); // return electron source id 
+               Int_t GetSource(TParticle * const mcpart); // return source id 
+               Int_t GetElecSource(TParticle * const mcpart); // return electron source id 
                Int_t GetSource(AliAODMCParticle * const mcpart); // return electron source id for AOD
 
   protected:
index 41d4365..f50c2bc 100644 (file)
@@ -49,8 +49,9 @@ ClassImp(AliHFEpid)
 
 //____________________________________________________________
 AliHFEpid::AliHFEpid():
+  TNamed(),
   fEnabledDetectors(0),
-  fPIDstrategy(0),
+  fPIDstrategy(kUndefined),
   fQAlist(0x0),
   fDebugLevel(0),
   fCommonObjects(NULL)
@@ -62,11 +63,32 @@ AliHFEpid::AliHFEpid():
 }
 
 //____________________________________________________________
+AliHFEpid::AliHFEpid(const Char_t *name):
+  TNamed(name, ""),
+  fEnabledDetectors(0),
+  fPIDstrategy(kUndefined),
+  fQAlist(NULL),
+  fDebugLevel(0),
+  fCommonObjects(NULL)
+{
+  //
+  // Default constructor
+  // Create PID objects for all detectors
+  //
+  memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
+  fDetectorPID[kMCpid] = new AliHFEpidMC("MCPID");
+  fDetectorPID[kTPCpid] = new AliHFEpidTPC("TRDPID");
+  fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRDPID");
+  fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOFPID");
+
+}
+
+//____________________________________________________________
 AliHFEpid::AliHFEpid(const AliHFEpid &c):
-  TObject(c),
+  TNamed(c),
   fEnabledDetectors(c.fEnabledDetectors),
-  fPIDstrategy(0),
-  fQAlist(0x0),
+  fPIDstrategy(kUndefined),
+  fQAlist(NULL),
   fDebugLevel(c.fDebugLevel),
   fCommonObjects(NULL)
 {
@@ -74,20 +96,7 @@ AliHFEpid::AliHFEpid(const AliHFEpid &c):
   // Copy Constructor
   //
   memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
-  if(c.fDetectorPID[kMCpid])
-    fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
-  if(c.fDetectorPID[kTPCpid])
-    fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
-  if(c.fDetectorPID[kTRDpid])
-    fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
-  if(c.fDetectorPID[kTOFpid])
-    fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
-  if(c.IsQAOn()) SetQAOn();
-  if(c.HasMCData()) SetHasMCData(kTRUE);
-  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
-    if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
-    if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData(kTRUE);
-  }
+  c.Copy(*this);
 }
 
 //____________________________________________________________
@@ -95,31 +104,8 @@ AliHFEpid& AliHFEpid::operator=(const AliHFEpid &c){
   //
   // Assignment operator
   //
-  TObject::operator=(c);
-
-  if(this != &c){
-    fEnabledDetectors = c.fEnabledDetectors;
-    fPIDstrategy = c.fPIDstrategy;
-    fQAlist = 0x0;
-    fDebugLevel = c.fDebugLevel;
-  
-    memset(fDetectorPID, 0, sizeof(AliHFEpidBase *) * kNdetectorPID);
-    if(c.fDetectorPID[kMCpid])
-      fDetectorPID[kMCpid] = new AliHFEpidMC(*(dynamic_cast<AliHFEpidMC *>(c.fDetectorPID[kMCpid])));
-    if(c.fDetectorPID[kTPCpid])
-      fDetectorPID[kTPCpid] = new AliHFEpidTPC(*(dynamic_cast<AliHFEpidTPC *>(c.fDetectorPID[kTPCpid])));
-    if(c.fDetectorPID[kTRDpid])
-      fDetectorPID[kTRDpid] = new AliHFEpidTRD(*(dynamic_cast<AliHFEpidTRD *>(c.fDetectorPID[kTOFpid])));
-    if(c.fDetectorPID[kTOFpid])
-      fDetectorPID[kTOFpid] = new AliHFEpidTOF(*(dynamic_cast<AliHFEpidTOF *>(c.fDetectorPID[kTOFpid])));
-    if(c.IsQAOn()) SetQAOn();
-    if(c.HasMCData()) SetHasMCData(kTRUE);
-    for(Int_t idet = 0; idet < kNdetectorPID; idet++){
-      if(c.IsQAOn() && fDetectorPID[idet]) fDetectorPID[idet]->SetQAOn(fQAlist);
-      if(c.HasMCData() && fDetectorPID[idet]) fDetectorPID[idet]->SetHasMCData();
-    }
-  }
-  return *this; 
+  if(&c != this) c.Copy(*this);
+  return *this;
 }
 
 //____________________________________________________________
@@ -127,15 +113,42 @@ AliHFEpid::~AliHFEpid(){
   //
   // Destructor
   //
-  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
-    if(fDetectorPID[idet])
-      delete fDetectorPID[idet];
-  } 
-  if(fQAlist) delete fQAlist; fQAlist = 0x0;  // Each detector has to care about its Histograms
+  for(Int_t idet = 0; idet < kNdetectorPID; idet++)
+    if(fDetectorPID[idet]) delete fDetectorPID[idet];
+  if(fQAlist) delete fQAlist; fQAlist = NULL;  // Each detector has to care about its Histograms
   ClearCommonObjects();
 }
 
 //____________________________________________________________
+void AliHFEpid::Copy(TObject &o) const{
+  //
+  // Make copy
+  //
+  
+  TNamed::Copy(o);
+  AliHFEpid &target = dynamic_cast<AliHFEpid &>(o);
+  target.ClearCommonObjects();
+
+  target.fEnabledDetectors = fEnabledDetectors;
+  target.fPIDstrategy = fPIDstrategy;
+  if(target.fQAlist){
+    delete target.fQAlist; target.fQAlist = NULL;
+  }
+  if(fQAlist) target.fQAlist = new TList;
+  target.fQAlist = 0x0;
+  target.fDebugLevel = fDebugLevel;
+  // Copy detector PIDs
+  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
+    //Cleanup pointers in case of assignment
+    if(target.fDetectorPID[idet])  
+      delete target.fDetectorPID[idet];     
+    if(fDetectorPID[idet]) 
+      target.fDetectorPID[idet] = dynamic_cast<AliHFEpidBase *>(fDetectorPID[idet]->Clone());
+  }
+}
+
+//____________________________________________________________
 void AliHFEpid::AddCommonObject(TObject * const o){
   //
   // Add common object to the garbage collection
@@ -164,16 +177,15 @@ Bool_t AliHFEpid::InitializePID(TString arg){
   // + Initializes Detector PID objects
   // + Handles QA
   //
-  fDetectorPID[kMCpid] = new AliHFEpidMC("Monte Carlo PID"); // Always there
-  SETBIT(fEnabledDetectors, kMCpid);
-  // Initialize detector PIDs according to PID Strategies
+  
+  Bool_t initFail = kFALSE;
   if(arg.BeginsWith("Strategy")){
+    // Initialize detector PIDs according to PID Strategies
     arg.ReplaceAll("Strategy", "");
     fPIDstrategy = arg.Atoi();
-    AliInfo(Form("PID Strategy %d enabled", fPIDstrategy));
-    Int_t strategyStatus = kTRUE;
+    AliDebug(1, Form("%s - PID Strategy %d enabled", GetName(), fPIDstrategy));
     switch(fPIDstrategy){
-      case 0: break;    // Pure MC PID - only valid in MC mode
+      case 0: SwitchOnDetector(kMCpid); break;    // Pure MC PID - only valid in MC mode
       case 1: InitStrategy1(); break;
       case 2: InitStrategy2(); break;
       case 3: InitStrategy3(); break;
@@ -181,43 +193,53 @@ Bool_t AliHFEpid::InitializePID(TString arg){
       case 5: InitStrategy5(); break;
       case 6: InitStrategy6(); break;
       case 7: InitStrategy7(); break;
-      default: strategyStatus = kFALSE;
+      case 8: InitStrategy8(); break;
+      default: initFail = kFALSE;
     }
-    return strategyStatus;
-  }
-  // No Strategy defined, Initialize according to detectors specified
-  AliInfo(Form("Doing InitializePID for Detectors %s end", arg.Data()));
-  fDetectorPID[kITSpid] = new AliHFEpidITS("ITS development PID");  // Development version of the ITS pid, for the moment always there
-  SETBIT(fEnabledDetectors, kITSpid);
+  } else {
+    // No Strategy defined, Initialize according to detectors specified
+    AliDebug(1, Form("%s - Doing InitializePID for Detectors %s end", GetName(), arg.Data()));
   
-  TObjArray *detsEnabled = arg.Tokenize(":");
-  TIterator *detIterator = detsEnabled->MakeIterator();
-  TObjString *det = 0x0;
-  while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
-    if(det->String().CompareTo("TPC") == 0){
-      AliInfo("Doing TPC PID");
-      fDetectorPID[kTPCpid] = new AliHFEpidTPC("TPC PID");
-      SETBIT(fEnabledDetectors, kTPCpid);
-    } else if(det->String().CompareTo("TRD") == 0){
-      fDetectorPID[kTRDpid] = new AliHFEpidTRD("TRD PID");
-      SETBIT(fEnabledDetectors, kTRDpid);
-    } else if(det->String().CompareTo("TOF") == 0){
-      AliInfo("Doing TOF PID");
-      fDetectorPID[kTOFpid] = new AliHFEpidTOF("TOF PID");
-      SETBIT(fEnabledDetectors, kTOFpid);
+    TObjArray *detsEnabled = arg.Tokenize(":");
+    TIterator *detIterator = detsEnabled->MakeIterator();
+    TObjString *det = NULL;
+    Int_t detector = -1;
+    TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF"};
+    Int_t nDetectors = 0;
+    while((det = dynamic_cast<TObjString *>(detIterator->Next()))){
+      TString &detstring = det->String();
+      detector = -1;
+      for(Int_t idet = 0; idet < kNdetectorPID; idet++){
+        if(!detstring.CompareTo(detectors[idet])){
+          detector = idet;
+          break;
+        }
+      }
+      if(detector > -1){
+        SwitchOnDetector(detector);
+        nDetectors++;
+      } else AliError(Form("Detector %s not implemented (yet)", detstring.Data()));
     }
-    // Here is still space for ESD PID
+    if(!nDetectors) initFail = kTRUE;
   }
+  if(initFail){
+    AliError("Initializaion of the PID Failed");
+    return kFALSE;
+  }
+
   // Initialize PID Objects
   Bool_t status = kTRUE;
   for(Int_t idet = 0; idet < kNdetectorPID; idet++){
+    if(!IsDetectorOn(idet)) continue;
     if(fDetectorPID[idet]){ 
       status &= fDetectorPID[idet]->InitializePID();
       if(IsQAOn() && status) fDetectorPID[idet]->SetQAOn(fQAlist);
       if(HasMCData() && status) fDetectorPID[idet]->SetHasMCData();
     }
   }
+  PrintStatus();
   return status;
+  AliDebug(1, Form("%s - Done", GetName()));
 }
 
 //____________________________________________________________
@@ -230,7 +252,8 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
     // MC Event
     return (TMath::Abs(fDetectorPID[kMCpid]->IsSelected(track)) == 11);
   }
-  if(fPIDstrategy < 8){
+  if(fPIDstrategy < 9){
+    AliDebug(1, Form("%s - PID Strategy %d", GetName(), fPIDstrategy));
     Int_t pid = 0;
     switch(fPIDstrategy){
       case 0: pid = IdentifyStrategy0(track); break;
@@ -241,6 +264,7 @@ Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track){
       case 5: pid = IdentifyStrategy5(track); break;
       case 6: pid = IdentifyStrategy6(track); break;
       case 7: pid = IdentifyStrategy7(track); break;
+      case 8: pid = IdentifyStrategy8(track); break;
       default: break;
     }
     return pid;
@@ -273,10 +297,6 @@ Bool_t AliHFEpid::MakePidTpcTof(AliHFEpidObject *track){
   // Combines TPC and TOF PID decision
   //
   if(track->fAnalysisType != AliHFEpidObject::kESDanalysis) return kFALSE;
-  AliESDtrack *esdTrack = dynamic_cast<AliESDtrack *>(track->fRecTrack);
-  if(!esdTrack) return kFALSE;
-  // return if TOF information not available
-  if(!(esdTrack->GetStatus() & AliESDtrack::kTOFpid)) return kFALSE;
 
   AliHFEpidTOF *tofPID = dynamic_cast<AliHFEpidTOF*>(fDetectorPID[kTOFpid]);
   if(!tofPID){
@@ -448,12 +468,9 @@ void AliHFEpid::InitStrategy1(){
   //
   // TPC alone, 3-sigma cut
   //
-  AliHFEpidTPC *pid = new AliHFEpidTPC("strat1TPCpid");
+  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   pid->SetTPCnSigma(1);
-  Bool_t status = pid->InitializePID();
-  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
-  if(HasMCData() && status) pid->SetHasMCData();
-  fDetectorPID[kTPCpid] = pid;
+  SwitchOnDetector(kTPCpid);
 }
 
 //____________________________________________________________
@@ -461,13 +478,10 @@ void AliHFEpid::InitStrategy2(){
   //
   // TPC alone, symmetric 3 sigma cut and asymmetric sigma cut in the momentum region between 2GeV/c and 10 GeV/c and sigma between -1 and 100
   //
-  AliHFEpidTPC *pid = new AliHFEpidTPC("strat2TPCpid");
+  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   pid->SetTPCnSigma(3);
   pid->SetAsymmetricTPCsigmaCut(2., 10., 0., 4.);
-  Bool_t status = pid->InitializePID();
-  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
-  if(HasMCData() && status) pid->SetHasMCData();
-  fDetectorPID[kTPCpid] = pid;
+  SwitchOnDetector(kTPCpid);
 }
 
 //____________________________________________________________
@@ -475,13 +489,10 @@ void AliHFEpid::InitStrategy3(){
   //
   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
   //   
-  AliHFEpidTPC *pid = new AliHFEpidTPC("strat3TPCpid");
+  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   pid->SetTPCnSigma(3);
   pid->SetRejectParticle(AliPID::kPion, 0., -100., 10., 1.);
-  Bool_t status = pid->InitializePID();
-  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
-  if(HasMCData() && status) pid->SetHasMCData();
-  fDetectorPID[kTPCpid] = pid;
+  SwitchOnDetector(kTPCpid);
 }
 
 //____________________________________________________________
@@ -490,11 +501,9 @@ void AliHFEpid::InitStrategy4(){
   // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
   //
   InitStrategy1();
-  AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat4TRDpid");
-  Bool_t status = trdpid->InitializePID();
-  if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
-  if(HasMCData() && status) trdpid->SetHasMCData();
-  fDetectorPID[kTRDpid] = trdpid;
+  AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fDetectorPID[kTRDpid]);
+  trdpid->SetPIDMethod(AliHFEpidTRD::kLQ);
+  SwitchOnDetector(kTRDpid);
 }
 
 //____________________________________________________________
@@ -503,11 +512,7 @@ void AliHFEpid::InitStrategy5(){
   // TPC and TRD combined, TPC 3 sigma cut and TRD NN 90% el efficiency level above 2 GeV/c
   //
   InitStrategy1();
-  AliHFEpidTRD *trdpid = new AliHFEpidTRD("strat5TRDpid");
-  Bool_t status = trdpid->InitializePID();
-  if(IsQAOn() && status) trdpid->SetQAOn(fQAlist);
-  if(HasMCData() && status) trdpid->SetHasMCData();
-  fDetectorPID[kTRDpid] = trdpid;
+  SwitchOnDetector(kTRDpid);
 }
 
 //____________________________________________________________
@@ -515,39 +520,24 @@ void AliHFEpid::InitStrategy6(){
   //
   // Combined TPC-TOF PID, combination is discribed in the funtion MakePidTpcTof
   //
-  AliHFEpidTPC *tpcpid = new AliHFEpidTPC("strat6TPCpid");
-  AliHFEpidTOF *tofpid = new AliHFEpidTOF("strat6TOFpid");
+  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
+  AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
   tpcpid->SetTPCnSigma(2);
   tofpid->SetTOFnSigma(3);
-  Bool_t status = tpcpid->InitializePID();
   //TF1 *upperCut = new TF1("upperCut", "[0] * TMath::Exp([1]*x)", 0, 20);
   TF1 *upperCut = new TF1("upperCut", "[0]", 0, 20); // Use constant upper cut
   TF1 *lowerCut = new TF1("lowerCut", "[0] * TMath::Exp([1]*x)", 0, 20);
-  upperCut->SetParameter(0, 3.);
+  upperCut->SetParameter(0, 5.);
   //upperCut->SetParameter(0, 2.7);
   //upperCut->SetParameter(1, -0.4357);
-  lowerCut->SetParameter(0, -2.7);
-  lowerCut->SetParameter(1, -0.4357);
+  lowerCut->SetParameter(0, -2.65);
+  lowerCut->SetParameter(1, -0.6757);
   tpcpid->SetUpperSigmaCut(upperCut);
   tpcpid->SetLowerSigmaCut(lowerCut);
   AddCommonObject(upperCut);
   AddCommonObject(lowerCut);
-  if(!status)
-    AliError("Initialization of TPC PID failed");
-  Bool_t status1 = tofpid->InitializePID();
-  if(!status1)
-    AliError("Initialization of TOF PID failed");
-  status &= status1;
-  if(IsQAOn() && status){
-    tpcpid->SetQAOn(fQAlist);
-    tofpid->SetQAOn(fQAlist);
-  }
-  if(HasMCData() && status){
-    tpcpid->SetHasMCData();
-    tofpid->SetHasMCData();
-  }
-  fDetectorPID[kTPCpid] = tpcpid;
-  fDetectorPID[kTOFpid] = tofpid;
+  SwitchOnDetector(kTPCpid);
+  SwitchOnDetector(kTOFpid);
 }
 
 //____________________________________________________________
@@ -555,14 +545,28 @@ void AliHFEpid::InitStrategy7(){
   //
   // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
   //   
-  AliHFEpidTPC *pid = new AliHFEpidTPC("strat7TPCpid");
+  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   pid->SetTPCnSigma(2);
   pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
   pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
-  Bool_t status = pid->InitializePID();
-  if(IsQAOn() && status) pid->SetQAOn(fQAlist);
-  if(HasMCData() && status) pid->SetHasMCData();
-  fDetectorPID[kTPCpid] = pid;
+  SwitchOnDetector(kTPCpid);
+}
+
+//____________________________________________________________
+void AliHFEpid::InitStrategy8(){
+  //
+  // TOF, TRD and TPC together
+  // 
+  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
+  AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
+  AliHFEpidTRD *trdpid = dynamic_cast<AliHFEpidTRD *>(fDetectorPID[kTRDpid]);
+
+  tpcpid->SetTPCnSigma(3);
+  tofpid->SetTOFnSigma(3);
+  trdpid->SetPIDMethod(AliHFEpidTRD::kLQ);
+  SwitchOnDetector(kTPCpid);
+  SwitchOnDetector(kTOFpid);
+  SwitchOnDetector(kTRDpid);
 }
 
 
@@ -613,3 +617,33 @@ Bool_t AliHFEpid::IdentifyStrategy7(AliHFEpidObject *track){
   return kTRUE;
 }
 
+//____________________________________________________________
+Bool_t AliHFEpid::IdentifyStrategy8(AliHFEpidObject *track){
+  // 
+  // Identify TPC, TRD, TOF
+  //
+  if(TMath::Abs(fDetectorPID[kTOFpid]->IsSelected(track)) != 11) return kFALSE;
+  Int_t trdpid = TMath::Abs(fDetectorPID[kTRDpid]->IsSelected(track));
+  return (trdpid == 0 || trdpid == 11) && (TMath::Abs(fDetectorPID[kTPCpid]->IsSelected(track)) == 11);
+}
+
+//____________________________________________________________
+void AliHFEpid::PrintStatus() const {
+  //
+  // Print the PID configuration
+  //
+  printf("\n%s: Printing configuration\n", GetName());
+  printf("===============================================\n");
+  printf("PID Strategy: %d\n", fPIDstrategy);
+  printf("PID Detectors: \n");
+  Int_t npid = 0;
+  TString detectors[kNdetectorPID] = {"MC", "ESD", "ITS", "TPC", "TRD", "TOF"};
+  for(Int_t idet = 0; idet < kNdetectorPID; idet++){
+    if(IsDetectorOn(idet)){
+      printf("\t%s\n", detectors[idet].Data());
+      npid++;
+    }
+  }
+  if(!npid) printf("\tNone\n");
+  printf("\n");
+}
index bce1e20..6dd0d8c 100644 (file)
 #ifndef ALIHFEPID_H
 #define ALIHFEPID_H
 
-#ifndef ROOT_TObject
-#include <TObject.h>
+#ifndef ROOT_TNamed
+#include <TNamed.h>
 #endif
 
 #ifndef ALIHFEPIDBASE_H
 #include "AliHFEpidBase.h"
 #endif
 
+#include <climits>
+
 class AliESDpid;
 class AliESDtrack;
 class AliHFEpidBase;
@@ -36,11 +38,25 @@ class AliMCParticle;
 
 class TList;
 
-class AliHFEpid : public TObject{
+class AliHFEpid : public TNamed{
  public:
+    enum{
+      kUndefined = UINT_MAX 
+    };
+    enum DETtype_t {
+      kMCpid = 0,
+      kESDpid = 1,
+      kITSpid = 2,
+      kTPCpid = 3,
+      kTRDpid = 4,
+      kTOFpid = 5,
+      kNdetectorPID = 6
+    };
     AliHFEpid();
+    AliHFEpid(const Char_t *name);
     AliHFEpid(const AliHFEpid &c);
     AliHFEpid &operator=(const AliHFEpid &c);
+    void Copy(TObject &o) const;
     ~AliHFEpid();
     
     Bool_t InitializePID(TString argument);
@@ -53,6 +69,8 @@ class AliHFEpid : public TObject{
     void SetQAOn();
     void SetHasMCData(Bool_t hasMCdata = kTRUE) { SetBit(kHasMCData, hasMCdata); };
     TList *GetQAhistograms() const { return fQAlist; };
+    AliHFEpidBase *GetDetPID(DETtype_t det) const { return det < kNdetectorPID ? fDetectorPID[det] : NULL; }
+    void PrintStatus() const;
 
   protected:
     Bool_t MakePidTpcTof(AliHFEpidObject *track);
@@ -67,6 +85,7 @@ class AliHFEpid : public TObject{
     void InitStrategy5();
     void InitStrategy6();
     void InitStrategy7();
+    void InitStrategy8();
     Bool_t IdentifyStrategy0(AliHFEpidObject *track);
     Bool_t IdentifyStrategy1(AliHFEpidObject *track);
     Bool_t IdentifyStrategy2(AliHFEpidObject *track);
@@ -75,21 +94,13 @@ class AliHFEpid : public TObject{
     Bool_t IdentifyStrategy5(AliHFEpidObject *track);
     Bool_t IdentifyStrategy6(AliHFEpidObject *track);
     Bool_t IdentifyStrategy7(AliHFEpidObject *track);
+    Bool_t IdentifyStrategy8(AliHFEpidObject *track);
   private:
     enum{
       kIsQAOn = BIT(14),
       kHasMCData = BIT(15)
     };
     enum{
-      kMCpid = 0,
-      kESDpid = 1,
-      kITSpid = 2,
-      kTPCpid = 3,
-      kTRDpid = 4,
-      kTOFpid = 5,
-      kNdetectorPID = 6
-    };
-    enum{
       kCombinedTPCTRD=0
     };
     enum{
@@ -99,6 +110,17 @@ class AliHFEpid : public TObject{
 
     void AddCommonObject(TObject * const o);
     void ClearCommonObjects();
+    //-----Switch on/off detectors in PID sequence------
+    void SwitchOnDetector(UInt_t det){ 
+      if(det < kNdetectorPID) SETBIT(fEnabledDetectors, det);
+    }
+    void SwitchOffDetector(UInt_t det){
+      if(det < kNdetectorPID) CLRBIT(fEnabledDetectors, det);
+    }
+    Bool_t IsDetectorOn(UInt_t det) const {
+      return det < kNdetectorPID ? TESTBIT(fEnabledDetectors, det): kFALSE;
+    }
+    //--------------------------------------------------
 
     AliHFEpidBase *fDetectorPID[kNdetectorPID];     //! Detector PID classes
     UInt_t fEnabledDetectors;                       //  Enabled Detectors
index 38bedf3..a109d0e 100644 (file)
@@ -238,10 +238,10 @@ void AliHFEpidQA::Init(){
   //
   // event based histograms
   //
-  Int_t n_T0[2] = {10000, 100};
-  Double_t min_T0[2] = {114500, -1000};
-  Double_t max_T0[2] = {124500, 1000};
-  fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, n_T0, min_T0, max_T0);
+  Int_t nT0[2] = {10000, 100};
+  Double_t minT0[2] = {114500, -1000};
+  Double_t maxT0[2] = {124500, 1000};
+  fOutput->CreateTHnSparse("hEvent_T0", "T0 as a function of run number; run number; T0 (ns)", 2, nT0, minT0, maxT0);
 
   //
   // test the tender V0 supply
@@ -517,7 +517,7 @@ void AliHFEpidQA::FillElectronLikelihoods(TObjArray * const particles, Int_t spe
   // pion and proton efficiency and the thresholds
   //
   Long_t status = 0;
-  const Char_t *detname[4] = {"ITS", "TPC", "TRD", "TOF"};
+  Char_t *detname[4] = {"ITS", "TPC", "TRD", "TOF"};
   Char_t specname[256];
 
   switch(species){
@@ -1019,14 +1019,19 @@ TObjArray * AliHFEpidQA::MakeCleanListElectrons(TObjArray *electrons) const {
   TIter candidates(electrons);
   AliESDEvent *esd; AliAODEvent *aod;
   AliHFEV0info *hfetrack;
+  // const Double_t kSigmaTight = 1.;
+  // const Double_t kSigmaLoose = 4.;
+  const Double_t kSigmaTight = 0.5;
+  const Double_t kSigmaLoose = 0.5;
+  const Double_t shift = 1.;
   if((esd = dynamic_cast<AliESDEvent *>(fEvent))){
     AliESDtrack *track = NULL, *partnerTrack = NULL;
     while((hfetrack = dynamic_cast<AliHFEV0info *>(candidates()))){
       track = dynamic_cast<AliESDtrack *>(hfetrack->GetTrack());
       partnerTrack = esd->GetTrack(hfetrack->GetPartnerID());
-      Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron));
-      Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron));
-      if((nSigmaTrack < 1 && nSigmaPartner < 4) || (nSigmaTrack < 4 && nSigmaPartner < 1))
+      Double_t nSigmaTrack = TMath::Abs(fESDpid->NumberOfSigmasTPC(track, AliPID::kElectron) - shift);
+      Double_t nSigmaPartner = TMath::Abs(fESDpid->NumberOfSigmasTPC(partnerTrack, AliPID::kElectron) - shift);
+      if((nSigmaTrack < kSigmaTight && nSigmaPartner < kSigmaLoose) || (nSigmaTrack < kSigmaLoose && nSigmaPartner < kSigmaTight))
         tracks->Add(track);
     }
   } else {
index a29990a..5004d1b 100644 (file)
@@ -97,6 +97,7 @@ void AliHFEpidTOF::Copy(TObject &ref) const {
 
   target.fPID = fPID;          
   target.fQAList = fQAList;
+  target.fNsigmaTOF = fNsigmaTOF;
 
   AliHFEpidBase::Copy(ref);
 }
@@ -326,6 +327,7 @@ Int_t AliHFEpidTOF::MakePIDesdV3(AliESDtrack *track, AliMCParticle * /*mctrack*/
     return 0;
   }
   if(!(track->GetStatus() & AliESDtrack::kTOFpid)) return 0;
+
   Double_t t0 = fESDpid->GetTOFResponse().GetTimeZero();
   Double_t p = track->GetOuterParam() ? track->GetOuterParam()->P() : track->P();
 
index 41cae94..d5397f8 100644 (file)
@@ -166,6 +166,7 @@ Int_t AliHFEpidTPC::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mctrack){
     AliError("No ESD PID object available");
     return kFALSE;
   }
+  AliDebug(1, "Doing TPC PID based on n-Sigma cut approach");
   Float_t nsigma = fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron);
   if(IsQAon()) FillTPChistograms(esdTrack, mctrack, kFALSE);
   // exclude crossing points:
index 663b5a9..fc63ad5 100644 (file)
 // 
 #include <TH1F.h>
 #include <TH2F.h>
+#include <THnSparse.h>
 #include <TList.h>
 #include <TString.h>
 
 #include "AliAODTrack.h"
 #include "AliAODMCParticle.h"
 #include "AliESDtrack.h"
+#include "AliESDpid.h"
 #include "AliLog.h"
 #include "AliMCParticle.h"
 #include "AliPID.h"
 
+#include "AliHFEcollection.h"
 #include "AliHFEpidTRD.h"
 
 ClassImp(AliHFEpidTRD)
@@ -42,6 +45,7 @@ const Double_t AliHFEpidTRD::fgkVerySmall = 1e-12;
 //___________________________________________________________________
 AliHFEpidTRD::AliHFEpidTRD() :
     AliHFEpidBase()
+  , fMinP(1.)
   , fPIDMethod(kNN)
   , fContainer(0x0)
 {
@@ -54,6 +58,7 @@ AliHFEpidTRD::AliHFEpidTRD() :
 //___________________________________________________________________
 AliHFEpidTRD::AliHFEpidTRD(const char* name) :
     AliHFEpidBase(name)
+  , fMinP(1.)
   , fPIDMethod(kNN)
   , fContainer(0x0)
 {
@@ -66,6 +71,7 @@ AliHFEpidTRD::AliHFEpidTRD(const char* name) :
 //___________________________________________________________________
 AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
     AliHFEpidBase("")
+  , fMinP(1.)
   , fPIDMethod(kLQ)
   , fContainer(0x0)
 {
@@ -94,9 +100,10 @@ void AliHFEpidTRD::Copy(TObject &ref) const {
   //
   AliHFEpidTRD &target = dynamic_cast<AliHFEpidTRD &>(ref);
 
+  target.fMinP = fMinP;
   target.fPIDMethod = fPIDMethod;
   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
-  if(fContainer) target.fContainer = dynamic_cast<TList *>(fContainer->Clone());
+  if(fContainer) target.fContainer = dynamic_cast<AliHFEcollection *>(fContainer->Clone());
   AliHFEpidBase::Copy(ref);
 }
 
@@ -105,10 +112,8 @@ AliHFEpidTRD::~AliHFEpidTRD(){
   //
   // Destructor
   //
-  if(fContainer){
-    fContainer->Clear();
+  if(fContainer)
     delete fContainer;
-  }
 }
 
 //______________________________________________________
@@ -117,7 +122,10 @@ Bool_t AliHFEpidTRD::InitializePID(){
   // InitializePID: Load TRD thresholds and create the electron efficiency axis
   // to navigate 
   //
-  InitParameters();
+  if(fPIDMethod == kLQ)
+    InitParameters1DLQ();
+  else
+    InitParameters();
   return kTRUE;
 }
 
@@ -156,16 +164,18 @@ Int_t AliHFEpidTRD::MakePIDesd(AliESDtrack *esdTrack, AliMCParticle * /*mcTrack*
   // Does PID decision for ESD tracks as discribed above
   //
   Double_t p = esdTrack->GetOuterParam() ? esdTrack->GetOuterParam()->P() : esdTrack->P();
-  if(p < 2.0) return 0;
+  if(IsQAon())
+    FillStandardQA(0, esdTrack);
+  if(p < fMinP) return 0;
 
   Double_t pidProbs[AliPID::kSPECIES];
   esdTrack->GetTRDpid(pidProbs);
-  if(IsQAon())FillHistogramsLikelihood(0, p, pidProbs[AliPID::kElectron]);
-  Double_t threshold = GetTRDthresholds(0.91, p);
+  Double_t threshold = GetTRDthresholds(0.71, p);
   AliDebug(1, Form("Threshold: %f\n", threshold));
-  if(IsQAon()) (dynamic_cast<TH2F *>(fContainer->At(kHistTRDthresholds)))->Fill(p, threshold);
+  if(IsQAon()) fContainer->Fill("fTRDthresholds", p, threshold);
   if(pidProbs[AliPID::kElectron] > threshold){
-    if(IsQAon()) FillHistogramsLikelihood(1, p, pidProbs[AliPID::kElectron]);
+    if(IsQAon()) 
+      FillStandardQA(1, esdTrack);
     return 11;
   }
   return 211;
@@ -188,6 +198,7 @@ void AliHFEpidTRD::InitParameters(){
   // Fill the Parameters into an array
   //
 
+  AliDebug(2, "Loading threshold Parameter");
   // Parameters for 6 Layers
   fThreshParams[0] = -0.001839; // 0.7 electron eff
   fThreshParams[1] = 0.000276;
@@ -216,6 +227,41 @@ void AliHFEpidTRD::InitParameters(){
 }
 
 //___________________________________________________________________
+void AliHFEpidTRD::InitParameters1DLQ(){
+  //
+  // Init Parameters for 1DLQ PID (M. Fasel, Sept. 6th, 2010)
+  //
+
+  // Parameters for 6 Layers
+  AliDebug(2, Form("Loading threshold parameter for Method 1DLQ"));
+  fThreshParams[0] = -0.02241; // 0.7 electron eff
+  fThreshParams[1] = 0.05043;
+  fThreshParams[2] = 0.7925; 
+  fThreshParams[3] = 2.625;
+  fThreshParams[4] = 0.07438; // 0.75 electron eff
+  fThreshParams[5] = 0.05158;
+  fThreshParams[6] = 2.864;
+  fThreshParams[7] = 4.356;
+  fThreshParams[8] = 0.1977; // 0.8 electron eff
+  fThreshParams[9] = 0.05956;
+  fThreshParams[10] = 2.853;
+  fThreshParams[11] = 3.713;
+  fThreshParams[12] = 0.5206; // 0.85 electron eff
+  fThreshParams[13] = 0.03077;
+  fThreshParams[14] = 2.966;
+  fThreshParams[15] = 4.07;
+  fThreshParams[16] = 0.8808; // 0.9 electron eff
+  fThreshParams[17] = 0.002092;
+  fThreshParams[18] = 1.17;
+  fThreshParams[19] = 4.506;
+  fThreshParams[20] = 1.; // 0.95 electron eff
+  fThreshParams[21] = 0.;
+  fThreshParams[22] = 0.;
+  fThreshParams[23] = 0.;
+
+}
+
+//___________________________________________________________________
 void AliHFEpidTRD::GetParameters(Double_t electronEff, Double_t *parameters){
   //
   // return parameter set for the given efficiency bin
@@ -235,7 +281,7 @@ Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
   // Calculate mean over the last 2/3 slices
   //
   const Int_t kNSlices = 48;
-  AliDebug(2, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
+  AliDebug(3, Form("Number of Tracklets: %d\n", track->GetTRDpidQuality()));
   Double_t trdSlices[kNSlices], tmp[kNSlices];
   Int_t indices[48];
   Int_t icnt = 0;
@@ -253,8 +299,8 @@ Double_t AliHFEpidTRD::GetTRDSignalV1(AliESDtrack *track, Int_t mcPID){
     trdSlices[ien] = tmp[indices[ien]];
   Double_t trdSignal = TMath::Mean(static_cast<Int_t>(static_cast<Double_t>(icnt) * 2./3.), trdSlices);
   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
-  AliDebug(1, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
-  if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV1(trdSignal, mom, mcPID);
+  AliDebug(3, Form("PID Meth. 1: p[%f], TRDSignal[%f]", mom, trdSignal));
+  if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 1);
   return trdSignal;
 }
 
@@ -274,10 +320,10 @@ Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
     for(Int_t islice = 0; islice < 8; islice++){
       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
       if(islice < 5){
-        AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
+        AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
         trdSlicesLowTime[cntLowTime++] = track->GetTRDslice(idet, islice);
       } else{
-        AliDebug(2, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
+        AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
         trdSlicesRemaining[cntRemaining++] = track->GetTRDslice(idet, islice);
       }
     }
@@ -288,29 +334,23 @@ Double_t AliHFEpidTRD::GetTRDSignalV2(AliESDtrack *track, Int_t mcPID){
     trdSlicesRemaining[cntRemaining++] = trdSlicesLowTime[indices[ien]];
   Double_t trdSignal = TMath::Mean(cntRemaining, trdSlicesRemaining);
   Double_t mom = track->GetOuterParam() ? track->GetOuterParam()->P() : -1;
-  AliDebug(1, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
-  if(IsQAon() && mom > 0.) FillHistogramsTRDSignalV2(trdSignal, mom, mcPID);
+  AliDebug(3, Form("PID Meth. 2: p[%f], TRDSignal[%f]", mom, trdSignal));
+  if(IsQAon() && mom > 0.) FillHistogramsTRDSignal(trdSignal, mom, mcPID, 2);
   return trdSignal;
 }
 
 //___________________________________________________________________
-void AliHFEpidTRD::FillHistogramsTRDSignalV1(Double_t signal, Double_t mom, Int_t species){
-  //
-  // Fill histograms for TRD Signal from Method 1 vs. p for different particle species
-  //
-  (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV1)))->Fill(mom, signal);
-  if(species < 0 || species >= AliPID::kSPECIES) return;
-  (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + species)))->Fill(mom, signal);
-}
-
-//___________________________________________________________________
-void AliHFEpidTRD::FillHistogramsTRDSignalV2(Double_t signal, Double_t mom, Int_t species){
+void AliHFEpidTRD::FillHistogramsTRDSignal(Double_t signal, Double_t p, Int_t species, UInt_t version){
   //
   // Fill histograms for TRD Signal from Method 2 vs. p for different particle species
-  //
-  (dynamic_cast<TH2F *>(fContainer->At(kHistTRDSigV2)))->Fill(mom, signal);  
-  if(species < 0 || species >= AliPID::kSPECIES) return;
-  (dynamic_cast<TH2F *>(fContainer->At(kHistOverallSpecies + AliPID::kSPECIES + species)))->Fill(mom, signal);
+  // Non-standard QA content
+  //
+  if(version == 0 || version > 2) return;
+  if(species >= AliPID::kSPECIES || species < -1) species = -1;
+  Double_t content[3] = {species, p, signal};
+  Char_t histname[256]; sprintf(histname, "fTRDsignalV%d", version);
+  THnSparseF *hsig = dynamic_cast<THnSparseF *>(fContainer->Get(histname));
+  if(hsig) hsig->Fill(content);
 }
 
 //___________________________________________________________________
@@ -324,48 +364,54 @@ void AliHFEpidTRD::AddQAhistograms(TList *l){
   //    - TRD Signal from Meth. 1 vs p
   //    - TRD Signal from Meth. 2 vs p
   //
-  const Int_t kMomentumBins = 41;
+  const Int_t kMomentumBins = 100;
   const Double_t kPtMin = 0.1;
-  const Double_t kPtMax = 10.;
-  const Int_t kSigBinsMeth1 = 100; 
-  const Int_t kSigBinsMeth2 = 100;
-  const Double_t kMinSig = 0.;
-  const Double_t kMaxSigMeth1 = 10000.;
-  const Double_t kMaxSigMeth2 = 10000.;
+  const Double_t kPtMax = 20.;
   
-  if(!fContainer) fContainer = new TList;
-  fContainer->SetName("fTRDqaHistograms");
+  if(!fContainer) fContainer = new AliHFEcollection("fQAhistosTRD", "TRD QA histos");
+  fContainer->CreateTH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
+  fContainer->CreateTH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut; p [GeV/c]; TRD Electron Likelihood", kMomentumBins, kPtMin, kPtMax, 1000, 0, 1);
+  fContainer->CreateTH2F("fTRDthresholds", "TRD Electron Thresholds; p [GeV/c]; Thresholds", kMomentumBins, kPtMin, kPtMax, 1000, 0., 1.);
+  fContainer->CreateTH2F("fTPCsignalBefore", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
+  fContainer->CreateTH2F("fTPCsignalAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; TPC Signal [a.u.]", kMomentumBins, kPtMin, kPtMax, 100, 0., 100.);
+  fContainer->CreateTH2F("fTPCsigmaBefore", "TPC dEdx before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);
+  fContainer->CreateTH2F("fTPCsigmaAfter", "TPC dEdx Spectrum before TRD PID; p [GeV/c]; Normalized TPC distance to the electron line [n#sigma]", kMomentumBins, kPtMin, kPtMax, 100, -10., 10.);
 
-  Double_t momentumBins[kMomentumBins];
-  for(Int_t ibin = 0; ibin < kMomentumBins; ibin++)
-    momentumBins[ibin] = static_cast<Double_t>(TMath::Power(10,TMath::Log10(kPtMin) + (TMath::Log10(kPtMax)-TMath::Log10(kPtMin))/(kMomentumBins-1)*static_cast<Double_t>(ibin)));
-  // Likelihood Histograms
-  fContainer->AddAt(new TH2F("fTRDlikeBefore", "TRD Electron Likelihood before cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeBefore);
-  fContainer->AddAt(new TH2F("fTRDlikeAfter", "TRD Electron Likelihood after cut", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDlikeAfter);
-  fContainer->AddAt(new TH2F("fTRDthesholds", "TRD Electron thresholds", kMomentumBins - 1, momentumBins, 1000, 0, 1), kHistTRDthresholds);
-  // Signal Histograms
-  fContainer->AddAt(new TH2F("fTRDSigV1all", "TRD Signal (all particles, Method 1)", kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistTRDSigV1);
-  fContainer->AddAt(new TH2F("fTRDSigV2all", "TRD Signal (all particles, Method 2)", kMomentumBins - 1, momentumBins, kSigBinsMeth2, kMinSig, kMaxSigMeth2), kHistTRDSigV2);
-  for(Int_t ispec = 0; ispec < AliPID::kSPECIES; ispec++){
-    fContainer->AddAt(new TH2F(Form("fTRDSigV1%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 1)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth1), kHistOverallSpecies + ispec);
-    fContainer->AddAt(new TH2F(Form("fTRDSigV2%s", AliPID::ParticleName(ispec)), Form("TRD Signal (%s, Method 2)", AliPID::ParticleName(ispec)), kMomentumBins - 1, momentumBins, kSigBinsMeth1, kMinSig, kMaxSigMeth2), kHistOverallSpecies + AliPID::kSPECIES + ispec);
-  }
-  l->AddLast(fContainer);
+  // Monitor THnSparse for TRD Signal
+  const Int_t kBinsTRDsignal = 3; 
+  Int_t nBins[kBinsTRDsignal] = {AliPID::kSPECIES +1, kMomentumBins, 100};
+  Double_t binMin[kBinsTRDsignal] = {-1, kPtMin, 0};
+  Double_t binMax[kBinsTRDsignal] = {AliPID::kSPECIES, kPtMax, 1000};
+  fContainer->CreateTHnSparse("fTRDsignalV1", "TRD Signal V1", kBinsTRDsignal, nBins, binMin, binMax);
+  fContainer->CreateTHnSparse("fTRDsignalV2", "TRD Signal V2", kBinsTRDsignal, nBins, binMin, binMax);
+
+  // Make logatithmic binning
+  TString hnames[7] = {"fTRDlikeBefore", "fTRDlikeAfter", "fTRDthresholds", "fTRDsignalBefore", "fTRDsignalAfter", "fTPCsigmaBefore", "fTPCsigmaAfter"};
+  for(Int_t ihist = 0; ihist < 7; ihist++)
+    fContainer->BinLogAxis(hnames[ihist].Data(), 0);
+  fContainer->BinLogAxis("fTRDsignalV1", 1);
+  fContainer->BinLogAxis("fTRDsignalV2", 1);
+
+  l->Add(fContainer->GetList());
 }
 
 //___________________________________________________________________
-void AliHFEpidTRD::FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb){
+void AliHFEpidTRD::FillStandardQA(Int_t whenFilled, AliESDtrack *esdTrack){
   //
   // Fill Likelihood Histogram before respectively after decision
   //
-  TH2F *histo = NULL;
+  Double_t p =  esdTrack->P();
+  Double_t like[AliPID::kSPECIES];
+  esdTrack->GetTRDpid(like);
+  TString step;
   if(whenFilled)
-    histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeAfter));
+    step = "After";
   else
-    histo = dynamic_cast<TH2F *>(fContainer->At(kHistTRDlikeBefore));
-  if(!histo){
-    AliError("QA histograms not found");
-    return;
-  }
-  histo->Fill(p, elProb);
+    step = "Before";
+  TString histos[3] = {"fTRDlike", "fTPCsignal", "fTPCsigma"};
+  for(Int_t ihist = 0; ihist < 3; ihist++) histos[ihist] += step;
+  fContainer->Fill(histos[0].Data(), esdTrack->P(), like[AliPID::kElectron]);
+  fContainer->Fill(histos[1].Data(), p, esdTrack->GetTPCsignal());
+  if(fESDpid)
+    fContainer->Fill(histos[2].Data(), p, fESDpid->NumberOfSigmasTPC(esdTrack, AliPID::kElectron));
 }
index 60296b3..1c5d66e 100644 (file)
@@ -27,6 +27,7 @@
 class AliAODTrack;
 class AliAODMCParticle;
 class AliESDtrack;
+class AliHFEcollection;
 class AliMCParticle;
 class AliVParticle;
 class TList;
@@ -62,26 +63,33 @@ class AliHFEpidTRD : public AliHFEpidBase{
     Double_t GetTRDSignalV1(AliESDtrack *track, Int_t mcPID);
     Double_t GetTRDSignalV2(AliESDtrack *track, Int_t mcPID);
 
+    Bool_t IsCalculateTRDSignals() const { return TestBit(kTRDsignals); }
     void SetPIDMethod(PIDMethodTRD_t method) { fPIDMethod = method; };
+    void SetMinP(Double_t p) { fMinP = p; }
+    void CalculateTRDSignals(Bool_t docalc) { SetBit(kTRDsignals, docalc); } 
 
     Double_t GetTRDthresholds(Double_t electronEff, Double_t p);
   protected:
+    enum{
+      kTRDsignals = BIT(16)
+    };
     void Copy(TObject &ref) const;
     Int_t MakePIDesd(AliESDtrack *esdTrack, AliMCParticle *mcTrack);
     Int_t MakePIDaod(AliAODTrack *aofTrack, AliAODMCParticle *aodMC);
     Int_t GetMCpid(AliESDtrack *track);
     void InitParameters();
+    void InitParameters1DLQ();
     virtual void AddQAhistograms(TList *l);
     void GetParameters(Double_t electronEff, Double_t *parameters);
 
-    void FillHistogramsLikelihood(Int_t whenFilled, Float_t p, Float_t elProb);
-    void FillHistogramsTRDSignalV1(Double_t signal, Double_t p, Int_t species);
-    void FillHistogramsTRDSignalV2(Double_t signal, Double_t p, Int_t species);
+    void FillStandardQA(Int_t whenFilled, AliESDtrack *esdTrack);
+    void FillHistogramsTRDSignal(Double_t signal, Double_t p, Int_t species, UInt_t version);
   private:
     static const Double_t fgkVerySmall;                       // Check for 0
+    Double_t fMinP;                                         // Minimum momentum above which TRD PID is applied
     PIDMethodTRD_t fPIDMethod;                              // PID Method: 2D Likelihood or Neural Network
     Double_t fThreshParams[kThreshParams];                  // Threshold parametrisation
-    TList *fContainer;                                      // QA  Histogram Container
+    AliHFEcollection *fContainer;                                      // QA  Histogram Container
   ClassDef(AliHFEpidTRD, 1)     // TRD electron ID class
 };
 
index 3c73c1a..4cc3dc3 100644 (file)
@@ -239,13 +239,15 @@ void AliHFEsecVtx::Process(AliVTrack *signalTrack){
 }
 
 //_______________________________________________________________________________________________
-void AliHFEsecVtx::CreateHistograms(TList *qaList)
+void AliHFEsecVtx::CreateHistograms(TList * const qaList)
 { 
   //
   // create histograms
   //
+  
+  if(!qaList) return;
 
-  fSecVtxList = new TList;
+  fSecVtxList = qaList;
   fSecVtxList->SetName("SecVtx");
 
   MakeContainer();
@@ -267,7 +269,7 @@ void AliHFEsecVtx::CreateHistograms(TList *qaList)
   }
   */
 
-  qaList->AddLast(fSecVtxList);
+  //qaList->AddLast(fSecVtxList);
 }
 
 //_______________________________________________________________________________________________
index 45dde23..cfac67b 100644 (file)
@@ -52,7 +52,7 @@ class AliHFEsecVtx : public TObject {
     AliHFEsecVtx &operator=(const AliHFEsecVtx &); // assignment operator
     virtual ~AliHFEsecVtx();
 
-    void CreateHistograms(TList *qaList);
+    void CreateHistograms(TList * const qaList);
 
     void Process(AliVTrack *track);
 
index f7632e2..91af78a 100644 (file)
@@ -127,7 +127,7 @@ Bool_t AliHFEtools::BinLogAxis(TObject *o, Int_t dim){
     newBins[i] = factor * newBins[i-1];
   }
   axis->Set(bins, newBins);
-  delete [] newBins;
+  delete newBins;
 
   return kTRUE;
 }