Update of the HFE package
authorsma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jul 2011 10:31:55 +0000 (10:31 +0000)
committersma <sma@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 14 Jul 2011 10:31:55 +0000 (10:31 +0000)
48 files changed:
PWG3/hfe/AliAnalysisTaskDCA.cxx
PWG3/hfe/AliAnalysisTaskDisplacedElectrons.cxx
PWG3/hfe/AliAnalysisTaskHFE.cxx
PWG3/hfe/AliAnalysisTaskHFE.h
PWG3/hfe/AliAnalysisTaskHFEpidQA.cxx
PWG3/hfe/AliAnalysisTaskHFEpidQA.h
PWG3/hfe/AliHFEV0cuts.cxx
PWG3/hfe/AliHFEV0cuts.h
PWG3/hfe/AliHFEV0pid.cxx
PWG3/hfe/AliHFEV0pid.h
PWG3/hfe/AliHFEV0pidMC.cxx
PWG3/hfe/AliHFEV0pidMC.h
PWG3/hfe/AliHFEcollection.cxx
PWG3/hfe/AliHFEcollection.h
PWG3/hfe/AliHFEcuts.cxx
PWG3/hfe/AliHFEcuts.h
PWG3/hfe/AliHFEemcalPIDqa.cxx
PWG3/hfe/AliHFEemcalPIDqa.h
PWG3/hfe/AliHFEextraCuts.cxx
PWG3/hfe/AliHFEextraCuts.h
PWG3/hfe/AliHFEmcQA.cxx
PWG3/hfe/AliHFEpid.cxx
PWG3/hfe/AliHFEpid.h
PWG3/hfe/AliHFEpidBase.h
PWG3/hfe/AliHFEpidEMCAL.cxx
PWG3/hfe/AliHFEpidEMCAL.h
PWG3/hfe/AliHFEpidESD.h
PWG3/hfe/AliHFEpidQA.cxx
PWG3/hfe/AliHFEpidQA.h
PWG3/hfe/AliHFEpidTPC.cxx
PWG3/hfe/AliHFEpidTPC.h
PWG3/hfe/AliHFEpidTRD.cxx
PWG3/hfe/AliHFEpidTRD.h
PWG3/hfe/AliHFEsecVtx.cxx
PWG3/hfe/AliHFEsecVtx.h
PWG3/hfe/AliHFEsignalCuts.cxx
PWG3/hfe/AliHFEspectrum.cxx
PWG3/hfe/AliHFEspectrum.h
PWG3/hfe/AliHFEtaggedTrackAnalysis.cxx
PWG3/hfe/AliHFEtaggedTrackAnalysis.h
PWG3/hfe/AliHFEtpcPIDqa.cxx
PWG3/hfe/AliHFEtpcPIDqa.h
PWG3/hfe/AliHFEtrackFilter.cxx
PWG3/hfe/AliHFEtrackFilter.h
PWG3/hfe/AliHFEtrdPIDqa.cxx
PWG3/hfe/AliHFEtrdPIDqa.h
PWG3/hfe/AliHFEtrdPIDqaV1.cxx
PWG3/hfe/AliHFEtrdPIDqaV1.h

index 2480a44..96f1e6c 100644 (file)
@@ -336,7 +336,8 @@ void AliAnalysisTaskDCA::UserCreateOutputObjects(){
     fHFEpid->AddDetector("TOF", 0);
     fHFEpid->AddDetector("TPC", 1);
     cout<<endl<<" ---> TPC and TOF added to the PID"<<endl;
-    fHFEpid->ConfigureTPCrejection();
+               fHFEpid->ConfigureTOF();
+    fHFEpid->ConfigureTPCdefaultCut();
     fHFEpid->InitializePID();
   }
 
index 832a7d8..2e25a13 100644 (file)
@@ -255,7 +255,8 @@ void AliAnalysisTaskDisplacedElectrons::UserCreateOutputObjects(){
   fDePID->SetHasMCData(HasMCData());
   fDePID->AddDetector("TPC", 0);
   fDePID->AddDetector("TOF", 1);
-  fDePID->ConfigureTPCrejection();
+       fDePID->ConfigureTOF();
+  fDePID->ConfigureTPCdefaultCut();
   fDePID->InitializePID();     // Only restrictions to TPC allowed   
 
 
index 4dd843b..9eba095 100644 (file)
@@ -57,6 +57,7 @@
 #include "AliMCEvent.h"
 #include "AliMCEventHandler.h"
 #include "AliMCParticle.h"
+#include "AliMultiplicity.h"
 #include "AliPID.h"
 #include "AliStack.h"
 #include "AliTriggerAnalysis.h"
@@ -87,6 +88,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(kTRUE)
+  , fFillNoCuts(kFALSE)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -95,7 +97,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE():
   , fHasSpecialTriggerSelection(kFALSE)
   , fRejectKinkMother(kTRUE)
   , fSpecialTrigger("NONE")
-  , fCentralityF(99.0)
+  , fCentralityF(-1)
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
@@ -134,6 +136,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(kTRUE)
+  , fFillNoCuts(kFALSE)
   , fBackGroundFactorApply(kFALSE)
   , fRemovePileUp(kFALSE)
   , fIdentifiedAsPileUp(kFALSE)
@@ -142,7 +145,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const char * name):
   , fHasSpecialTriggerSelection(kFALSE)
   , fRejectKinkMother(kTRUE)
   , fSpecialTrigger("NONE")
-  , fCentralityF(99.0)
+  , fCentralityF(-1)
   , fContributors(0.5)
   , fWeightBackGround(0.)
   , fVz(0.0)
@@ -190,6 +193,7 @@ AliAnalysisTaskHFE::AliAnalysisTaskHFE(const AliAnalysisTaskHFE &ref):
   , fQAlevel(0)
   , fPlugins(0)
   , fFillSignalOnly(ref.fFillSignalOnly)
+  , fFillNoCuts(ref.fFillNoCuts)
   , fBackGroundFactorApply(ref.fBackGroundFactorApply)
   , fRemovePileUp(ref.fRemovePileUp)
   , fIdentifiedAsPileUp(ref.fIdentifiedAsPileUp)
@@ -251,6 +255,7 @@ void AliAnalysisTaskHFE::Copy(TObject &o) const {
   target.fQAlevel = fQAlevel;
   target.fPlugins = fPlugins;
   target.fFillSignalOnly = fFillSignalOnly;
+  target.fFillNoCuts = fFillNoCuts;
   target.fBackGroundFactorApply = fBackGroundFactorApply;
   target.fRemovePileUp = fRemovePileUp;
   target.fIdentifiedAsPileUp = fIdentifiedAsPileUp;
@@ -347,28 +352,8 @@ void AliAnalysisTaskHFE::UserCreateOutputObjects(){
   // First Part: Make QA histograms
   fQACollection = new AliHFEcollection("TaskQA", "QA histos from the Electron Task");
   fQACollection->CreateTH1F("nElectronTracksEvent", "Number of Electron Candidates", 100, 0, 100);
-  fQACollection->CreateProfile("conr", "Electron PID contamination", 20, 0, 20);
-  fQACollection->CreateTH1F("alpha_rec", "Alpha from reconstructed tracks with TRD hits", 36, -TMath::Pi(), TMath::Pi());
-  fQACollection->CreateTH1F("alpha_sim", "Alpha from simulated electron tracks", 36, -TMath::Pi(), TMath::Pi());
   fQACollection->CreateTH1F("nElectron", "Number of electrons", 100, 0, 100);
-  fQACollection->CreateProfile("pidquality", "TRD PID quality as function of momentum", 20, 0, 20);
-  fQACollection->CreateProfile("ntrdclusters", "Number of TRD clusters as function of momentum", 20, 0, 20);
-  fQACollection->CreateTH1F("chi2TRD","#chi2 per TRD cluster", 20, 0, 20);
-  fQACollection->CreateTH1F("mccharge", "MC Charge", 200, -100, 100);
   fQACollection->CreateTH2F("radius", "Production Vertex", 100, 0.0, 5.0, 100, 0.0, 5.0);
-  // Temporary histograms for TPC number of clusters for all signal tracks (MC true electrons) and for selected tracks (Markus Fasel)
-  fQACollection->CreateTH2F("TPCclusters2_1_Signal", "TPCclusterInfo for findable clusters for 2 neighbors for signal tracks", 30, 0.1, 10., 162, 0., 161.);
-  fQACollection->CreateTH2F("TPCclusters2_0_Signal", "TPCclusterInfo for the ratio for 2 neighbors for signal tracks", 30, 0.1, 10., 100, 0., 1.);
-  fQACollection->CreateTH2F("TPCclusters2_1_Selected", "TPCclusterInfo for findable clusters for 2 neighbors for selected tracks", 30, 0.1, 10., 162, 0., 161.);
-  fQACollection->CreateTH2F("TPCclusters2_0_Selected", "TPCclusterInfo for the ratio for 2 neighbors for selected tracks", 30, 0.1, 10., 110, 0., 1.1);
-  fQACollection->CreateTH2F("TPCncls_Signal", "TPC Number of clusters for signal tracks", 30, 0.1, 10., 162, 0., 161.);
-  fQACollection->CreateTH2F("TPCclr_Signal", "TPC cluster ratio for signal tracks", 30, 0.1, 10., 110, 0., 1.1);
-  fQACollection->BinLogAxis("TPCclusters2_1_Signal", 0); 
-  fQACollection->BinLogAxis("TPCclusters2_0_Signal", 0);
-  fQACollection->BinLogAxis("TPCclusters2_1_Selected", 0); 
-  fQACollection->BinLogAxis("TPCclusters2_0_Selected", 0);
-  fQACollection->BinLogAxis("TPCncls_Signal", 0); 
-  fQACollection->BinLogAxis("TPCclr_Signal", 0);
   // Temporary histograms for TRD efficiency maps (Markus Fasel)
   fQACollection->CreateTH2F("TRDmapPosBefore", "Pos. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
   fQACollection->CreateTH2F("TRDmapNegBefore", "Neg. charged tracks before TRD; #eta; #phi", 100, -1., 1., 180, 0., 2*TMath::Pi());
@@ -520,8 +505,8 @@ void AliAnalysisTaskHFE::UserExec(Option_t *){
   }
 
   // need the centrality for everything (MC also)
-  fCentralityF = -100.0;
-  if(!ReadCentrality()) fCentralityF = -100.0;
+  fCentralityF = -1;
+  if(!ReadCentrality()) fCentralityF = -1;
   //printf("pass centrality\n");
   //printf("Reading fCentralityF %f\n",fCentralityF);
   
@@ -657,7 +642,7 @@ void AliAnalysisTaskHFE::ProcessMC(){
   fCFM->GetEventContainer()->Fill(eventContainer,AliHFEcuts::kEventStepGenerated);
   Int_t nElectrons = 0;
   if(IsESDanalysis()){
-   if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
+   if(!((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0))){ //kStepMCGeneratedZOutNoPileUpCentralityFine
     if (HasMCData() && IsQAOn(kMCqa)) {
       AliDebug(2, "Running MC QA");
 
@@ -732,6 +717,8 @@ void AliAnalysisTaskHFE::ProcessESD(){
   if(fTaggedTrackAnalysis) {
     fTaggedTrackAnalysis->SetMagneticField(fESD->GetMagneticField());
     fTaggedTrackAnalysis->SetCentrality(fCentralityF);
+    if(IsPbPb()) fTaggedTrackAnalysis->SetPbPb();
+    else fTaggedTrackAnalysis->SetPP();
   }
 
   // Do event Normalization
@@ -762,7 +749,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecNoPileUp);
 
   //
-  if(TMath::Abs(fCentralityF+100.0) < 0.01) return; 
+  if(TMath::Abs(fCentralityF) < 0) return; 
   fCFM->GetEventContainer()->Fill(eventContainer, AliHFEcuts::kEventStepRecCentralityOk);
   //printf("In ProcessESD %f\n",fCentralityF);
 
@@ -803,6 +790,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
   Int_t nElectronCandidates = 0;
   AliESDtrack *track = NULL, *htrack = NULL;
   AliMCParticle *mctrack = NULL;
+  TParticle* mctrack4QA = NULL;
   Int_t pid = 0;
 
   Bool_t signal = kTRUE;
@@ -852,6 +840,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
     if(HasMCData()){
       // Check if it is electrons near the vertex
       if(!(mctrack = dynamic_cast<AliMCParticle *>(fMCEvent->GetTrack(TMath::Abs(track->GetLabel()))))) continue;
+      mctrack4QA = mctrack->Particle();
 
       if(fFillSignalOnly && !fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, mctrack)) signal = kFALSE; 
       else AliDebug(3, "Signal Electron");
@@ -859,30 +848,16 @@ void AliAnalysisTaskHFE::ProcessESD(){
     // Cache new Track information inside the var manager
     fVarManager->NewTrack(track, mctrack, fCentralityF, -1, signal);
 
-    if(signal) {
-      fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
-      fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
-      if((track->GetStatus() & AliESDtrack::kTPCout) 
-          && (TMath::Abs(track->Eta()) < 0.8) 
-          && (track->GetKinkIndex(0) == 0)){
-        fQACollection->Fill("TPCclusters2_1_Signal", track->Pt(), track->GetTPCClusterInfo(2,1));
-        fQACollection->Fill("TPCclusters2_0_Signal", track->Pt(), track->GetTPCNclsF() > 0 ?  track->GetTPCClusterInfo(2,1)/track->GetTPCNclsF() : 0.);
-        fQACollection->Fill("TPCncls_Signal", track->Pt(), track->GetTPCNcls());
-        fQACollection->Fill("TPCclr_Signal", track->Pt(), track->GetTPCNclsF() > 0 ? static_cast<Double_t>(track->GetTPCNcls())/static_cast<Double_t>(track->GetTPCNclsF()) : 0.);
+    if(fFillNoCuts) {
+      if(signal || !fFillSignalOnly){
+        fVarManager->FillContainer(fContainer, "recTrackContReco", AliHFEcuts::kStepRecNoCut, kFALSE);
+        fVarManager->FillContainer(fContainer, "recTrackContMC", AliHFEcuts::kStepRecNoCut, kTRUE);
       }
     }
 
     // RecKine: ITSTPC cuts  
     if(!ProcessCutStep(AliHFEcuts::kStepRecKineITSTPC, track)) continue;
     
-    // Check TRD criterions (outside the correction framework)
-    if(track->GetTRDncls()){
-      fQACollection->Fill("chi2TRD", track->GetTRDchi2()/track->GetTRDncls());
-      fQACollection->Fill("alpha_rec", track->GetAlpha());
-      fQACollection->Fill("pidquality", container[0], track->GetTRDpidQuality());
-      fQACollection->Fill("ntrdclusters", container[0], track->GetTRDncls());
-    }
-
     
     // RecPrim
     if(fRejectKinkMother) { 
@@ -928,11 +903,11 @@ void AliAnalysisTaskHFE::ProcessESD(){
     hfetrack.SetRecTrack(track);
     if(HasMCData()) hfetrack.SetMCTrack(mctrack);
     hfetrack.SetCentrality(fCentralityF);
+    if(IsPbPb()) hfetrack.SetPbPb();
+    else hfetrack.SetPP();
     fPID->SetVarManager(fVarManager);
     if(!fPID->IsSelected(&hfetrack, fContainer, "recTrackCont", fPIDqa)) continue;
     nElectronCandidates++;
-    fQACollection->Fill("TPCclusters2_1_Selected", track->Pt(), track->GetTPCClusterInfo(2,1));
-    fQACollection->Fill("TPCclusters2_0_Selected", track->Pt(), track->GetTPCClusterInfo(2,0));
 
     // Fill Histogram for Hadronic Background
     if(HasMCData()){
@@ -944,7 +919,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
     if(signal) {
       // Apply weight for background contamination
            if(fBackGroundFactorApply==kTRUE) {
-             if(IsPbPb()) fWeightBackGround =  fBackGroundFactorArray[(Int_t)fCentralityF]->Eval(TMath::Abs(track->P()));
+             if(IsPbPb() && fCentralityF >= 0) fWeightBackGround =  fBackGroundFactorArray[fCentralityF]->Eval(TMath::Abs(track->P()));
              else    fWeightBackGround =  fBackGroundFactorArray[0]->Eval(TMath::Abs(track->P())); // pp case
 
              if(fWeightBackGround < 0.0) fWeightBackGround = 0.0;
@@ -961,7 +936,7 @@ void AliAnalysisTaskHFE::ProcessESD(){
       if(fSecVtx->Process(track) && signal) {
         fVarManager->FillContainer(fContainer, "recTrackContSecvtxReco", AliHFEcuts::kStepHFEcutsSecvtx, kFALSE);
         fVarManager->FillContainer(fContainer, "recTrackContSecvtxMC", AliHFEcuts::kStepHFEcutsSecvtx, kTRUE);
-       bTagged=kTRUE;
+        bTagged=kTRUE;
       }
     }
 
@@ -1075,7 +1050,7 @@ void AliAnalysisTaskHFE::ProcessAOD(){
   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
   if(!fAOD){
     AliError("AOD Event required for AOD Analysis");
-    return;
+      return;
   }
   
   //
@@ -1208,7 +1183,6 @@ Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
   }
 
   if(!fCFM->CheckParticleCuts(AliHFEcuts::kStepMCGenerated, track)) return kFALSE;
-  fQACollection->Fill("mccharge", signalContainer[3]);
   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGenerated, kFALSE);
   signalContainer[4] = 0;
   if(fSignalCuts->IsSelected(track)){
@@ -1233,10 +1207,9 @@ Bool_t AliAnalysisTaskHFE::ProcessMCtrack(AliVParticle *track){
     signalContainer[5] = 2;
   }
   fQACollection->Fill("SignalToBackgroundMC", signalContainer);
-  fQACollection->Fill("alpha_sim", track->Phi() - TMath::Pi());
 
   // Step GeneratedZOutNoPileUp
-  if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (TMath::Abs(fCentralityF+100.0) < 0.01)) return kFALSE;
+  if((fIdentifiedAsPileUp) || (TMath::Abs(fVz) > fCuts->GetVertexRange()) || (fCentralityF < 0)) return kFALSE;
   fVarManager->FillContainer(fContainer, "MCTrackCont", AliHFEcuts::kStepMCGeneratedZOutNoPileUpCentralityFine, kFALSE);
   //printf("In ProcessMCtrack %f\n",fCentralityF);
 
@@ -1544,153 +1517,91 @@ Bool_t AliAnalysisTaskHFE::ReadCentrality() {
   //
   // Recover the centrality of the event from ESD or AOD
   //
- if(IsAODanalysis()){
-
-   AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
-   if(!fAOD){
-     AliError("AOD Event required for AOD Analysis");
-     return kFALSE;
-   }
-
-   if(IsPbPb()) {
-     // Centrality
-     AliCentrality *aodCentrality = fAOD->GetCentrality();
-     Float_t fCentralityFtemp = aodCentrality->GetCentralityPercentile("V0M");
-     
-     if( fCentralityFtemp >=  0. && fCentralityFtemp <  10.) fCentralityF =  0;
-     else if ( fCentralityFtemp >=  10. && fCentralityFtemp <  20.) fCentralityF =  1;
-     else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
-     else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
-     else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
-     else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
-     else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  90.) fCentralityF = 6;
-     else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 7;
-     //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <  95.) fCentralityF = 8;
-     //else if ( fCentralityF_temp >= 95. && fCentralityF_temp <  90.) fCentralityF = 9;
-     //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
-     else return kFALSE;
+  Int_t bin = -1;
+  if(IsPbPb()) {
+    // Centrality
+    AliCentrality *centrality = fInputEvent->GetCentrality();
+    Float_t fCentralityFtemp = centrality->GetCentralityPercentile("V0M");
+    Float_t centralityLimits[12] = {0.,5.,10., 20., 30., 40., 50., 60.,70.,80., 90., 100.};
+    for(Int_t ibin = 0; ibin < 11; ibin++){
+      if(fCentralityFtemp >= centralityLimits[ibin] && fCentralityFtemp < centralityLimits[ibin+1]){
+        bin = ibin;
+        break;
+      }
+    } 
+    if(bin == -1) bin = 11; // Overflow
+  } else {
+    // PP: Tracklet multiplicity, use common definition
+    Int_t itsMultiplicity = GetITSMultiplicity(fInputEvent);
+    Int_t multiplicityLimits[8] = {0, 1, 9, 17, 25, 36, 60, 500};
+    for(Int_t ibin = 0; ibin < 7; ibin++){  
+      if(itsMultiplicity >= multiplicityLimits[ibin] && itsMultiplicity < multiplicityLimits[ibin + 1]){
+        bin = ibin;
+        break;
+      }
+    }
+    if(bin == -1) bin = 7;  // Overflow
+  }
+  fCentralityF = bin;
+  AliDebug(2, Form("Centrality class %d\n", fCentralityF));
  
-     // contributors
-     fContributors = 0.5;
-     Int_t contributorstemp = 0;
-     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
-     if(vtxAOD) contributorstemp = vtxAOD->GetNContributors();
-     
-     //printf("PbPb contributors_temp %d\n",contributors_temp);
-     
-     if( contributorstemp <=  0) fContributors =  0.5;
-     else fContributors = 1.5;   
-   
-
-
-   }
-   else {
-     fCentralityF = 0;
-     Int_t centralityFtemp = 0;
-     const AliAODVertex *vtxAOD = fAOD->GetPrimaryVertex();
-     if(vtxAOD) centralityFtemp = vtxAOD->GetNContributors();
-     
-     //printf("pp centralityF_temp %d\n",centralityF_temp);
-     
-     if( centralityFtemp <=  0) fCentralityF =  0;
-     else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
-     else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
-     else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
-     else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
-     else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
-     else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
-     else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
-     else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
-     else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
-     else if ( centralityFtemp >= 50) fCentralityF = 10;
-     
-   }
-
-   return kTRUE;
-
-   
- } else {
-   
-   AliDebug(3, "Processing ESD Centrality");
-   AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
-   if(!fESD){
-     AliError("ESD Event required for ESD Analysis");
-     return kFALSE;
-   }
-   const char* type = fESD->GetBeamType();
-
-   if (strstr(type,"Pb-Pb")) {
-   
-   // Centrality
-   AliCentrality *esdCentrality = fESD->GetCentrality();
-   Float_t fCentralityFtemp = esdCentrality->GetCentralityPercentile("V0M");
-   //printf("PbPb fCentralityF_temp %f\n",fCentralityF_temp);
-
-   if( fCentralityFtemp >=  0. && fCentralityFtemp <   10.) fCentralityF =  0;
-   else if ( fCentralityFtemp >= 10. && fCentralityFtemp <  20.) fCentralityF =  1;
-   else if ( fCentralityFtemp >= 20. && fCentralityFtemp <  30.) fCentralityF = 2;
-   else if ( fCentralityFtemp >= 30. && fCentralityFtemp <  40.) fCentralityF = 3;
-   else if ( fCentralityFtemp >= 40. && fCentralityFtemp <  50.) fCentralityF = 4;
-   else if ( fCentralityFtemp >= 50. && fCentralityFtemp <  60.) fCentralityF = 5;
-   else if ( fCentralityFtemp >= 60. && fCentralityFtemp <  80.) fCentralityF = 6;
-   else if ( fCentralityFtemp >= 80. && fCentralityFtemp <  90.) fCentralityF = 7;
-   else if ( fCentralityFtemp >= 90. && fCentralityFtemp <=  100.) fCentralityF = 8;
-   //else if ( fCentralityF_temp >= 80. && fCentralityF_temp <  90.) fCentralityF = 9;
-   //else if ( fCentralityF_temp >= 90. && fCentralityF_temp <=100.) fCentralityF = 10;
-   else return kFALSE;
-
-   //   Float_t fCentralityF_temp10 = esdCentrality->GetCentralityClass10("V0M");
-   //   printf("PbPb fCentralityF_temp %f %f %f \n",fCentralityF_temp, fCentralityF_temp10, fCentralityF);
-
-   // contributors
-   fContributors = 0.5;
-   Int_t contributorstemp = 0;
-   const AliESDVertex *vtxESD = fESD->GetPrimaryVertexSPD();
-   if(vtxESD) contributorstemp = vtxESD->GetNContributors();
-   
-   //printf("PbPb contributors_temp %d\n",contributors_temp);
-   
-   if( contributorstemp <=  0) fContributors =  0.5;
-   else fContributors = 1.5;   
-   
-   return kTRUE;
-
-   }
-
-   
-   if (strstr(type,"p-p")) {
-     fCentralityF = 0;
-     Int_t centralityFtemp = 0;
-     const AliESDVertex *vtxESD = fESD->GetPrimaryVertexTracks();
-     if(vtxESD && vtxESD->GetStatus()) centralityFtemp = vtxESD->GetNContributors();
-     
-     //printf("pp centralityF_temp %d\n",centralityF_temp);
-     
-     if( centralityFtemp <=  0) fCentralityF =  0;
-     else if ( centralityFtemp >  0 && centralityFtemp <  2) fCentralityF = 1;
-     else if ( centralityFtemp >=  2 && centralityFtemp <  3) fCentralityF = 2;
-     else if ( centralityFtemp >= 3 && centralityFtemp <  4) fCentralityF = 3;
-     else if ( centralityFtemp >= 4 && centralityFtemp <  5) fCentralityF = 4;
-     else if ( centralityFtemp >= 5 && centralityFtemp <  10) fCentralityF = 5;
-     else if ( centralityFtemp >= 10 && centralityFtemp <  20) fCentralityF = 6;
-     else if ( centralityFtemp >= 20 && centralityFtemp <  30) fCentralityF = 7;
-     else if ( centralityFtemp >= 30 && centralityFtemp <  40) fCentralityF = 8;
-     else if ( centralityFtemp >= 40 && centralityFtemp <  50) fCentralityF = 9;
-     else if ( centralityFtemp >= 50) fCentralityF = 10;
-    
-     return kTRUE; 
-     
-   }
-
-   return kFALSE;
-
-  //printf("centrality %f\n",fCentralityF);
-   
- }
+  // contributors, to be outsourced
+  const AliVVertex *vtx;
+  if(IsAODanalysis()){
+    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
+    if(!fAOD){
+      AliError("AOD Event required for AOD Analysis");
+      return kFALSE;
+    }
+    vtx = fAOD->GetPrimaryVertex();
+  } else {
+    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
+    if(!fESD){
+      AliError("ESD Event required for ESD Analysis");
+      return kFALSE;
+    }
+    vtx = fESD->GetPrimaryVertexSPD();
+  }
+  if(!vtx){ 
+    fContributors = 0.5;
+    return kFALSE;
+  }
+  else {
+    Int_t contributorstemp = vtx->GetNContributors();
+    if( contributorstemp <=  0) fContributors =  0.5;
+    else fContributors = 1.5;
+  }
+  return kTRUE;
+}
 
- //printf("centrality %f\n",fCentralityF);
+//___________________________________________________
+Int_t AliAnalysisTaskHFE::GetITSMultiplicity(AliVEvent *ev){
+  //
+  // Definition of the Multiplicity according to the JPSI group (F. Kramer)
+  //
+  Int_t nTracklets = 0;
+  Int_t nAcc = 0;
+  Double_t etaRange = 1.6;
+
+  if (ev->IsA() == AliAODEvent::Class()) {
+    AliAODTracklets *tracklets = ((AliAODEvent*)ev)->GetTracklets();
+    nTracklets = tracklets->GetNumberOfTracklets();
+    for (Int_t nn = 0; nn < nTracklets; nn++) {
+      Double_t theta = tracklets->GetTheta(nn);
+      Double_t eta = -TMath::Log(TMath::Tan(theta/2.0));
+      if (TMath::Abs(eta) < etaRange) nAcc++;
+    }
+  } else if (ev->IsA() == AliESDEvent::Class()) {
+    nTracklets = ((AliESDEvent*)ev)->GetMultiplicity()->GetNumberOfTracklets();
+    for (Int_t nn = 0; nn < nTracklets; nn++) {
+       Double_t eta = ((AliESDEvent*)ev)->GetMultiplicity()->GetEta(nn);
+      if (TMath::Abs(eta) < etaRange) nAcc++;
+    }
+  } else return -1;
 
+  return nAcc;
 }
+
 //___________________________________________________
 void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
   //
@@ -1701,7 +1612,7 @@ void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
    AliAODEvent *fAOD = dynamic_cast<AliAODEvent *>(fInputEvent);
    if(!fAOD){
      AliError("AOD Event required for AOD Analysis");
-     return;
+       return;
    }
    // PileUp
    if(fRemovePileUp && fAOD->IsPileupFromSPD()) fIdentifiedAsPileUp = kTRUE; 
@@ -1718,7 +1629,7 @@ void AliAnalysisTaskHFE::RejectionPileUpVertexRangeEventCut() {
    AliESDEvent *fESD = dynamic_cast<AliESDEvent *>(fInputEvent);
    if(!fESD){
      AliError("ESD Event required for ESD Analysis");
-     return;
+       return;
    }
    // PileUp
    fIdentifiedAsPileUp = kFALSE;
index c3818ad..59258ad 100644 (file)
@@ -107,6 +107,7 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     void SwitchOnPlugin(Int_t plug);
     void SetHasMCData(Bool_t hasMC = kTRUE) { SetBit(kHasMCdata, hasMC); };
     void SetFillSignalOnly(Bool_t signalOnly) { fFillSignalOnly = signalOnly; }
+    void SetFillNoCuts(Bool_t fillNoCuts) { fFillNoCuts = fillNoCuts; }
     void SetRemovePileUp(Bool_t removePileUp) { fRemovePileUp = removePileUp; }
     void SetPIDPreselect(AliHFEpid * const cuts) { fPIDpreselect = cuts; };
     void SetAODAnalysis() { SetBit(kAODanalysis, kTRUE); };
@@ -138,12 +139,14 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     void ProcessMC();
     void ProcessESD();
     void ProcessAOD();
+    Int_t GetITSMultiplicity(AliVEvent *ev);
     Bool_t PreSelectTrack(AliESDtrack *track) const;
     Bool_t ProcessMCtrack(AliVParticle *track);
     Bool_t ProcessCutStep(Int_t cutStep, AliVParticle *track);
     ULong_t fQAlevel;                     // QA level
     UShort_t fPlugins;                    // Enabled Plugins
     Bool_t fFillSignalOnly;               // Fill container only with MC Signal Tracks
+    Bool_t fFillNoCuts;                   // Fill container before any cut
     Bool_t fBackGroundFactorApply;        // Apply Background Function Subtraction
     Bool_t fRemovePileUp;                 // Remove Pile Up
     Bool_t fIdentifiedAsPileUp;           // Identified as pile-up
@@ -152,7 +155,7 @@ class AliAnalysisTaskHFE : public AliAnalysisTaskSE{
     Bool_t fHasSpecialTriggerSelection;   // Select special triggered events
     Bool_t fRejectKinkMother;             // Reject Kink Mother
     TString fSpecialTrigger;              // Special trigger selection
-    Float_t fCentralityF;                 // Centrality
+    Int_t   fCentralityF;                 // Centrality
     Float_t fContributors;                // Contributors
     Double_t fWeightBackGround;            // weight background function
     Double_t fVz;                         // z position of the primary vertex
index e5573ff..aedf4cd 100644 (file)
@@ -43,6 +43,7 @@ AliAnalysisTaskHFEpidQA::AliAnalysisTaskHFEpidQA():
   , fOutput(NULL)
   , fEvents(NULL)
   , fNNref(NULL)
+  , fTRDTotalChargeInSlice0(kFALSE)
 {
   //
   // Default Constructor
@@ -55,6 +56,7 @@ AliAnalysisTaskHFEpidQA::AliAnalysisTaskHFEpidQA(const Char_t *name):
   , fOutput(NULL)
   , fEvents(NULL)
   , fNNref(NULL)
+  , fTRDTotalChargeInSlice0(kFALSE)
 {
   //
   // Default Constructor
@@ -77,11 +79,13 @@ void AliAnalysisTaskHFEpidQA::UserCreateOutputObjects(){
   // Initialize PID QA
   //
   fOutput = new TList;
+  fOutput->SetOwner();
 
   // Counter for number of events
   fOutput->Add((fEvents = new TH1I("nEvents", "NumberOfEvents", 1, 1, 2)));
 
   fPIDqa = new AliHFEpidQA;
+  if(fTRDTotalChargeInSlice0) fPIDqa->SetTRDTotalChargeInSlice0();
   if(HasV0pidQA()) fPIDqa->SetV0pidQA();
   if(HasRecalculateTRDpid()) fPIDqa->SetRecalculateTRDpid();
   if(fNNref) fPIDqa->SetNNref(fNNref);
index 2b08443..dcb1df4 100644 (file)
@@ -48,6 +48,7 @@ class AliAnalysisTaskHFEpidQA : public AliAnalysisTaskSE{
     Bool_t HasRecalculateTRDpid() const { return TestBit(kRecalculateTRDpid); };
     void SetV0pidQA(Bool_t v0pidQA = kTRUE) { SetBit(kV0pidQA, v0pidQA); };
     void SetRecalculateTRDpid(Bool_t recal = kTRUE) { SetBit(kRecalculateTRDpid, recal); };
+    void SetTRDTotalChargeInSlice0() { fTRDTotalChargeInSlice0 = kTRUE; }
 
     void SetNNref(TFile *f) { fNNref = f; };
 
@@ -62,6 +63,7 @@ class AliAnalysisTaskHFEpidQA : public AliAnalysisTaskSE{
     TList *fOutput;         //! Container for output histos
     TH1 *fEvents;           //! Number of Events
     TFile  *fNNref;         //  reference file for NN
+    Bool_t fTRDTotalChargeInSlice0;   // Fix for Foreware/Backward compatibility
 
     ClassDef(AliAnalysisTaskHFEpidQA, 1)
 };
index ff6dfec..b35a361 100644 (file)
@@ -30,8 +30,6 @@
 #include "AliLog.h"
 #include "AliExternalTrackParam.h"
 
-#include "AliHFEcollection.h"
-
 #include "AliHFEV0cuts.h"
 
 ClassImp(AliHFEV0cuts)
@@ -46,6 +44,7 @@ AliHFEV0cuts::AliHFEV0cuts():
   , fCurrentV0id(0)
   , fPdaughterPDG(0)
   , fNdaughterPDG(0)
+  , fDestBits(0)
 {
  
   //
@@ -60,8 +59,9 @@ AliHFEV0cuts::~AliHFEV0cuts()
   //
   // destructor
   //
-  if (fQA) delete fQA;
-  if (fQAmc) delete fQAmc;
+  
+  if (fQA && TESTBIT(fDestBits, kBitQA)) delete fQA;
+  if (fQAmc && TESTBIT(fDestBits, kBitQAmc)) delete fQAmc;
 }
 
 //________________________________________________________________
@@ -75,6 +75,7 @@ AliHFEV0cuts::AliHFEV0cuts(const AliHFEV0cuts &ref):
   , fCurrentV0id(0)
   , fPdaughterPDG(0)
   , fNdaughterPDG(0)
+  , fDestBits(0)
 {
   //
   // Copy constructor
@@ -122,6 +123,12 @@ void AliHFEV0cuts::Init(const char* name){
   // [1] jus before the cut on given variable was applied, but after all the previous cuts
   //
 
+  memset(&fDestBits, 0, sizeof(UInt_t));
+  // now set the first two bits to 1
+  SETBIT(fDestBits, kBitQA);
+  SETBIT(fDestBits, kBitQAmc);
+  
+
   fQA = new AliHFEcollection("fQA", name);
 
   fQAmc = new AliHFEcollection("fQAmc", name);
@@ -409,9 +416,6 @@ Bool_t AliHFEV0cuts::GammaCuts(AliESDv0 *v0){
     }
   }
 
-  // loose cuts first
-  //if(LooseRejectK0(v0) || LooseRejectLambda(v0)) return kFALSE;
-
   AliVTrack* daughter[2];
   Int_t pIndex = 0, nIndex = 0;
   if(CheckSigns(v0)){
@@ -453,7 +457,7 @@ Bool_t AliHFEV0cuts::GammaCuts(AliESDv0 *v0){
   // possible new cuts
   //
   // separation cut at the entrance to the TPC
-  const Double_t cutSeparation = 0.;          // ORG 3.0 cm
+  const Double_t cutSeparation = 1.;          // ORG 3.0 cm
 
 
 
@@ -1024,9 +1028,7 @@ Bool_t AliHFEV0cuts::LambdaCuts(AliESDv0 *v0, Bool_t &isLambda ){
       fQAmc->Fill("h_Mass_L_as_K0",  v0->GetEffMass(2, 0));
     }
   }
-  // loose cuts first
-  //if(LooseRejectK0(v0) || LooseRejectGamma(v0)) return kFALSE;
-  
+
   const Double_t cL0mass=TDatabasePDG::Instance()->GetParticle(kLambda0)->Mass();  // PDG lambda mass
 
   AliVTrack* daughter[2];
@@ -1345,19 +1347,20 @@ Bool_t AliHFEV0cuts::LambdaCuts(AliESDv0 *v0, Bool_t &isLambda ){
     fQAmc->Fill("h_alambda_p_B", iP);
   }
   //
-  if(4 == TMath::Abs(fCurrentV0id)){
-    fQAmc->Fill("h_ProtonL_P_S", 5, p[ixMC[0]]);
-    fQAmc->Fill("h_PionL_P_S", 5, p[ixMC[1]]);
-  }
-  else if(-2 != fCurrentV0id){
-    fQAmc->Fill("h_ProtonL_P_B", 5, p[ixMC[0]]);
-    fQAmc->Fill("h_PionL_P_B", 5, p[ixMC[1]]);
+  if(fMCEvent){
+    if(4 == TMath::Abs(fCurrentV0id)){
+      fQAmc->Fill("h_ProtonL_P_S", 5, p[ixMC[0]]);
+      fQAmc->Fill("h_PionL_P_S", 5, p[ixMC[1]]);
+    }
+    else if(-2 != fCurrentV0id){
+      fQAmc->Fill("h_ProtonL_P_B", 5, p[ixMC[0]]);
+      fQAmc->Fill("h_PionL_P_B", 5, p[ixMC[1]]);
+    }
   }
-  
   return kTRUE;
 }
 //________________________________________________________________
-Double_t AliHFEV0cuts::OpenAngle(AliESDv0 *v0) const {
+Double_t AliHFEV0cuts::OpenAngle(AliESDv0 const *v0){
   //
   // Opening angle between two daughter tracks
   //
@@ -1374,7 +1377,7 @@ Double_t AliHFEV0cuts::OpenAngle(AliESDv0 *v0) const {
   return TMath::Abs(openAngle);
 }
 //________________________________________________________________
-Double_t AliHFEV0cuts::PsiPair(AliESDv0 *v0) {
+Double_t AliHFEV0cuts::PsiPair(const AliESDv0 *v0) {
   //
   // Angle between daughter momentum plane and plane 
   // 
@@ -1445,7 +1448,7 @@ Double_t AliHFEV0cuts::PsiPair(AliESDv0 *v0) {
   return psiPair; 
 }
 //________________________________________________________________
-AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec){
+AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack const *pdaughter, AliVTrack const *ndaughter, Int_t pspec, Int_t nspec){
   //
   // Creates a mother particle
   //
@@ -1483,42 +1486,8 @@ AliKFParticle *AliHFEV0cuts::CreateMotherParticle(AliVTrack* const pdaughter, Al
 
   return m;
 }
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectK0(AliESDv0 * const v0) const {
-  //
-  // Reject K0 based on loose cuts
-  //
-  Double_t mass = v0->GetEffMass(AliPID::kPion, AliPID::kPion);
-  if(mass > 0.494 && mass < 0.501) return kTRUE;
-  return kFALSE;
-}
-
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectLambda(AliESDv0 * const v0) const {
-  //
-  // Reject Lambda based on loose cuts
-  //
-  Double_t mass1 = v0->GetEffMass(AliPID::kPion, AliPID::kProton);
-  Double_t mass2 = v0->GetEffMass(AliPID::kProton, AliPID::kPion);
-  
-  if(mass1 > 1.1 && mass1 < 1.12) return kTRUE;
-  if(mass2 > 1.1 && mass2 < 1.12) return kTRUE;
-  return kFALSE;
-}
-
-//_________________________________________________
-Bool_t AliHFEV0cuts::LooseRejectGamma(AliESDv0 * const v0) const {
-  //
-  // Reject Lambda based on loose cuts
-  //
-  Double_t mass = v0->GetEffMass(AliPID::kElectron, AliPID::kElectron);
-  
-  if(mass < 0.02) return kTRUE;
-  return kFALSE;
-}
 //___________________________________________________________________
-void  AliHFEV0cuts::Armenteros(AliESDv0 *v0, Float_t val[2]){
+void  AliHFEV0cuts::Armenteros(const AliESDv0 *v0, Float_t val[2]){
   //
   // computes the Armenteros variables for given V0
   // fills the histogram
@@ -1555,7 +1524,7 @@ void  AliHFEV0cuts::Armenteros(AliESDv0 *v0, Float_t val[2]){
 
 }
 //___________________________________________________________________
-Bool_t AliHFEV0cuts::CheckSigns(AliESDv0* const v0){
+Bool_t AliHFEV0cuts::CheckSigns(AliESDv0 const *v0){
   //
   // check wheter the sign was correctly applied to 
   // V0 daughter tracks
index 8651d4d..48f2dcd 100644 (file)
@@ -34,6 +34,10 @@ class AliVTrack;
 
 class AliHFEV0cuts : public TObject {
  public:
+  enum{
+    kBitQA = 1,
+      kBitQAmc = 2
+      };
   enum{ // Reconstructed V0
     kUndef = 0,
       kRecoGamma = 1,
@@ -62,8 +66,14 @@ class AliHFEV0cuts : public TObject {
   void SetInputEvent(AliVEvent* const e)      { fInputEvent = e; };
   void SetPrimaryVertex(AliKFVertex* const v) { fPrimaryVertex = v; };
   
-  TList* GetList()    { return fQA->GetList(); };
-  TList* GetListMC()  { return fQAmc->GetList(); };
+  TList* GetList()    { 
+    CLRBIT(fDestBits, kBitQA);
+    return fQA->GetList(); 
+  };
+  TList* GetListMC()  { 
+    CLRBIT(fDestBits, kBitQAmc);
+    return fQAmc->GetList(); 
+  };
   
   Bool_t   TrackCutsCommon(AliESDtrack* track);
   Bool_t   V0CutsCommon(AliESDv0 *v0);
@@ -71,16 +81,12 @@ class AliHFEV0cuts : public TObject {
   Bool_t   K0Cuts(AliESDv0 *v0);
   Bool_t   LambdaCuts(AliESDv0 *v0, Bool_t &isLambda);
  
-  Bool_t LooseRejectK0(AliESDv0 * const v0) const;
-  Bool_t LooseRejectLambda(AliESDv0 * const v0) const;
-  Bool_t LooseRejectGamma(AliESDv0 * const v0) const;
-
-  void     Armenteros(AliESDv0 *v0, Float_t val[2]);
+  void     Armenteros(const AliESDv0 *v0, Float_t val[2]);
 
-  Double_t OpenAngle(AliESDv0 *v0) const;//opening angle between V0 daughters; close to zero for conversions
-  Double_t PsiPair(AliESDv0 *v0);
+  Double_t OpenAngle(AliESDv0 const *v0);//opening angle between V0 daughters; close to zero for conversions
+  Double_t PsiPair(const AliESDv0 *v0);
   
-  Bool_t   CheckSigns(AliESDv0* const v0);
+  Bool_t   CheckSigns(AliESDv0 const *v0);
   Bool_t   GetConvPosXY(AliESDtrack * const ptrack, AliESDtrack * const ntrack, Double_t convpos[2]);
   Bool_t   GetHelixCenter(AliESDtrack * const track, Double_t b,Int_t charge, Double_t center[2]);
 
@@ -89,7 +95,7 @@ class AliHFEV0cuts : public TObject {
   void     SetCurrentV0id(Int_t id) { fCurrentV0id = id; };
   void     SetDaughtersID(Int_t d[2]) {fPdaughterPDG = d[0]; fNdaughterPDG = d[1]; };
   
-  AliKFParticle *CreateMotherParticle(AliVTrack* const pdaughter, AliVTrack* const ndaughter, Int_t pspec, Int_t nspec);
+  AliKFParticle *CreateMotherParticle(AliVTrack const *pdaughter, AliVTrack const *ndaughter, Int_t pspec, Int_t nspec);
 
  private:
   void Copy(TObject &ref) const;
@@ -106,7 +112,8 @@ class AliHFEV0cuts : public TObject {
   Int_t                fPdaughterPDG;   // MC id of the positive daugeter
   Int_t                fNdaughterPDG;   // MC id of the negative daugeter
 
-
+  UInt_t               fDestBits;           // status bits for destructor
   ClassDef(AliHFEV0cuts, 1)
 };
     
index 8b2006d..d196ad2 100644 (file)
@@ -45,9 +45,8 @@
 #include "AliHFEV0pid.h"
 ClassImp(AliHFEV0pid)
 
-//____________________________________________________________
 AliHFEV0pid::AliHFEV0pid():
-  TObject()
+  TNamed()
   , fInputEvent(NULL)
   , fNtracks(0)
   , fMCEvent(NULL)
@@ -66,10 +65,39 @@ AliHFEV0pid::AliHFEV0pid():
   , fQA(NULL)
   , fV0cuts(NULL)
   , fOutput(NULL)
+  , fDestBits(0)
+
 {
   //
   // Default constructor
   //
+}
+//____________________________________________________________
+AliHFEV0pid::AliHFEV0pid(const char *name):
+  TNamed(name, "")
+  , fInputEvent(NULL)
+  , fNtracks(0)
+  , fMCEvent(NULL)
+  , fMCon(kFALSE)
+  , fPrimaryVertex(NULL)
+  , fElectrons(NULL)
+  , fPionsK0(NULL)
+  , fPionsL(NULL)
+  , fKaons(NULL)
+  , fProtons(NULL)
+  , fGammas(NULL)
+  , fK0s(NULL)
+  , fLambdas(NULL)
+  , fAntiLambdas(NULL)
+  , fIndices(NULL)
+  , fQA(NULL)
+  , fV0cuts(NULL)
+  , fOutput(NULL)
+  , fDestBits(0)
+{
+  //
+  // Standard constructor
+  //
   fElectrons = new TObjArray();
   fPionsK0 = new TObjArray();
   fPionsL = new TObjArray();
@@ -90,8 +118,8 @@ AliHFEV0pid::AliHFEV0pid():
   fIndices = new AliHFEV0pidTrackIndex();
   
 }
-
 //____________________________________________________________
+
 AliHFEV0pid::~AliHFEV0pid(){
   //
   // Destructor
@@ -109,9 +137,12 @@ AliHFEV0pid::~AliHFEV0pid(){
   if(fAntiLambdas) delete fAntiLambdas;
 
   if(fIndices) delete fIndices;
-  if(fQA) delete fQA;
   if(fV0cuts) delete fV0cuts;
-  if(fOutput) delete fOutput;
+
+  if(TESTBIT(fDestBits, 1)){
+    if(fQA) delete fQA;
+    if(fOutput) delete fOutput;
+  }
 }
 
 //____________________________________________________________
@@ -120,7 +151,11 @@ void AliHFEV0pid::InitQA(){
   // Initialize QA histograms
   //
   
+  memset(&fDestBits, 0, sizeof(fDestBits));
+  SETBIT(fDestBits, 1);
+
   fOutput = new TList();
+  fOutput->SetOwner();
 
   fV0cuts = new AliHFEV0cuts();
   fV0cuts->Init("V0cuts");
@@ -798,7 +833,7 @@ AliHFEV0pid::AliHFEV0pidTrackIndex::~AliHFEV0pidTrackIndex(){
 }
 
 //____________________________________________________________
-void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
+const void AliHFEV0pid::AliHFEV0pidTrackIndex::Flush(){
   //
   // Reset containers
   //
@@ -907,6 +942,8 @@ TList *AliHFEV0pid::GetListOfQAhistograms(){
   // Getter for V0 PID QA histograms
   //
 
+  CLRBIT(fDestBits, 1);
+
   TList *tmp = fV0cuts->GetList();
   tmp->SetName("V0cuts");
   fOutput->Add(tmp);
index a56638a..755015d 100644 (file)
@@ -20,8 +20,8 @@
 #ifndef ALIHFEV0PID_H
 #define ALIHFEV0PID_H
 
-#ifndef ROOT_TObject
-#include <TObject.h>
+#ifndef ROOT_TNamed
+#include <TNamed.h>
 #endif
 
 class TObjArray;
@@ -39,9 +39,10 @@ class AliMCEvent;
 class AliHFEV0cuts;
 class AliHFEcollection;
 
-class AliHFEV0pid : public TObject{
+class AliHFEV0pid : public TNamed{
   public:
     AliHFEV0pid();
+    AliHFEV0pid(const char *name);
     ~AliHFEV0pid();
 
     void  Process(AliVEvent * const inputEvent);
@@ -96,7 +97,7 @@ class AliHFEV0pid : public TObject{
       Int_t GetNumberOfPionsL() const { return fNPionsL; };
       Int_t GetNumberOfKaons() const { return fNKaons; };
       Int_t GetNumberOfProtons() const { return fNProtons; };
-      void Flush();
+      const void Flush();
       
     private:
       AliHFEV0pidTrackIndex(const AliHFEV0pidTrackIndex &ref);
@@ -136,6 +137,8 @@ class AliHFEV0pid : public TObject{
     AliHFEV0cuts     *fV0cuts;       // separate class for studying and applying the V0 cuts
     TList       *fOutput;            // collection list
 
+    UInt_t       fDestBits;              // logical bits for destructor
+
     ClassDef(AliHFEV0pid, 1)          // V0 PID Class
 
 };
index 51f7f11..1e9c083 100644 (file)
@@ -37,7 +37,7 @@ ClassImp(AliHFEV0pidMC)
   AliHFEV0pidMC::AliHFEV0pidMC():
     fMC(0x0)
     , fColl(NULL)
-    , fV0cuts(NULL)
+    , fDestBits(0)
 {
   //
   // default constructor
@@ -50,13 +50,15 @@ AliHFEV0pidMC::~AliHFEV0pidMC(){
   // destructor
   //
   if (fColl) delete fColl;
-  //if (fV0cuts) delete fV0cuts;
 }
 //____________________________________________________________
 void AliHFEV0pidMC::Init(){
   //
   // initialize objects
   //
+  
+  memset(&fDestBits, 0, sizeof(UInt_t));
+  SETBIT(fDestBits, 1);
 
   fColl = new AliHFEcollection("V0pidMC", "MC based V0 benchmarking");
   // QA
@@ -218,6 +220,8 @@ TList *AliHFEV0pidMC::GetListOfQAhistograms(){
   //
   // Get QA histograms
   //
+  
+  CLRBIT(fDestBits, 1);
   if(fColl)
     return fColl->GetList();
   return NULL;
index 77c2adc..260f8eb 100644 (file)
@@ -29,7 +29,6 @@ class TList;
 
 class AliMCEvent;
 
-class AliHFEV0cuts;
 class AliHFEcollection;
 
 class AliHFEV0pidMC : public TObject {
@@ -54,7 +53,8 @@ class AliHFEV0pidMC : public TObject {
 
   AliMCEvent*         fMC;      // MC event
   AliHFEcollection*   fColl;    // Histogram collection
-  AliHFEV0cuts*       fV0cuts;  // V0 cut class
+
+  UInt_t              fDestBits;    // status bits for the destructor
 
    ClassDef(AliHFEV0pidMC, 1)   // QA class for V0 PID
 };
index 0a9050f..e0af341 100644 (file)
@@ -20,8 +20,6 @@
 // Matus Kalisky  <matus.kalisky@cern.ch>
 //
 
-//#include <iostream>
-
 #include <TH1F.h>
 #include <TH2F.h>
 #include <THnSparse.h>
@@ -31,6 +29,7 @@
 #include <TMath.h>
 
 #include "AliLog.h"
+
 #include "AliHFEcollection.h"
 
 using namespace std;
@@ -99,6 +98,7 @@ void AliHFEcollection::Copy(TObject &ref) const {
 
   // Clone List Content
   target.fList = new THashList();          
+  target.fList->SetOwner();
   for(Int_t ien = 0; ien < fList->GetEntries(); ien++)
     target.fList->Add(fList->At(ien)->Clone());
 }
@@ -563,3 +563,21 @@ void AliHFEcollection::Browse(TBrowser *b)
       }
    }
 }
+//____________________________________________________________________
+void AliHFEcollection::Print(Option_t *) const{
+  //
+  // Print content of the collection
+  //
+  TIter histIter(fList);
+  TObject *o = NULL;
+  Int_t nHistos = 0;
+  printf("Collection %s\n", GetName());
+  printf("Content of the collection:\n=========================================\n");
+  while((o = histIter())){
+    printf("Histo %s, Type %s\n", o->GetName(), o->IsA()->GetName());
+    nHistos++;
+  }
+  printf("Number of histos in the collection: %d\n", nHistos);
+  printf("\n");
+}
+
index d1e3150..307d3eb 100644 (file)
@@ -67,6 +67,7 @@ class AliHFEcollection : public TNamed{
     
 
   Long64_t Merge(TCollection *list);
+  virtual void Print(Option_t *) const;
 
   // Get functions
   TList* GetList() const { return fList; }
index 63ee2d2..beb9513 100644 (file)
@@ -366,6 +366,11 @@ void AliHFEcuts::Initialize(){
   SetHFElectronTRDCuts();
   SetHFElectronDcaCuts();
 
+  // Connect the event cuts
+  SetEventCutList(kEventStepGenerated);
+  SetEventCutList(kEventStepReconstructed);
+
+
 }
 
 //__________________________________________________________________
@@ -418,9 +423,9 @@ void AliHFEcuts::SetEventCutList(Int_t istep){
 void AliHFEcuts::SetParticleGenCutList(){
   //
   // Initialize Particle Cuts for Monte Carlo Tracks
-  // Production Vertex: < 1cm (Beampipe)
+  // Production Vertex Radius: < 3cm
   // Particle Species: Electrons
-  // Eta: < 0.9 (TRD-TOF acceptance)
+  // Eta: < 0.8
   //
   
   TObjArray *mcCuts = new TObjArray;
@@ -649,7 +654,7 @@ Bool_t AliHFEcuts::CheckParticleCuts(UInt_t step, TObject *o){
   // Checks the cuts without using the correction framework manager
   // 
   AliDebug(2, "Called\n");
-  TString stepnames[kNcutStepsMCTrack + kNcutStepsRecTrack + kNcutStepsDETrack + kNcutStepsSecvtxTrack + 1] = {"fPartGenCuts","fPartEvCutPileupZ","fPartEvCut","fPartAccCuts","fPartRecNoCuts","fPartRecKineITSTPCCuts", "fPartPrimCuts", "fPartHFECutsITS","fPartHFECutsTRD","fPartHFECutsDca", "fPartHFECutsSecvtx"};
+  TString stepnames[kNcutStepsMCTrack + kNcutStepsRecTrack + kNcutStepsDETrack + kNcutStepsSecvtxTrack + 1] = {"fPartGenCuts","fPartEvCutPileupZ","fPartEvCut","fPartAccCuts","fPartRecNoCuts","fPartRecKineITSTPCCuts", "fPartPrimCuts", "fPartHFECutsITS","fPartHFECutsTOF","fPartHFECutsTRD","fPartHFECutsDca", "fPartHFECutsSecvtx"};
   AliDebug(2, Form("Doing cut %s", stepnames[step].Data()));
  TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(stepnames[step].Data()));
   if(!cuts) return kTRUE;
@@ -661,3 +666,21 @@ Bool_t AliHFEcuts::CheckParticleCuts(UInt_t step, TObject *o){
   }
   return status;
 }
+
+
+//__________________________________________________________________
+Bool_t AliHFEcuts::CheckEventCuts(const char*namestep, TObject *o){
+  //
+  // Checks the cuts without using the correction framework manager
+  // 
+  AliDebug(2, "Called\n");
+  TObjArray *cuts = dynamic_cast<TObjArray *>(fCutList->FindObject(namestep));
+  if(!cuts) return kTRUE;
+  TIter it(cuts);
+  AliCFCutBase *mycut;
+  Bool_t status = kTRUE;
+  while((mycut = dynamic_cast<AliCFCutBase *>(it()))){
+    status &= mycut->IsSelected(o);
+  }
+  return status;
+}
index c464eec..cd35792 100644 (file)
@@ -82,6 +82,7 @@ class AliHFEcuts : public TNamed{
     void Initialize();
 
     Bool_t CheckParticleCuts(UInt_t step, TObject *o);
+    Bool_t CheckEventCuts(const char*namestep, TObject *o);
   
     TList *GetQAhistograms() const { return fHistQA; }
     
index 82334b9..0cb061e 100644 (file)
 #include <THnSparse.h>
 #include <TString.h>
 
+#include <TMath.h>
+#include "AliESDInputHandler.h"
+//#include "AliVCluster.h"
+//#include "AliVCaloCells.h"
+//#include "AliVEvent.h"
+#include "AliMagF.h"
+
 #include "AliLog.h"
 #include "AliPID.h"
 #include "AliVParticle.h"
-#include "AliVTrack.h"
-#include "AliESDtrack.h"
+//#include "AliVTrack.h"
+//#include "AliESDtrack.h"
 #include "AliHFEcollection.h"
 #include "AliHFEpid.h"
 #include "AliHFEpidBase.h"
@@ -45,7 +52,7 @@
 #include "AliHFEtools.h"
 #include "AliHFEemcalPIDqa.h"
 
-ClassImp(AliHFEpidEMCAL)
+ClassImp(AliHFEemcalPIDqa)
 
 //_________________________________________________________
 AliHFEemcalPIDqa::AliHFEemcalPIDqa():
@@ -149,30 +156,59 @@ void AliHFEemcalPIDqa::Initialize(){
   const Double_t kTPCSigMax = 140.;
 
   // 1st histogram: TPC dEdx with/without EMCAL (p, pT, TPC Signal, Centrality)
-  Int_t nBins[4] = {20, 20, 400, kCentralityBins};
-  Double_t min[4] = {kMinP, kMinP, kTPCSigMim,  0};
-  Double_t max[4] = {kMaxP, kMaxP, kTPCSigMax,  11.};
-  fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality", 4, nBins, min, max);
-  
+  Int_t nBins[6] = {AliPID::kSPECIES + 1, 20, 20, 400, kCentralityBins, 2};
+  Double_t min[6] = {-1, kMinP, kMinP, kTPCSigMim,  0, 0.};
+  Double_t max[6] = {AliPID::kSPECIES, kMaxP, kMaxP, kTPCSigMax,  11., 2.};
+  fHistos->CreateTHnSparse("EMCAL_TPCdedx", "EMCAL signal; species; p [GeV/c]; pT [GeV/c] ; TPC signal [a.u.]; Centrality; PID Step", 6, nBins, min, max);
+
+  //2nd histogram: EMCAL signal - E/p and Rmatch
+  Int_t nBins2[6] = {AliPID::kSPECIES + 1, 40, 40, 1000, 250, 2};
+  Double_t min2[6] = {-1, kMinP, kMinP, 0,  0, 0.};
+  Double_t max2[6] = {AliPID::kSPECIES, kMaxP, kMaxP, 10,  1., 2.};
+  fHistos->CreateTHnSparse("EMCAL_Signal", "EMCAL true signal; species; p [GeV/c]; pT [GeV/c] ; E/p; Rmatch; PID Step", 6, nBins2, min2, max2);
+    
 }
 
+
+
+
 //_________________________________________________________
-void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t /*step*/){
+void AliHFEemcalPIDqa::ProcessTrack(const AliHFEpidObject *track,AliHFEdetPIDqa::EStep_t step){
   //
   // Fill TPC histograms
   //
   //AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
   Float_t centrality = track->GetCentrality();
 
+  //const AliVTrack *vtrack = dynamic_cast<const AliVTrack *>(track->GetRecTrack());
+  //const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(vtrack);
   const AliESDtrack *esdtrack = static_cast<const AliESDtrack *>(track->GetRecTrack());
 
-  Double_t contentSignal[4];
-  contentSignal[0] = track->GetRecTrack()->P();
-  contentSignal[1] = track->GetRecTrack()->Pt();
-  contentSignal[2] = esdtrack->GetTPCsignal(); //?
-  contentSignal[3] = centrality;
+  Int_t species = track->GetAbInitioPID();
+  if(species >= AliPID::kSPECIES) species = -1;
+
+  Double_t contentSignal[6];
+  contentSignal[0] = species;
+  contentSignal[1] = track->GetRecTrack()->P();
+  contentSignal[2] = track->GetRecTrack()->Pt();
+  contentSignal[3] = esdtrack->GetTPCsignal(); //?
+  contentSignal[4] = centrality;
+  contentSignal[5] = step == AliHFEdetPIDqa::kBeforePID ? 0. : 1.;
+
+  TVector3 emcsignal = MomentumEnergyMatchV2(esdtrack);
+
+
+  Double_t contentSignal2[6];
+  contentSignal2[0] = species;
+  contentSignal2[1] = track->GetRecTrack()->P();
+  contentSignal2[2] = track->GetRecTrack()->Pt();
+  contentSignal2[3] = emcsignal(0);
+  contentSignal2[4] = emcsignal(1);
+  contentSignal2[5] = contentSignal[5];
+
   //printf("ProcessTrack ; Print Content %g; %g; %g; %g \n",contentSignal[0],contentSignal[1],contentSignal[2],contentSignal[3]); 
   fHistos->Fill("EMCAL_TPCdedx", contentSignal);
+  fHistos->Fill("EMCAL_Signal", contentSignal2);
 }
 
 //_________________________________________________________
@@ -186,3 +222,116 @@ TH1 *AliHFEemcalPIDqa::GetHistogram(const char *name){
   return dynamic_cast<TH1 *>(histo);
 }
 
+
+//___________________________________________________________________________
+TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const
+{ // Returns e/p & Rmatch
+
+  Float_t  clsPos[3];
+  Double_t trkPos[3];
+  Double_t matchclsE = 9999.9;
+  TVector3 refVec(-9999,-9999,-9999);
+
+  AliESDEvent *evt = esdtrack->GetESDEvent();
+
+   //Int_t icl = esdtrack->GetEMCALcluster();
+   Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();
+
+   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+   if(!cluster->IsEMCAL()) return refVec;
+
+   cluster->GetPosition(clsPos);
+   esdtrack->GetXYZ(trkPos);
+
+   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+
+   Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
+   Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
+
+   double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
+
+   matchclsE = cluster->E();
+
+   //double feop = -9999.9;
+   //if(matchclsE<9999) 
+   double feop = matchclsE/esdtrack->P();
+
+   //   if(feop!=-9999.9)printf("%f\n",feop) ; 
+   TVector3 emcsignal(feop,rmatch,0);
+
+   return emcsignal;
+
+}
+
+
+
+
+/*
+//___________________________________________________________________________
+TVector3 AliHFEemcalPIDqa::MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const
+{ // Returns e/p & Rmatch
+
+  Float_t  clsPos[3];
+  Double_t trkPos[3];
+  Double_t Rmatch;
+  Double_t matchclsE = 9999.9;
+  TVector3 refVec(-9999,-9999,-9999);
+
+  AliESDEvent *evt = esdtrack->GetESDEvent();
+  Double_t  magF = evt->GetMagneticField();
+  Double_t magSign = 1.0;
+  if(magF<0)magSign = -1.0;
+  //printf("magF ; %g ; %g \n", magF,magSign);
+
+  if (!TGeoGlobalMagField::Instance()->GetField()) {
+          printf("Loading field map...\n");
+          //AliMagF* field = new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG);
+          AliMagF* field = new AliMagF("Maps","Maps", magSign, magSign, AliMagF::k5kG); // for 10d
+          TGeoGlobalMagField::Instance()->SetField(field);
+                    }
+
+  AliEMCALTrack *emctrack = new AliEMCALTrack(*esdtrack);
+  Double_t fieldB[3];
+  emctrack->GetBxByBz(fieldB);
+  //printf("%g %g %g \n", fieldB[0], fieldB[1], fieldB[2]);
+
+  for(Int_t icl=0; icl<evt->GetNumberOfCaloClusters(); icl++)
+   {
+
+   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+   if(!cluster->IsEMCAL()) return refVec;
+
+   cluster->GetPosition(clsPos);
+   if(!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) )  return refVec;
+   emctrack->GetXYZ(trkPos);
+
+   TVector3 clsPosVec(clsPos[0],clsPos[1],clsPos[2]);
+   TVector3 trkPosVec(trkPos[0],trkPos[1],trkPos[2]);
+
+   Double_t delEmcphi = clsPosVec.Phi()-trkPosVec.Phi();  // track cluster matching
+   Double_t delEmceta = clsPosVec.Eta()-trkPosVec.Eta();  // track cluster matching
+
+   double rmatch = sqrt(pow(delEmcphi,2)+pow(delEmceta,2));
+
+   if(rmatch<0.02)
+      {
+       matchclsE = cluster->E();
+       Rmatch = rmatch;
+      }
+  }
+   delete emctrack;
+
+   //double feop = -9999.9;
+   //if(matchclsE<9999) 
+   double feop = matchclsE/esdtrack->P();
+
+   //   if(feop!=-9999.9)printf("%f\n",feop) ; 
+   TVector3 emcsignal(feop, Rmatch, 0);
+   return emcsignal;
+
+}
+*/
+
index 1b3df4e..69cd74a 100644 (file)
@@ -31,6 +31,9 @@ class AliHFEpidObject;
 class AliESDtrack;
 class AliAODTrack;
 
+class TVector3;
+
+
 class AliHFEemcalPIDqa : public AliHFEdetPIDqa{
   public:
     AliHFEemcalPIDqa();
@@ -49,6 +52,10 @@ class AliHFEemcalPIDqa : public AliHFEdetPIDqa{
     TH1 *GetHistogram(const char *name); 
     AliHFEcollection *GetHistoCollection() const { return fHistos; }
 
+    TVector3 MomentumEnergyMatchV2(const AliESDtrack *esdtrack) const;
+    //TVector3 MomentumEnergyMatchV1(const AliESDtrack *esdtrack) const;
+
+
   protected:
     void ProcessESDtrack(const AliESDtrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species);
     void ProcessAODtrack(const AliAODTrack *track, AliHFEdetPIDqa::EStep_t step, Int_t species);
index 8dc68e0..0bf37e8 100644 (file)
@@ -300,6 +300,16 @@ Bool_t AliHFEextraCuts::CheckRecCuts(AliVTrack *track){
                if(fCheck && !(statusL0 || statusL1))
                    SETBIT(survivedCut, kPixelITS);
                    break;
+      case kExclusiveSecond:
+        AliDebug(2, "Exlusive second");
+        if(fCheck){ // Cut out tracks which pass a dead ITS Layer 0
+          if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0) && statusL0)
+            SETBIT(survivedCut, kPixelITS);
+        } else {
+          if(TESTBIT(itsPixel,1) && !TESTBIT(itsPixel,0))
+            SETBIT(survivedCut, kPixelITS);
+        }
+        break;
       default: 
         AliDebug(2, "None");
         break;
index c556026..b586bea 100644 (file)
@@ -37,7 +37,8 @@ class AliHFEextraCuts : public AliCFCutBase{
       kSecond = 1,
       kBoth = 2,
       kNone = 3,
-      kAny = 4
+      kAny = 4,
+      kExclusiveSecond = 5
     } ITSPixel_t;
     typedef enum{
       kFound = 0,
index 14c3c0b..8b56f7f 100644 (file)
@@ -124,6 +124,9 @@ void AliHFEmcQA::PostAnalyze() const
 //_______________________________________________________________________________________________
 void AliHFEmcQA::SetBackgroundWeightFactor(Double_t *elecBackgroundFactor, Double_t *binLimit)
 {
+   //
+   // copy background weighting factors into data member
+   //
    memcpy(fElecBackgroundFactor,elecBackgroundFactor,sizeof(Double_t) * kElecBgSpecies * kBgPtBins);
    memcpy(fBinLimit,binLimit,sizeof(Double_t) * (kBgPtBins+1));
 /*
index 006d159..045eab0 100644 (file)
@@ -224,7 +224,7 @@ Bool_t AliHFEpid::InitializePID(){
 }
 
 //____________________________________________________________
-Bool_t AliHFEpid::IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
+Bool_t AliHFEpid::IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont, const Char_t *contname, AliHFEpidQAmanager *pidqa){
   //
   // Select Tracks
   //
@@ -307,78 +307,77 @@ void AliHFEpid::ConfigureTPCrejectionSimple(){
 }
 
 //____________________________________________________________
-void AliHFEpid::ConfigureTPCrejection(const char *lowerCutParam, Double_t * const params, Float_t upperTPCCut, Float_t TOFCut){
+void AliHFEpid::ConfigureTOF(Float_t TOFCut){
   //
-  // Combined TPC-TOF PID
-  // if no function parameterizaion is given, then the default one (exponential) is chosen
+  // Set Number of sigmas for TOF PID
   //
-  if(HasMCData()) AliInfo("Configuring TPC for MC\n");
-  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   AliHFEpidTOF *tofpid = dynamic_cast<AliHFEpidTOF *>(fDetectorPID[kTOFpid]);
   if(tofpid) tofpid->SetTOFnSigma(TOFCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure centrality dependent cut function for TPC PID 
+  //
+  ConfigureTPCcut(centralityBin, lowerCutParam, params, upperTPCCut);
+}
 
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCdefaultCut(const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure default cut function for TPC PID 
+  //
+  ConfigureTPCcut(-1, lowerCutParam, params, upperTPCCut);
+}
+
+//____________________________________________________________
+void AliHFEpid::ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut){
+  //
+  // Cofigure cut function for TPC PID 
+  // if no function parameterizaion is given, then the default one (exponential) is chosen
+  //
+  
+  if(HasMCData()) AliInfo("Configuring TPC for MC\n");
+  AliHFEpidTPC *tpcpid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
   //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", lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
+  TF1 *upperCut = new TF1(Form("upperCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), "[0]", 0, 20); // Use constant upper cut
+  TF1 *lowerCut = new TF1(Form("lowerCut%s", centralityBin <  0 ? "Default" : Form("Bin%d", centralityBin)), lowerCutParam == NULL ? "[0] * TMath::Exp([1]*x) + [2]": lowerCutParam, 0, 20);
 
   upperCut->SetParameter(0, upperTPCCut); // pp
 
-//  upperCut->SetParameter(0, 3.5); //PbPb
-//  printf("upper %f \n",upperCut->GetParameter(0));
-  //upperCut->SetParameter(0, 2.7);
-  //upperCut->SetParameter(1, -0.4357);
   if(params){
-    for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++) lowerCut->SetParameter(ipar, params[ipar]);
+      for(Int_t ipar = 0; ipar < lowerCut->GetNpar(); ipar++)
+      {
+         lowerCut->SetParameter(ipar, params[ipar]);
+        //  printf("printout %i %s %f \n", centralityBin, lowerCutParam, params[ipar]);
+      }
   } else {
     // Set default parameterization
-      if(HasMCData()) lowerCut->SetParameter(0, -2.5);
-      else lowerCut->SetParameter(0, -4.03);  //pp
-//      else lowerCut->SetParameter(0, -3.83);  //PbPb
-  
-     
-    //else lowerCut->SetParameter(0, -3.71769);
-    //else lowerCut->SetParameter(0, -3.7);
-
+    if(HasMCData()) lowerCut->SetParameter(0, -2.5);
+    else lowerCut->SetParameter(0, -4.03);  //pp
     lowerCut->SetParameter(1, -0.22); // pp
-//     lowerCut->SetParameter(1,-0.36 ); // PbPb
-    //lowerCut->SetParameter(1, -0.40263);
-    //lowerCut->SetParameter(1, -0.8);
     
     if(HasMCData()) lowerCut->SetParameter(2, -2.2);
     else lowerCut->SetParameter(2, 0.92); //pp
-//     else lowerCut->SetParameter(2, 0.27); //PbPb
-   // else lowerCut->SetParameter(2, 0.92); //pp
-//     printf("lower0 %f \n",lowerCut->GetParameter(0));
-//     printf("lower1 %f \n",lowerCut->GetParameter(1));
-//     printf("lower2 %f \n",lowerCut->GetParameter(2));
-    //else lowerCut->SetParameter(2, 0.267857);
-    //else lowerCut->SetParameter(2, -0.35);
   }
 
 
   if(tpcpid){
     tpcpid->SetTPCnSigma(2);
-    tpcpid->SetUpperSigmaCut(upperCut);
-    tpcpid->SetLowerSigmaCut(lowerCut);
+    if(centralityBin < 0){
+      tpcpid->SetUpperSigmaCutDefault(upperCut);
+      tpcpid->SetLowerSigmaCutDefault(lowerCut);
+    } else {
+      tpcpid->SetUpperSigmaCutCentrality(upperCut, centralityBin);
+      tpcpid->SetLowerSigmaCutCentrality(lowerCut, centralityBin);
+    }
   }
   AddCommonObject(upperCut);
   AddCommonObject(lowerCut);
 }
 
 //____________________________________________________________
-void AliHFEpid::ConfigureTPCstrategyParis(){
-  //
-  // TPC alone, symmetric 3 sigma cut and 2 - -100 sigma pion rejection
-  //   
-  AliHFEpidTPC *pid = dynamic_cast<AliHFEpidTPC *>(fDetectorPID[kTPCpid]);
-  if(pid){
-    pid->SetTPCnSigma(2);
-    pid->SetRejectParticle(AliPID::kProton, 0., -3., 10., 3.);
-    pid->SetRejectParticle(AliPID::kKaon, 0., -3., 10., 3.);
-  }
-}
-
-//____________________________________________________________
 void AliHFEpid::PrintStatus() const {
   //
   // Print the PID configuration
index 1edfeb2..b290954 100644 (file)
@@ -62,7 +62,7 @@ class AliHFEpid : public TNamed{
     ~AliHFEpid();
     
     Bool_t InitializePID();
-    Bool_t IsSelected(AliHFEpidObject *track, AliHFEcontainer *cont = NULL, const Char_t *contname = "trackContainer", AliHFEpidQAmanager *qa = NULL);
+    Bool_t IsSelected(const AliHFEpidObject * const track, AliHFEcontainer *cont = NULL, const Char_t *contname = "trackContainer", AliHFEpidQAmanager *qa = NULL);
 
     Bool_t HasMCData() const { return TestBit(kHasMCData); };
 
@@ -82,15 +82,15 @@ class AliHFEpid : public TNamed{
       else return fgkDetectorName[kNdetectorPID];
     }    
     //-----Configure PID detectors with predefined stettings------
+    void ConfigureTOF(Float_t TOFcut = 3.);
     void ConfigureTPCasymmetric(Double_t pmin = 0.1, Double_t pmax = 20., Double_t sigmamin = -0.2, Double_t sigmamax = 5.);
     void ConfigureTPCrejectionSimple();
-    void ConfigureTPCrejection(const char *lowerCutParam = NULL, Double_t * const params = NULL, Float_t upperTPCCut=3.0, Float_t TOFCut=3.0);
-    void ConfigureTPCstrategyParis();
+    void ConfigureTPCcentralityCut(Int_t centralityBin, const char *lowerCutParam = NULL, const Double_t * const params = NULL, Float_t upperTPCCut=3.0);
+    void ConfigureTPCdefaultCut(const char *lowerCutParam = NULL, const Double_t * const params = NULL, Float_t upperTPCCut=3.0);
     //------------------------------------------------------------
 
   protected:
     Bool_t MakePidTpcTof(AliHFEpidObject *track);
-
   private:
     enum{
       kHasMCData = BIT(14)
@@ -116,6 +116,7 @@ class AliHFEpid : public TNamed{
       return det < kNdetectorPID ? TESTBIT(fEnabledDetectors, det): kFALSE;
     }
     //--------------------------------------------------
+    void ConfigureTPCcut(Int_t centralityBin, const char *lowerCutParam, const Double_t * const params, Float_t upperTPCCut);
 
     static const Char_t *fgkDetectorName[kNdetectorPID + 1]; // PID Detector Names
     AliHFEpidBase *fDetectorPID[kNdetectorPID];     //   Detector PID classes
index 27c36d0..a74386b 100644 (file)
@@ -40,14 +40,16 @@ class AliHFEpidObject{
       fkRecTrack(NULL), 
       fAnalysisType(kESDanalysis),
       fAbInitioPID(-1),
-      fCentrality(99.)
+      fCentrality(99),
+      fIsPbPb(kFALSE)         // Default: pp
       {
       }
     AliHFEpidObject(const AliHFEpidObject &ref):
       fkRecTrack(ref.fkRecTrack), 
       fAnalysisType(ref.fAnalysisType),
       fAbInitioPID(ref.fAbInitioPID),
-      fCentrality(ref.fCentrality)
+      fCentrality(ref.fCentrality),
+      fIsPbPb(ref.fIsPbPb)
       {
       }
     AliHFEpidObject &operator=(const AliHFEpidObject &ref);
@@ -57,19 +59,23 @@ class AliHFEpidObject{
     void SetMCTrack(const AliVParticle * mcTrack);
     void SetAnalysisType(AnalysisType_t type) { fAnalysisType = type; }
     void SetAbInitioPID(Int_t abInitioPID) { fAbInitioPID = abInitioPID; }
-    void SetCentrality(Float_t centrality) { fCentrality = centrality; }
+    void SetCentrality(Int_t centrality) { fCentrality = centrality; }
+    void SetPbPb() { fIsPbPb = kTRUE; }
+    void SetPP() { fIsPbPb = kFALSE; }
 
     const AliVParticle *GetRecTrack() const { return fkRecTrack; }
     Int_t GetAbInitioPID() const { return fAbInitioPID; }
-    Float_t GetCentrality() const { return fCentrality; }
+    Int_t GetCentrality() const { return fCentrality; }
     Bool_t IsAODanalysis() const { return fAnalysisType == static_cast<UChar_t>(kAODanalysis); }
     Bool_t IsESDanalysis() const { return fAnalysisType == static_cast<UChar_t>(kESDanalysis); }
+    Bool_t IsPbPb() const { return fIsPbPb; }
 
   private:
-    const AliVParticle *fkRecTrack;    // Reconstructed track
-    UChar_t fAnalysisType;      // Analysis Mode (ESD or AOD)
-    Int_t fAbInitioPID;         // AbInitio PID
-    Float_t fCentrality;        // Centrality Information
+    const AliVParticle *fkRecTrack;     // Reconstructed track
+    UChar_t fAnalysisType;              // Analysis Mode (ESD or AOD)
+    Int_t fAbInitioPID;                 // AbInitio PID
+    Int_t fCentrality;                  // Centrality Information
+    Bool_t fIsPbPb;                     // Collision type
 };
 
 class AliHFEpidBase : public TNamed{
index 48bbfc8..61038e6 100644 (file)
 #include "AliESDpid.h"
 
 #include "AliHFEdetPIDqa.h"
-//#include "AliHFEpidEMCAL.h"
+#include "AliHFEpidEMCAL.h"
 #include "AliHFEpidQAmanager.h"
 
-//#include "AliEMCALRecoUtils.h"
-//#include "AliEMCALGeometry.h"
-#include "AliVCluster.h"
-#include "AliVCaloCells.h"
-#include "AliVEvent.h"
+#include "AliHFEemcalPIDqa.h"
+//#include "AliVCluster.h"
+//#include "AliVCaloCells.h"
+//#include "AliVEvent.h"
 #include "AliLog.h"
-#include "AliEMCALPIDUtils.h"
 #include "AliPID.h"
-#include "AliESDEvent.h"
-#include "AliESDtrack.h"
-//#include "AliEMCALTrack.h"
-//#include "AliEMCALTracker.h"
+//#include "AliESDEvent.h"
+//#include "AliESDtrack.h"
 #include "AliHFEpidEMCAL.h"
 
 ClassImp(AliHFEpidEMCAL)
@@ -127,7 +123,7 @@ Bool_t AliHFEpidEMCAL::InitializePID(){
 }
 
 //___________________________________________________________________
-Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager */*pidqa*/) const
+Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *pidqa) const
 { //Function to a return a code indicating whether or not an electron candidate is selected
   //
   //
@@ -135,32 +131,62 @@ Int_t AliHFEpidEMCAL::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanage
   AliDebug(2, "PID object available");
   // EMCal not fESDpid  (s.s Feb. 11)
 
-  /*
-  AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
-  if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALpid)) return 0;
+  
+  //AliHFEpidObject::AnalysisType_t anaType = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
+  if(!((dynamic_cast<const AliVTrack *>(track->GetRecTrack()))->GetStatus() & AliESDtrack::kEMCALmatch)) return 0;
   AliDebug(2, "Track Has EMCAL PID");
-  */
-  //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
+  
+  //if(pidqa) 
+  pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kBeforePID);
   // not QA for EMCal, will be added (s.s Feb. 11)
 
-  Double_t eop = 0.;//MomentumEnergyMatch(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
+  Double_t eop = MomentumEnergyMatchV2(track->GetRecTrack()); // get eop (What is GetRecTrack ?)
   AliDebug(2, Form("Energy - Momentum Matching e/p : %f", eop));
   Int_t pdg = 0;
   if(eop>feopMim && eop<feopMax){
     pdg = 11;
-    //if(pidqa) pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
+    //if(pidqa) 
+    pidqa->ProcessTrack(track, AliHFEpid::kEMCALpid, AliHFEdetPIDqa::kAfterPID);
   }
   else
   {
    pdg = 211; // EMCal doesn't separate pi,k.p by e/p. return pion code as same as TRD
   } 
 
-    printf("eID %g ; %d \n",eop,pdg);  
+    AliDebug(1, Form("eID %g ; %d \n",eop,pdg));  
 
   return pdg;
 }
+
+
+//___________________________________________________________________________
+Double_t AliHFEpidEMCAL::MomentumEnergyMatchV2(const AliVParticle *const track) const
+{ // Returns e/p & Rmatch
+
+  Double_t matchclsE = 9999.9;
+  double feop = -9999.9;
+
+  const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
+  AliESDEvent *evt = esdtrack->GetESDEvent();
+
+   Int_t icl = (const_cast<AliESDtrack *>(esdtrack))->GetEMCALcluster();
+
+   AliVCluster *cluster = (AliVCluster*) evt->GetCaloCluster(icl);
+   if(!cluster->IsEMCAL()) {return feop;}
+
+       matchclsE = cluster->E();
+  
+
+   if(matchclsE<9999.0) feop = matchclsE/esdtrack->P();
+
+   return feop;
+
+}
+
+
 /*
-Double_t AliHFEpidEMCAL::MomentumEnergyMatch(const AliVParticle *track) const
+//____________________________________________________________________________________
+Double_t AliHFEpidEMCAL::MomentumEnergyMatchV1(const AliVParticle *track) const
 { // Returns e/p if an electron is matched
 
   Float_t  clsPos[3];
index a792438..e1848b3 100644 (file)
@@ -17,17 +17,6 @@ class AliVParticle;
 class AliPID;
 
 class AliHFEpidQAmanager;
-//class AliEMCALGeoUtils;
-//class AliEMCALGeometry;
-//class AliEMCALRecoUtils;
-/*
-class AliEMCALGeoUtils;
-class AliEMCALGeometry;
-#include "AliEMCALGeoParams.h"
-class AliEMCALRecoUtils;
-class AliEMCALPIDUtils;
-*/
-
 
 class AliHFEpidEMCAL : public AliHFEpidBase{
   public:
@@ -40,7 +29,8 @@ class AliHFEpidEMCAL : public AliHFEpidBase{
     virtual Bool_t    InitializePID();
     virtual Int_t     IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager *piqa) const;
       
-    //Double_t MomentumEnergyMatch(const AliVParticle *track) const;
+    //Double_t MomentumEnergyMatchV1(const AliVParticle *track) const;
+    Double_t MomentumEnergyMatchV2(const AliVParticle *track) const;
 
   protected:
     void Copy(TObject &ref) const;
index 660ad63..eeb0b64 100644 (file)
@@ -12,6 +12,7 @@ class AliHFEpidESD : public AliHFEpidBase{
     AliHFEpidESD();
     AliHFEpidESD(const AliHFEpidESD &ref);
     AliHFEpidESD &operator=(const AliHFEpidESD &ref);
+    virtual ~AliHFEpidESD();  
 
     virtual void InitializePID();
     virtual Int_t IsSelected(const AliVParticle *track);
index e686602..b743d06 100644 (file)
@@ -64,7 +64,8 @@ ClassImp(AliHFEpidQA)
     fTRDpidQA(NULL), 
     fOutput(NULL), 
     fESDpid(NULL),
-    fNNref(NULL)
+    fNNref(NULL),
+    fTRDTotalChargeInSlice0(kFALSE)
 {
   //
   // Default constructor
@@ -84,7 +85,8 @@ AliHFEpidQA::AliHFEpidQA(const AliHFEpidQA &ref):
   fTRDpidQA(NULL),
   fOutput(NULL), 
   fESDpid(NULL),
-  fNNref(NULL)
+  fNNref(NULL),
+  fTRDTotalChargeInSlice0(ref.fTRDTotalChargeInSlice0)
 {
   //
   // Copy constructor
@@ -131,6 +133,7 @@ void AliHFEpidQA::Copy(TObject &o) const {
 
   AliHFEpidQA &target = dynamic_cast<AliHFEpidQA &>(o);
   target.fMC = fMC;
+  target.fTRDTotalChargeInSlice0 = fTRDTotalChargeInSlice0;
 
   if(target.fESDpid) delete target.fESDpid;
   target.fESDpid = new AliESDpid;
@@ -176,6 +179,7 @@ void AliHFEpidQA::Init(){
   fV0pidMC->Init();
 
   fTRDpidQA = new AliHFEtrdPIDqa;
+  if(fTRDTotalChargeInSlice0) fTRDpidQA->SetTotalChargeInSlice0();
   fTRDpidQA->Init();
 
   fOutput = new AliHFEcollection("pidQA", "PID QA output");
@@ -937,7 +941,7 @@ void AliHFEpidQA::RecalculateTRDpid(AliAODTrack * /*track*/, Double_t * /*pidPro
   //  fTRDpidResponse->MakePID(track, pidProbs);
 }
 //___________________________________________________________________
-Float_t AliHFEpidQA::TOFbeta(AliESDtrack * const track) const {
+Float_t AliHFEpidQA::TOFbeta(const AliESDtrack * const track) const {
   // computes the TOF beta
   Double_t l = track->GetIntegratedLength();  // cm
   Double_t t = track->GetTOFsignal();
index 412ea4c..9793666 100644 (file)
@@ -69,9 +69,10 @@ class AliHFEpidQA : public TObject{
     void     SetMCEvent(AliMCEvent * const mc) { fMC = mc; };
     void     SetV0pidQA(Bool_t v0pidQA = kTRUE) { SetBit(kV0pidQA, v0pidQA); };
     void     SetRecalculateTRDpid(Bool_t recal = kTRUE) { SetBit(kRecalculateTRDpid, recal); };
+    void     SetTRDTotalChargeInSlice0() { fTRDTotalChargeInSlice0 = kTRUE; }
 
     void     SetESDpid(AliESDpid* const pid) { fESDpid = pid; }
-    Float_t  TOFbeta(AliESDtrack* const track) const;
+    Float_t  TOFbeta(const AliESDtrack* const track) const;
 
     void     CheckEvent();
     void     SetNNref(TFile *f) { fNNref = f; };
@@ -122,6 +123,7 @@ class AliHFEpidQA : public TObject{
  private:
     TFile             *fNNref;        // reference file for NN pid 
     TMultiLayerPerceptron *fNet[11];  //  reference networks
+    Bool_t fTRDTotalChargeInSlice0;     // Fix for Foreward/Backward compatibility
   
   ClassDef(AliHFEpidQA, 1)            // PID QA tool
 };
index 326613a..c0d376a 100644 (file)
@@ -49,16 +49,18 @@ AliHFEpidTPC::AliHFEpidTPC() :
   // add a list here
   AliHFEpidBase()
   , fLineCrossingsEnabled(0)
-  , fUpperSigmaCut(NULL)
-  , fLowerSigmaCut(NULL)
-  , fElectronMeanCorrection(NULL)
+  , fkElectronMeanCorrection(NULL)
+  , fHasCutModel(kFALSE)
   , fNsigmaTPC(3)
   , fRejectionEnabled(0)
-  , fPID(NULL)
 {
   //
   // default  constructor
   // 
+
+  memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
+  memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
+
   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
@@ -71,34 +73,32 @@ AliHFEpidTPC::AliHFEpidTPC(const char* name) :
   // add a list here
   AliHFEpidBase(name)
   , fLineCrossingsEnabled(0)
-  , fUpperSigmaCut(NULL)
-  , fLowerSigmaCut(NULL)
-  , fElectronMeanCorrection(NULL)
+  , fkElectronMeanCorrection(NULL)
+  , fHasCutModel(kFALSE)
   , fNsigmaTPC(3)
   , fRejectionEnabled(0)
-  , fPID(NULL)
 {
   //
   // default  constructor
   // 
   //
+  memset(fkUpperSigmaCut, 0, sizeof(const TF1 *) * 12);
+  memset(fkLowerSigmaCut, 0, sizeof(const TF1 *) * 12);
+
   memset(fRejection, 0, sizeof(Float_t) * 4 * AliPID::kSPECIES);
   memset(fLineCrossingSigma, 0, sizeof(Double_t) * AliPID::kSPECIES);
   memset(fPAsigCut, 0, sizeof(Float_t) * 2);
   memset(fNAsigmaTPC, 0, sizeof(Float_t) * 2);
-  fPID = new AliPID;
 }
 
 //___________________________________________________________________
 AliHFEpidTPC::AliHFEpidTPC(const AliHFEpidTPC &ref) :
   AliHFEpidBase("")
   , fLineCrossingsEnabled(0)
-  , fUpperSigmaCut(NULL)
-  , fLowerSigmaCut(NULL)
-  , fElectronMeanCorrection(NULL)
+  , fkElectronMeanCorrection(NULL)
+  , fHasCutModel(ref.fHasCutModel)
   , fNsigmaTPC(2)
   , fRejectionEnabled(0)
-  , fPID(NULL)
 {
   //
   // Copy constructor
@@ -125,12 +125,14 @@ void AliHFEpidTPC::Copy(TObject &o) const{
   AliHFEpidTPC &target = dynamic_cast<AliHFEpidTPC &>(o);
 
   target.fLineCrossingsEnabled = fLineCrossingsEnabled;
-  target.fUpperSigmaCut = fUpperSigmaCut;
-  target.fLowerSigmaCut = fLowerSigmaCut;
-  target.fElectronMeanCorrection = fElectronMeanCorrection;
+  target.fkElectronMeanCorrection = fkElectronMeanCorrection;
+  target.fHasCutModel = fHasCutModel;
   target.fNsigmaTPC = fNsigmaTPC;
   target.fRejectionEnabled = fRejectionEnabled;
-  target.fPID = new AliPID(*fPID);
+
+  memcpy(target.fkUpperSigmaCut, fkUpperSigmaCut, sizeof(const TF1 *) * 12);
+  memcpy(target.fkLowerSigmaCut, fkLowerSigmaCut, sizeof(const TF1 *) * 12);
+
   memcpy(target.fLineCrossingSigma, fLineCrossingSigma, sizeof(Double_t) * AliPID::kSPECIES);
   memcpy(target.fPAsigCut, fPAsigCut, sizeof(Float_t) * 2);
   memcpy(target.fNAsigmaTPC, fNAsigmaTPC, sizeof(Float_t) * 2);
@@ -143,7 +145,6 @@ AliHFEpidTPC::~AliHFEpidTPC(){
   //
   // Destructor
   //
-  if(fPID) delete fPID;
 }
 
 //___________________________________________________________________
@@ -195,8 +196,8 @@ Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
 
   // Check if we have an asymmetric sigma model set
   Int_t pdg = 0;
-  if(fUpperSigmaCut || fLowerSigmaCut){
-    pdg = CutSigmaModel(track->GetRecTrack(), anatype) ? 11 : 0;
+  if(fHasCutModel){
+    pdg = CutSigmaModel(track) ? 11 : 0;
   } else { 
     // Perform Asymmetric n-sigma cut if required, else perform symmetric TPC sigma cut
     Float_t p = 0.;
@@ -212,15 +213,19 @@ Int_t AliHFEpidTPC::IsSelected(const AliHFEpidObject *track, AliHFEpidQAmanager
 }
 
 //___________________________________________________________________
-Bool_t AliHFEpidTPC::CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const {
+Bool_t AliHFEpidTPC::CutSigmaModel(const AliHFEpidObject * const track) const {
   //
   // N SigmaCut using parametrization of the cuts
   //
   Bool_t isSelected = kTRUE;
-  Float_t nsigma = NumberOfSigmas(track, AliPID::kElectron, anaType);
-  Double_t p = GetP(track, anaType);
-  if(fUpperSigmaCut && nsigma > fUpperSigmaCut->Eval(p)) isSelected = kFALSE;
-  if(fLowerSigmaCut && nsigma < fLowerSigmaCut->Eval(p)) isSelected = kFALSE;
+  AliHFEpidObject::AnalysisType_t anatype = track->IsESDanalysis() ? AliHFEpidObject::kESDanalysis : AliHFEpidObject::kAODanalysis;
+  Float_t nsigma = NumberOfSigmas(track->GetRecTrack(), AliPID::kElectron, anatype);
+  Double_t p = GetP(track->GetRecTrack(), anatype);
+  Int_t centrality = track->IsPbPb() ? track->GetCentrality() + 1 : 0;
+  AliDebug(2, Form("Centrality: %d\n", centrality));
+  const TF1 *cutfunction;
+  if((cutfunction = fkUpperSigmaCut[centrality]) && nsigma > cutfunction->Eval(p)) isSelected = kFALSE;
+  if((cutfunction = fkLowerSigmaCut[centrality]) && nsigma < cutfunction->Eval(p)) isSelected = kFALSE;
   return isSelected;
 }
 
@@ -270,8 +275,8 @@ Double_t AliHFEpidTPC::NumberOfSigmas(const AliVParticle *track, AliPID::EPartic
     if(aodtrack && fAODpid) nSigmas = fAODpid->NumberOfSigmasTPC(aodtrack, species);
   }
   // Correct for the mean o
-  if(fElectronMeanCorrection)
-    nSigmas -= fElectronMeanCorrection->Eval(GetP(track, anaType));   
+  if(fkElectronMeanCorrection)
+    nSigmas -= fkElectronMeanCorrection->Eval(GetP(track, anaType));   
   return nSigmas;
 }
 
index e7bbe84..7a5f8af 100644 (file)
@@ -51,13 +51,15 @@ class AliHFEpidTPC : public AliHFEpidBase{
     void AddTPCdEdxLineCrossing(Int_t species, Double_t sigma);
     Bool_t HasAsymmetricSigmaCut() const { return TestBit(kAsymmetricSigmaCut);}
     Bool_t HasParticleRejection() const { return TestBit(kRejection); }
-    void SetElectronMeanCorrection(TF1 *electronLineCorrection) { fElectronMeanCorrection = electronLineCorrection; }
+    void SetElectronMeanCorrection(const TF1 * const electronLineCorrection) { fkElectronMeanCorrection = electronLineCorrection; }
     void SetTPCnSigma(Short_t nSigma) { fNsigmaTPC = nSigma; };
     inline void SetAsymmetricTPCsigmaCut(Float_t pmin, Float_t pmax, Float_t sigmaMin, Float_t sigmaMax);
     inline void SetRejectParticle(Int_t species, Float_t pmin, Float_t sigmaMin, Float_t pmax, Float_t sigmaMax);
 
-    void SetUpperSigmaCut(TF1 * const model) { fUpperSigmaCut = model; }
-    void SetLowerSigmaCut(TF1 * const model) { fLowerSigmaCut = model; }
+    void SetUpperSigmaCutDefault(const TF1 * const model) { fkUpperSigmaCut[0] = model; fHasCutModel = kTRUE; }
+    void SetUpperSigmaCutCentrality(const TF1 * const model, Int_t centralityBin) { if(centralityBin < 11) fkUpperSigmaCut[centralityBin+1] = model; fHasCutModel = kTRUE; }
+    void SetLowerSigmaCutDefault(const TF1 * const model) { fkLowerSigmaCut[0] = model; fHasCutModel = kTRUE; }
+    void SetLowerSigmaCutCentrality(const TF1 * const model, Int_t centralityBin) { if(centralityBin < 11) fkLowerSigmaCut[centralityBin+1] = model; fHasCutModel = kTRUE; }
 
     Double_t NumberOfSigmas(const AliVParticle *track, AliPID::EParticleType species, AliHFEpidObject::AnalysisType_t anatype) const;
     Double_t GetP(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
@@ -66,7 +68,7 @@ class AliHFEpidTPC : public AliHFEpidBase{
     void Copy(TObject &o) const;
     Int_t Reject(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
 
-    Bool_t CutSigmaModel(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anaType) const;
+    Bool_t CutSigmaModel(const AliHFEpidObject *anaType) const;
 
   private:
     enum{
@@ -75,15 +77,15 @@ class AliHFEpidTPC : public AliHFEpidBase{
     };
     Double_t fLineCrossingSigma[AliPID::kSPECIES];          // with of the exclusion point
     UChar_t fLineCrossingsEnabled;                          // Bitmap showing which line crossing is set
-    TF1 *fUpperSigmaCut;                                    // Upper Sigma Cut
-    TF1 *fLowerSigmaCut;                                    // Lower Sigma Cut
-    TF1 *fElectronMeanCorrection;                           // Correct the mean of the electron line position as function  of the momentum
+    const TF1 *fkUpperSigmaCut[12];                         // Upper Sigma Cut
+    const TF1 *fkLowerSigmaCut[12];                         // Lower Sigma Cut
+    const TF1 *fkElectronMeanCorrection;                    // Correct the mean of the electron line position as function  of the momentum
+    Bool_t fHasCutModel;                                    // Has cut model functions
     Float_t fPAsigCut[2];                                   // Momentum region where to perform asymmetric sigma cut
     Float_t fNAsigmaTPC[2];                                 // Asymmetric TPC Sigma band        
     Short_t fNsigmaTPC;                                     // TPC sigma band
     Float_t fRejection[4*AliPID::kSPECIES];                 // All informations for Particle Rejection, order pmin, sigmin, pmax, sigmax
     UChar_t fRejectionEnabled;                              // Bitmap for enabled particle rejection
-    AliPID *fPID;                                           //! PID Object
 
   ClassDef(AliHFEpidTPC, 1)   // TPC Electron ID class
 };
index e102d13..2720313 100644 (file)
@@ -48,6 +48,7 @@ AliHFEpidTRD::AliHFEpidTRD() :
   , fMinP(1.)
   , fElectronEfficiency(0.91)
   , fPIDMethod(kNN)
+  , fTotalChargeInSlice0(kFALSE)
 {
   //
   // default  constructor
@@ -62,6 +63,7 @@ AliHFEpidTRD::AliHFEpidTRD(const char* name) :
   , fMinP(1.)
   , fElectronEfficiency(0.91)
   , fPIDMethod(kNN)
+  , fTotalChargeInSlice0(kFALSE)
 {
   //
   // default  constructor
@@ -76,6 +78,7 @@ AliHFEpidTRD::AliHFEpidTRD(const AliHFEpidTRD &ref):
   , fMinP(ref.fMinP)
   , fElectronEfficiency(ref.fElectronEfficiency)
   , fPIDMethod(ref.fPIDMethod)
+  , fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
 {
   //
   // Copy constructor
@@ -106,6 +109,7 @@ void AliHFEpidTRD::Copy(TObject &ref) const {
   target.SetUseDefaultParameters(defaultParameters);
   target.fMinP = fMinP;
   target.fPIDMethod = fPIDMethod;
+  target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
   target.fElectronEfficiency = fElectronEfficiency;
   memcpy(target.fThreshParams, fThreshParams, sizeof(Double_t) * kThreshParams);
   AliHFEpidBase::Copy(ref);
@@ -322,13 +326,23 @@ Double_t AliHFEpidTRD::GetChargeLayer(const AliVParticle *track, UInt_t layer, A
   Double_t charge = 0.;
   if(anaType == AliHFEpidObject::kESDanalysis){
     const AliESDtrack *esdtrack = dynamic_cast<const AliESDtrack *>(track);
-    if(esdtrack)
-      for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
+    if(esdtrack){
+      // Distinction between old and new reconstruction: in the new reconstruction, the total charge is stored in slice 0, slices 1 to 8 are used for the slices for 
+      // the neural network. 
+      if(fTotalChargeInSlice0)
+        charge = esdtrack->GetTRDslice(static_cast<UInt_t>(layer), 0);
+      else
+       for(Int_t islice = 0; islice < esdtrack->GetNumberOfTRDslices(); islice++) charge += esdtrack->GetTRDslice(static_cast<UInt_t>(layer), islice);
+    }
   } else {
     const AliAODTrack *aodtrack = dynamic_cast<const AliAODTrack *>(track);
     AliAODPid *aoddetpid = aodtrack ? aodtrack->GetDetPid() : NULL;
-    if(aoddetpid)
-      for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
+    if(aoddetpid){
+      if(fTotalChargeInSlice0)
+        charge = aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices()];
+      else
+       for(Int_t islice = 0; islice < aoddetpid->GetTRDnSlices(); islice++) charge += aoddetpid->GetTRDsignal()[layer * aoddetpid->GetTRDnSlices() + islice];
+    }
   }
   return charge;
 }
@@ -375,7 +389,7 @@ Double_t AliHFEpidTRD::GetTRDSignalV1(const AliESDtrack *track, Float_t truncati
   Int_t indices[48];
   Int_t icnt = 0;
   for(Int_t idet = 0; idet < 6; idet++)
-    for(Int_t islice = 0; islice < kSlicePerLayer; islice++){
+    for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0 ; islice < kSlicePerLayer; islice++){
       AliDebug(2, Form("Chamber[%d], Slice[%d]: TRDSlice = %f", idet, islice, track->GetTRDslice(idet, islice)));
       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
       trdSlices[icnt++] = track->GetTRDslice(idet, islice) * kWeightSlice[islice];
@@ -409,7 +423,7 @@ Double_t AliHFEpidTRD::GetTRDSignalV2(const AliESDtrack *track, Float_t truncati
   Int_t indices[kLayers*kSlicesLow];
   Int_t cntLowTime=0, cntRemaining = 0;
   for(Int_t idet = 0; idet < 6; idet++)
-    for(Int_t islice = 0; islice < kSlicesLow+kSlicesHigh; islice++){
+    for(Int_t islice = fTotalChargeInSlice0 ? 1 : 0; islice < kSlicesLow+kSlicesHigh; islice++){
       if(TMath::Abs(track->GetTRDslice(idet, islice)) < fgkVerySmall) continue;;
       if(islice < kSlicesLow){
         AliDebug(3, Form("Part 1, Det[%d], Slice[%d], TRDSlice: %f", idet, islice, track->GetTRDslice(idet, islice)));
index 8043bcf..88a0d39 100644 (file)
@@ -65,6 +65,7 @@ class AliHFEpidTRD : public AliHFEpidBase{
     Bool_t IsCalculateTRDSignals() const { return TestBit(kTRDsignals); }
     Bool_t IsRenormalizeElPi() const { return TestBit(kTRDrenormalize); }
     void SetPIDMethod(PIDMethodTRD_t method) { fPIDMethod = method; };
+    void SetTotalChargeInSlice0() { fTotalChargeInSlice0 = kTRUE; }
     void SetRenormalizeElPi(Bool_t doRenorm = kTRUE) { if(doRenorm) SetBit(kTRDrenormalize, kTRUE); else SetBit(kTRDrenormalize, kFALSE);}
     void SetElectronEfficiency(Double_t electronEfficiency) { fElectronEfficiency = electronEfficiency; }
     void SetThresholdParameters(Double_t electronEff, Double_t *params);
@@ -96,6 +97,7 @@ class AliHFEpidTRD : public AliHFEpidBase{
     Double_t fElectronEfficiency;                           // Cut on electron efficiency
     PIDMethodTRD_t fPIDMethod;                              // PID Method: 2D Likelihood or Neural Network
     Double_t fThreshParams[kThreshParams];                  // Threshold parametrisation
+    Bool_t fTotalChargeInSlice0;                            // Flag for foreward/backward compatibility for the TRD total charge
   ClassDef(AliHFEpidTRD, 1)     // TRD electron ID class
 };
 #endif
index dd97f48..41e976c 100644 (file)
@@ -1401,7 +1401,7 @@ Int_t AliHFEsecVtx::GetElectronSource(Int_t iTrack)
 }
 
 //_______________________________________________________________________________________________
-Int_t AliHFEsecVtx::GetPDG(AliVTrack *track)
+Int_t AliHFEsecVtx::GetPDG(const AliVTrack *track)
 {
   //
   // get KF particle input pdg for mass hypothesis
@@ -1659,9 +1659,9 @@ void AliHFEsecVtx::MakeContainer(){
   fSecVtxList->AddAt(fSecvtxQA,1);
 
   for(Int_t ivar = 0; ivar < nDimPair; ivar++)
-    delete binEdgesPair[ivar];
+    delete[] binEdgesPair[ivar];
   for(Int_t ivar = 0; ivar < nDimSecvtx; ivar++)
-    delete binEdgesSecvtx[ivar];
+    delete[] binEdgesSecvtx[ivar];
 }
 
 //____________________________________________________________
@@ -1687,6 +1687,7 @@ void AliHFEsecVtx::MakeHistos(Int_t step){
   fSecVtxList->AddAt(new TH1F(hname+"ebgElec", "pT of e", iBin[0],binEdges[0]), step+4);
   fSecVtxList->AddAt(new TH1F(hname+"hcontaminElec", "pT of e", iBin[0],binEdges[0]), step+5);
   fSecVtxList->AddAt(new TH1F(hname+"elseElec", "pT of e", iBin[0],binEdges[0]), step+6);
+  delete[] binEdges[0];
 }
 
 //____________________________________________________________
index efe8a5f..3d2f6d4 100644 (file)
@@ -78,7 +78,7 @@ class AliHFEsecVtx : public TObject {
     Int_t GetPairOriginAOD(AliAODTrack* track1, AliAODTrack* track2); // return pair origin as a pdg code
     Int_t GetPairCode(const AliVTrack* const track1, const AliVTrack* const track2); // return corresponding pair code to pdg code
     Int_t GetElectronSource(Int_t mclabel); // return origin of the electron
-    Int_t GetPDG(AliVTrack *track);     // return pdg 
+    Int_t GetPDG(const AliVTrack *track);     // return pdg 
     void GetESDPID(const AliESDtrack *track, Int_t &recpid, Double_t &recprob); //return esd pid likelihood
     void GetPrimaryCondition();
     void RecalcPrimvtx(Int_t nkftrk, const Int_t * const, const AliKFParticle * const); //recalculate primary vertex
index 2080ca5..bddaf85 100644 (file)
@@ -25,7 +25,9 @@
 #include <TParticle.h>
 #include <TString.h>
 
+#include "AliAODTrack.h"
 #include "AliAODMCParticle.h"
+#include "AliESDtrack.h"
 #include "AliLog.h"
 #include "AliMCEvent.h"
 #include "AliMCParticle.h"
@@ -315,17 +317,16 @@ Int_t AliHFEsignalCuts::GetElecSource(const AliVParticle * const track) const {
     return 0;
   }
 
-  TString sourcetype = track->IsA()->GetName();
+  TClass *tracktype;
   const AliVParticle *mctrack = NULL;
   TParticle *mcpart = NULL;
-  if(!sourcetype.CompareTo("AliESDtrack") || !sourcetype.CompareTo("AliAODTrack")){
+  if((tracktype = track->IsA()) == AliESDtrack::Class() || tracktype == AliAODTrack::Class()){
     mctrack = fMC->GetTrack(TMath::Abs(track->GetLabel()));
   } else  mctrack = track;
   if(!mctrack) return 0;
 
-  TString mctype = mctrack->IsA()->GetName();
   Int_t eSource = 0;
-  if(!mctype.CompareTo("AliMCParticle")){
+  if(mctrack->IsA() == AliMCParticle::Class()){
     const AliMCParticle *esdmc = dynamic_cast<const AliMCParticle *>(mctrack);
     if(esdmc){
       mcpart = esdmc->Particle();
index 356e76a..1461984 100644 (file)
@@ -70,6 +70,7 @@ AliHFEspectrum::AliHFEspectrum(const char *name):
   fInclusiveSpectrum(kTRUE),
   fDumpToFile(kFALSE),
   fUnSetCorrelatedErrors(kTRUE),
+  fSetSmoothing(kFALSE),
   fIPanaHadronBgSubtract(kFALSE),
   fIPanaCharmBgSubtract(kFALSE),
   fIPanaConversionBgSubtract(kFALSE),
@@ -243,6 +244,8 @@ Bool_t AliHFEspectrum::Init(const AliHFEcontainer *datahfecontainer, const AliHF
     else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
     else if((fStepMC==(AliHFEcuts::kNcutStepsMCTrack + AliHFEcuts::kStepHFEcutsTRD - 1))) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepbeforePID");
     else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
+    
+    if(!mccorrelation) mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID");
   }
   else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterPID"); // we confirmed that we get same result by using it instead of correlationstepafterDE
   //else mccorrelation = mchfecontainer->GetCorrelationMatrix("correlationstepafterDE");
@@ -456,8 +459,8 @@ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
 
     if(fBeamType==1) {
 
-      TCanvas * ccorrected_allspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
-      ccorrected_allspectra->Divide(2,1);
+      TCanvas * ccorrectedallspectra = new TCanvas("correctedallspectra","correctedallspectra",1000,700);
+      ccorrectedallspectra->Divide(2,1);
       TLegend *legtotal = new TLegend(0.4,0.6,0.89,0.89);
       TLegend *legtotalg = new TLegend(0.4,0.6,0.89,0.89);
       
@@ -507,15 +510,15 @@ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
        titleenamednormalized += "[";
        Int_t nbEvents = 0;
        for(Int_t k = fLowBoundaryCentralityBinAtTheEnd[binc]; k < fHighBoundaryCentralityBinAtTheEnd[binc]; k++) {
-         //printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
+         printf("Number of events %d in the bin %d added!!!\n",fNEvents[k],k);
          nbEvents += fNEvents[k];
        }
-       //Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
-       //Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
-       //printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
-       //Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
-       //Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
-       //printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
+       Double_t lowedgega = cenaxisa->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+       Double_t upedgega = cenaxisa->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+       printf("Bin Low edge %f, up edge %f for a\n",lowedgega,upedgega);
+       Double_t lowedgegb = cenaxisb->GetBinLowEdge(fLowBoundaryCentralityBinAtTheEnd[binc]+1);
+       Double_t upedgegb = cenaxisb->GetBinUpEdge(fHighBoundaryCentralityBinAtTheEnd[binc]);
+       printf("Bin Low edge %f, up edge %f for b\n",lowedgegb,upedgegb);
        cenaxisa->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
        cenaxisb->SetRange(fLowBoundaryCentralityBinAtTheEnd[binc]+1,fHighBoundaryCentralityBinAtTheEnd[binc]);
        TCanvas * ccorrectedcheck = new TCanvas((const char*) titleec,(const char*) titleec,1000,700);
@@ -563,7 +566,7 @@ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
        legcorrectedud->AddEntry(aftersun,"Corrected","p");
        legcorrectedud->AddEntry(aftersdn,"Alltogether","p");
        legcorrectedud->Draw("same");
-       ccorrected_allspectra->cd(1);
+       ccorrectedallspectra->cd(1);
        gPad->SetLogy();
        TH1D *aftersunn = (TH1D *) aftersun->Clone();
        aftersunn->SetMarkerStyle(stylee[binc]);
@@ -571,7 +574,7 @@ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
        if(binc==0) aftersunn->Draw("AP");
        else aftersunn->Draw("P");
        legtotal->AddEntry(aftersunn,(const char*) titlee,"p");
-       ccorrected_allspectra->cd(2);
+       ccorrectedallspectra->cd(2);
        gPad->SetLogy();
        TH1D *aftersdnn = (TH1D *) aftersdn->Clone();
        aftersdnn->SetMarkerStyle(stylee[binc]);
@@ -592,9 +595,9 @@ Bool_t AliHFEspectrum::Correct(Bool_t subtractcontamination){
        ratiocorrectedbinc->Draw();
       }
 
-      ccorrected_allspectra->cd(1);
+      ccorrectedallspectra->cd(1);
       legtotal->Draw("same");
-      ccorrected_allspectra->cd(2);
+      ccorrectedallspectra->cd(2);
       legtotalg->Draw("same");
       
       cenaxisa->SetRange(0,nbbin);
@@ -1516,7 +1519,7 @@ TList *AliHFEspectrum::Unfold(AliCFDataGrid* const bgsubpectrum){
   //unfolding.SetUseCorrelatedErrors();
   if(fUnSetCorrelatedErrors) unfolding.UnsetCorrelatedErrors();
   unfolding.SetMaxNumberOfIterations(fNumberOfIterations);
-  unfolding.UseSmoothing();
+  if(fSetSmoothing) unfolding.UseSmoothing();
   unfolding.Unfold();
 
   // Results
@@ -1802,10 +1805,10 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1(TH1 *input,Int_t i) const {
       dN = input->GetBinError(ibin);
 
       // New point
-      nCorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents[i]) * 1/(2. * TMath::Pi() * p) * n;
+      nCorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * 1./(2. * TMath::Pi() * p) * n;
       errdN = 1./(2. * TMath::Pi() * p);
       errdp = 1./(2. * TMath::Pi() * p*p) * n;
-      dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+      dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(fNEvents[i]) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
       
       spectrumNormalized->SetPoint(point, p, nCorr);
       spectrumNormalized->SetPointError(point, dp, dNcorr);
@@ -1845,10 +1848,10 @@ TGraphErrors *AliHFEspectrum::NormalizeTH1N(TH1 *input,Int_t normalization) cons
       dN = input->GetBinError(ibin);
 
       // New point
-      nCorr = 0.5 * 1/1.6 * 1./(Double_t)(normalization) * 1/(2. * TMath::Pi() * p) * n;
+      nCorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * 1./(2. * TMath::Pi() * p) * n;
       errdN = 1./(2. * TMath::Pi() * p);
       errdp = 1./(2. * TMath::Pi() * p*p) * n;
-      dNcorr = 0.5 * 1/1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
+      dNcorr = 0.5 * 1./1.6 * 1./(Double_t)(normalization) * TMath::Sqrt(errdN * errdN * dN *dN + errdp *errdp * dp *dp);
       
       spectrumNormalized->SetPoint(point, p, nCorr);
       spectrumNormalized->SetPointError(point, dp, dNcorr);
index 71705fd..352a568 100644 (file)
@@ -86,6 +86,7 @@ class AliHFEspectrum : public TNamed{
     void SetNbDimensions(Int_t nbDimensions) { fNbDimensions = nbDimensions; };
     void SetChargeChoosen(Int_t chargechoosen) {fChargeChoosen = chargechoosen; };
     void SetUnSetCorrelatedErrors(Bool_t unsetcorrelatederrors) {fUnSetCorrelatedErrors = unsetcorrelatederrors;};
+    void SetSmoothing(Bool_t setSmoothing) {fSetSmoothing = setSmoothing;};
 
     void SetNCentralityBinAtTheEnd(Int_t nCentralityBinAtTheEnd) {fNCentralityBinAtTheEnd = nCentralityBinAtTheEnd; };
     void SetLowHighBoundaryCentralityBinAtTheEnd(Int_t low, Int_t high, Int_t i) { fLowBoundaryCentralityBinAtTheEnd[i] = low; fHighBoundaryCentralityBinAtTheEnd[i] = high;};
@@ -141,6 +142,7 @@ class AliHFEspectrum : public TNamed{
     Bool_t fDumpToFile;           // Write Result in a file
 
     Bool_t fUnSetCorrelatedErrors;    // Unset correlated errors
+    Bool_t fSetSmoothing;             // Set smoothing
 
     Bool_t fIPanaHadronBgSubtract;     // Hadron background subtraction
     Bool_t fIPanaCharmBgSubtract;      // Charm background subtraction 
index a0eaa43..6547af7 100644 (file)
@@ -51,10 +51,11 @@ AliHFEtaggedTrackAnalysis::AliHFEtaggedTrackAnalysis():
   , fCuts(NULL)
   , fCFM(NULL)
   , fQAhistos(NULL)
-  , fCentralityF(0.)
+  , fCentralityF(0)
   , fClean(kFALSE)
   , fMagneticField(0.0)
   , fVariablesTRD(kFALSE)
+  , fIsPbPb(kFALSE)
 {
   //
   // Dummy constructor
@@ -71,10 +72,11 @@ AliHFEtaggedTrackAnalysis::AliHFEtaggedTrackAnalysis(const char *name):
   , fCuts(NULL)
   , fCFM(NULL)
   , fQAhistos(NULL)
-  , fCentralityF(0.)
+  , fCentralityF(0)
   , fClean(kFALSE)
   , fMagneticField(0.0)
   , fVariablesTRD(kFALSE)
+  , fIsPbPb(kFALSE)
 {
   //
   // Default constructor
@@ -104,6 +106,7 @@ AliHFEtaggedTrackAnalysis::AliHFEtaggedTrackAnalysis(const AliHFEtaggedTrackAnal
   , fClean(ref.fClean)
   , fMagneticField(ref.fMagneticField)
   , fVariablesTRD(ref.fVariablesTRD)
+  , fIsPbPb(ref.fIsPbPb)
 {
   //
   // Copy constructor
@@ -131,6 +134,7 @@ AliHFEtaggedTrackAnalysis &AliHFEtaggedTrackAnalysis::operator=(const AliHFEtagg
     fClean = ref.fClean;
     fMagneticField = ref.fMagneticField;
     fVariablesTRD = ref.fVariablesTRD;
+    fIsPbPb = ref.fIsPbPb;
 
     if(ref.fContainer) InitContainer();
    
@@ -272,6 +276,8 @@ void AliHFEtaggedTrackAnalysis::ProcessTrack(AliVParticle *track, Int_t abinitio
      hfetrack.SetRecTrack(track);
      hfetrack.SetAbInitioPID(abinitioPID);
      hfetrack.SetCentrality(fCentralityF);
+     if(fIsPbPb) hfetrack.SetPbPb();
+     else hfetrack.SetPP();
      fPID->SetVarManager(fVarManager);
      fPID->IsSelected(&hfetrack, fContainer, "taggedTrackContainer", fPIDqa);
    }
index bcd394b..e2c5536 100644 (file)
@@ -49,12 +49,16 @@ class AliHFEtaggedTrackAnalysis : public TNamed{
     TList * GetCutQA() const;
     AliHFEcollection * GetQAcollection() const { return fQAhistos; }
     Bool_t  GetClean() const { return fClean; }; 
+    Bool_t IsPbPb() const { return fIsPbPb; }
+    Bool_t IsPP() const { return !fIsPbPb; }
     Double_t GetMagneticField() const { return fMagneticField; };
     AliHFEvarManager *GetVarManager() const { return fVarManager; }
 
     void SetCuts(AliHFEcuts *cuts);
     void SetPID(AliHFEpid *pid);
-    void SetCentrality(Float_t centrality) { fCentralityF = centrality; };
+    void SetCentrality(Int_t centrality) { fCentralityF = centrality; };
+    void SetPbPb(){ fIsPbPb = kTRUE; }
+    void SetPP() { fIsPbPb = kFALSE; }
     void SetClean(Bool_t clean) { fClean = clean; };
     void SetMagneticField(Double_t magneticField) { fMagneticField = magneticField; };
     void SetVariablesTRD(Bool_t variablesTRD) { fVariablesTRD = variablesTRD; };
@@ -71,10 +75,11 @@ class AliHFEtaggedTrackAnalysis : public TNamed{
     AliHFEcuts          *fCuts;         // Single track cuts
     AliCFManager        *fCFM;          // CF Manager used for the track filtering
     AliHFEcollection    *fQAhistos;     // QA histos
-    Float_t              fCentralityF;  // Centrality
+    Int_t              fCentralityF;  // Centrality
     Bool_t               fClean;        // Clean
     Double_t             fMagneticField; // Magnetic field
-    Bool_t               fVariablesTRD;  //  Use phi angle at the first plane of the TRD 
+    Bool_t               fVariablesTRD;  //  Use phi angle at the first plane of the TRD
+    Bool_t               fIsPbPb;        // Analysis Type: pp or PbPb 
     
   ClassDef(AliHFEtaggedTrackAnalysis, 0)
 };
index a2af7ba..bb616c7 100644 (file)
@@ -51,6 +51,7 @@ ClassImp(AliHFEtpcPIDqa)
 AliHFEtpcPIDqa::AliHFEtpcPIDqa():
     AliHFEdetPIDqa()
   , fHistos(NULL)
+  , fBrowseCentrality(-1)
 {
   //
   // Dummy constructor
@@ -61,6 +62,7 @@ AliHFEtpcPIDqa::AliHFEtpcPIDqa():
 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
     AliHFEdetPIDqa(name, "QA for TPC")
   , fHistos(NULL)
+  , fBrowseCentrality(-1)
 {
   //
   // Default constructor
@@ -71,6 +73,7 @@ AliHFEtpcPIDqa::AliHFEtpcPIDqa(const char* name):
 AliHFEtpcPIDqa::AliHFEtpcPIDqa(const AliHFEtpcPIDqa &o):
     AliHFEdetPIDqa(o)
   , fHistos(NULL)
+  , fBrowseCentrality(o.fBrowseCentrality)
 {
   //
   // Copy constructor
@@ -107,6 +110,7 @@ void AliHFEtpcPIDqa::Copy(TObject &o) const {
     target.fHistos = NULL;
   }
   if(fHistos) target.fHistos = new AliHFEcollection(*fHistos);
+  target.fBrowseCentrality = fBrowseCentrality;
 }
 
 //_________________________________________________________
@@ -153,11 +157,11 @@ void AliHFEtpcPIDqa::Browse(TBrowser *b){
       TH2 *hptr = NULL; 
       for(Int_t ispec = 0; ispec < AliPID::kSPECIES+1; ispec++){
         for(Int_t istep = 0; istep < 2; istep++){
-          hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          hptr->SetName(Form("hTPCdEdx%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+          hptr = MakeSpectrumdEdx(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
+          hptr->SetName(Form("hTPCdEdx%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After" , fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
           listdEdx->Add(hptr);
-          hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          hptr->SetName(Form("hTPCnsigma%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+          hptr = MakeSpectrumNSigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], fBrowseCentrality);
+          hptr->SetName(Form("hTPCnsigma%s%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After", fBrowseCentrality == -1 ? "MinBias" : Form("Cent%d", fBrowseCentrality)));
           listNsigma->Add(hptr);
         }
       }
@@ -245,7 +249,7 @@ Double_t AliHFEtpcPIDqa::GetTPCsignal(const AliVParticle *track, AliHFEpidObject
 }
 
 //_________________________________________________________
-TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species){
+TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
   //
   // Plot the Spectrum
   //
@@ -254,13 +258,34 @@ TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t speci
   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
   if(species >= 0 && species < AliPID::kSPECIES)
     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
-  TH2 *hTmp = hSignal->Projection(2,1);
   TString hname = Form("hTPCsignal%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
           htitle = Form("TPC dE/dx Spectrum %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
   if(species > -1){
     hname += AliPID::ParticleName(species);
     htitle += Form(" for %ss", AliPID::ParticleName(species));
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+    if(hSignal->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo= kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form(" Centrality %d", centralityClass);
+    }
+  } else {
+    centname = centtitle = "MinBias";
+    hasCentralityInfo= kFALSE;
+  }
+  hname += centtitle;
+  htitle += centtitle;
+
+  TH2 *hTmp = hSignal->Projection(2,1);
   hTmp->SetName(hname.Data());
   hTmp->SetTitle(htitle.Data());
   hTmp->SetStats(kFALSE);
@@ -268,11 +293,12 @@ TH2 *AliHFEtpcPIDqa::MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t istep, Int_t speci
   hTmp->GetYaxis()->SetTitle("TPC signal [a.u.]");
   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
+  if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
   return hTmp;
 }
 
 //_________________________________________________________
-TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species){
+TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t species, Int_t centralityClass){
   //
   // Plot the Spectrum
   //
@@ -281,13 +307,34 @@ TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t spe
   hSignal->GetAxis(3)->SetRange(istep + 1, istep + 1);
   if(species >= 0 && species < AliPID::kSPECIES)
     hSignal->GetAxis(0)->SetRange(2 + species, 2 + species);
-  TH2 *hTmp = hSignal->Projection(2,1);
   TString hname = Form("hTPCsigma%s", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after"), 
           htitle = Form("TPC dE/dx Spectrum[#sigma] %s selection", istep == AliHFEdetPIDqa::kBeforePID ? "before" : "after");
   if(species > -1){
     hname += AliPID::ParticleName(species);
     htitle += Form(" for %ss", AliPID::ParticleName(species));
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+    if(hSignal->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo= kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      hSignal->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form(" Centrality %d", centralityClass);
+    }
+  } else {
+    centname = centtitle = "MinBias";
+    hasCentralityInfo= kFALSE;
+  }
+  hname += centtitle;
+  htitle += centtitle;
+
+  TH2 *hTmp = hSignal->Projection(2,1);
   hTmp->SetName(hname.Data());
   hTmp->SetTitle(htitle.Data());
   hTmp->SetStats(kFALSE);
@@ -295,6 +342,7 @@ TH2 *AliHFEtpcPIDqa::MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t istep, Int_t spe
   hTmp->GetYaxis()->SetTitle("TPC dE/dx - <dE/dx>|_{el} [#sigma]");
   hSignal->GetAxis(3)->SetRange(0, hSignal->GetAxis(3)->GetNbins());
   hSignal->GetAxis(0)->SetRange(0, hSignal->GetAxis(0)->GetNbins());
+  if(hasCentralityInfo) hSignal->GetAxis(4)->SetRange(0, hSignal->GetAxis(4)->GetNbins());
   return hTmp;
 }
 
index dc30940..7e631a1 100644 (file)
@@ -48,15 +48,18 @@ class AliHFEtpcPIDqa : public AliHFEdetPIDqa{
     virtual void Initialize();
     virtual void ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step);
 
+    void SetBrowseCentrality(Int_t browseCentrality) { browseCentrality < 11  && browseCentrality >= -1 ? fBrowseCentrality = browseCentrality : -1;} // *MENU*
+
     AliHFEcollection *GetHistograms() const { return fHistos; }
-    TH2 *MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
-    TH2 *MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
+    TH2 *MakeSpectrumdEdx(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+    TH2 *MakeSpectrumNSigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
 
   protected:
     Double_t GetTPCsignal(const AliVParticle *track, AliHFEpidObject::AnalysisType_t anatype);
 
   private:
     AliHFEcollection *fHistos;        // Container for Histograms
+    Int_t fBrowseCentrality;          // Centrality Class for Browser
 
     ClassDef(AliHFEtpcPIDqa, 1);
 };
index 7a9e1b4..489a996 100644 (file)
@@ -511,7 +511,7 @@ void AliHFEtrackFilter::SetMC(AliMCEvent * const mc){
 }
 
 //__________________________________________________________________
-void AliHFEtrackFilter::SetRecEvent(AliVEvent *rec){
+void AliHFEtrackFilter::SetRecEvent(const AliVEvent *rec){
   //
   // Publish MC event to the single cut steps
   //
index 0017d7d..e0b7185 100644 (file)
@@ -46,7 +46,7 @@ class AliHFEtrackFilter : public TNamed{
     // Base Functionality
     void AddCutStep(AliHFEcutStep *cuts);
     void GenerateCutSteps();
-    void SetRecEvent(AliVEvent *rec);
+    void SetRecEvent(const AliVEvent *rec);
     AliHFEcutStep *GetCutStep(Int_t istep);
     AliHFEcutStep *GetCutStep(const Char_t *name);
     void FilterTracks(const AliESDEvent *const esd);
index 3b24486..97722c7 100644 (file)
@@ -79,7 +79,8 @@ AliHFEtrdPIDqa::AliHFEtrdPIDqa():
   fProtonEfficiencies(NULL),
   fKaonEfficiencies(NULL),
   fThresholds(NULL),
-  fShowMessage(kFALSE)
+  fShowMessage(kFALSE),
+  fTotalChargeInSlice0(kFALSE)
 {
   //
   // Default Constructor
@@ -95,7 +96,8 @@ AliHFEtrdPIDqa::AliHFEtrdPIDqa(const Char_t *name):
   fProtonEfficiencies(NULL),
   fKaonEfficiencies(NULL),
   fThresholds(NULL),
-  fShowMessage(kFALSE)
+  fShowMessage(kFALSE),
+  fTotalChargeInSlice0(kFALSE)
 {
   //
   // Main Constructor
@@ -111,7 +113,8 @@ AliHFEtrdPIDqa::AliHFEtrdPIDqa(const AliHFEtrdPIDqa &ref):
   fProtonEfficiencies(NULL),
   fKaonEfficiencies(NULL),
   fThresholds(NULL),
-  fShowMessage(kFALSE)
+  fShowMessage(kFALSE),
+  fTotalChargeInSlice0(ref.fTotalChargeInSlice0)
 {
   //
   // Copy constructor
@@ -150,7 +153,8 @@ void AliHFEtrdPIDqa::Copy(TObject &ref) const{
 
   AliHFEtrdPIDqa &target = dynamic_cast<AliHFEtrdPIDqa &>(ref);
   target.fTRDpid = fTRDpid;
-  target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());
+  target.fHistos = dynamic_cast<AliHFEcollection *>(fHistos->Clone());  
+  target.fTotalChargeInSlice0 = fTotalChargeInSlice0;
 }
 
 //__________________________________________________________________
@@ -289,8 +293,8 @@ void AliHFEtrdPIDqa::CreateHistoTruncatedMean(){
 
   fHistos->CreateTHnSparse("fTRDtruncMean","TRD TruncatedMean studies", kQuantitiesTruncMean, nbins, binMin, binMax);
   fHistos->BinLogAxis("fTRDtruncMean", kP);
-  fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 500, 0, 2000);
-  fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 500, 0, 2000);
+  fHistos->CreateTH2F("fTRDslicesPions","TRD dEdx per slice for Pions", 8, 0, 8, 2000, 0, 8000);
+  fHistos->CreateTH2F("fTRDslicesElectrons","TRD dEdx per slice for Electrons", 8, 0, 8, 2000, 0, 8000);
 }
 
 
@@ -420,7 +424,8 @@ void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t speci
     quantitiesdEdx[kP] = track->GetTRDmomentum(iplane);
     dEdxSum = 0.;
     for(Int_t islice = 0; islice < nSlices; islice++){
-      qSlice = track->GetTRDslice(iplane, islice);
+      if(fTotalChargeInSlice0 && islice >= 7) break;
+      qSlice = track->GetTRDslice(iplane, fTotalChargeInSlice0 ? islice + 1 : islice);  // hack by mfasel: For data with the new reconstruction, slice 0 is used to store the total charge, the total number of slices is 7 instead of 8
       if(qSlice > 1e-1){
         // cut out 0 slices
         nSlicesNonZero++;
@@ -434,7 +439,7 @@ void AliHFEtrdPIDqa::FillTRDQAplots(const AliESDtrack * const track, Int_t speci
       }
     }
     quantitiesdEdx[kNonZeroSlices] = nSlicesNonZero;
-    quantitiesdEdx[kdEdx] = dEdxSum;
+    quantitiesdEdx[kdEdx] = fTotalChargeInSlice0 ? track->GetTRDslice(iplane, 0) : dEdxSum; // hack by mfasel: In the new reconstruction, the total charge is stored in the first slice, in the old reconstruction it has to be calculated from the slice charges.     
     if(dEdxSum) ntrackletsNonZero++;
     // Fill dEdx histogram
     if(dEdxSum > 1e-1) fHistos->Fill("fQAdEdx", quantitiesdEdx); // Cut out 0 entries
@@ -465,14 +470,17 @@ void AliHFEtrdPIDqa::FinishAnalysis(){
   if(!fPionEfficiencies){
     fPionEfficiencies = new TList;
     fPionEfficiencies->SetName("pionEfficiencies");
+    fPionEfficiencies->SetOwner();
   }
   if(!fProtonEfficiencies){
     fProtonEfficiencies = new TList;
     fProtonEfficiencies->SetName("protonEfficiencies");
+    fProtonEfficiencies->SetOwner();
   }
   if(!fThresholds){
     fThresholds = new TList;
     fThresholds->SetName("thresholds");
+    fThresholds->SetOwner();
   }
 
   for(Int_t itr = 4; itr <= 6; itr++){
@@ -500,7 +508,7 @@ void AliHFEtrdPIDqa::StoreResults(const Char_t *filename){
 }
 
 //__________________________________________________________________
-void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
+void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name, Double_t lowerLimit, Double_t upperLimit){
   //
   // Fit the threshold histograms with the given parametrisation
   // and store the TF1 in the file
@@ -550,7 +558,7 @@ void AliHFEtrdPIDqa::SaveThresholdParameters(const Char_t *name){
 
       threshhist = dynamic_cast<TGraph *>(lHistos->FindObject(Form("eff%d", static_cast<Int_t>(fgkElectronEff[ieff] * 100))));
       if(!threshhist) continue;
-      threshparam = MakeThresholds(threshhist);
+      threshparam = MakeThresholds(threshhist, lowerLimit, upperLimit);
       threshparam->SetName(Form("thresh_%d_%d", itracklet, static_cast<Int_t>(fgkElectronEff[ieff] * 100)));
       lFormulas->Add(threshparam);
     }
@@ -600,9 +608,9 @@ void AliHFEtrdPIDqa::AnalyseNTracklets(Int_t nTracklets){
   hLikeTRD->GetAxis(kNTracklets)->SetRange(0, hLikeTRD->GetAxis(kNTracklets)->GetNbins());
 
   // Prepare List for output
-  TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets));
-  TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets));
-  TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets));
+  TList *listPions = new TList; listPions->SetName(Form("%dTracklets", nTracklets)); listPions->SetOwner();
+  TList *listProtons = new TList; listProtons->SetName(Form("%dTracklets", nTracklets)); listProtons->SetOwner();
+  TList *listThresholds = new TList; listThresholds->SetName(Form("%dTracklets", nTracklets)); listThresholds->SetOwner();
   fPionEfficiencies->Add(listPions);
   fProtonEfficiencies->Add(listProtons);
   fThresholds->Add(listThresholds);
@@ -776,7 +784,7 @@ void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax,
 
     // Optionally do Fit
     if(doFit){
-      threshfit = MakeThresholds(tr);
+      threshfit = MakeThresholds(tr, pmin, pmax);
       threshfit->SetLineColor(kBlack);
       threshfit->Draw("same");
     }
@@ -791,13 +799,13 @@ void AliHFEtrdPIDqa::DrawTracklet(Int_t itracklet, Double_t pmin, Double_t pmax,
 }
 
 //__________________________________________________________________
-TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist){
+TF1 *AliHFEtrdPIDqa::MakeThresholds(TGraph *threshist, Double_t lowerLimit, Double_t upperLimit){
   //
   // Create TF1 containing the threshold parametrisation
   //
 
   TF1 *threshparam = new TF1("thresh", "1-[0]-[1]*x-[2]*TMath::Exp(-[3]*x)", 0.1, 10);
-  threshist->Fit(threshparam, "NE", "", 0.1, 3.5);
+  threshist->Fit(threshparam, "NE", "", lowerLimit, upperLimit);
   return threshparam;
 }
 
index cbdef06..e0792c8 100644 (file)
@@ -69,8 +69,9 @@ class AliHFEtrdPIDqa : public TNamed{
     void Init();
     void FinishAnalysis();
     void ShowMessages() { fShowMessage = kTRUE; }
+    void SetTotalChargeInSlice0() { fTotalChargeInSlice0 = kTRUE; }
     void StoreResults(const Char_t *filename = "HFEtrdPIDqa.root");
-    void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root");
+    void SaveThresholdParameters(const Char_t * filename = "TRD.Thresholds.root", Double_t lowerLimit = 0.5, Double_t upperLimit = 3.5);
 
     void DrawTracklet(Int_t tracklet, Double_t pmin = 0., Double_t pmax = 0., Bool_t doFit = kFALSE);
     void ClearLists();
@@ -134,7 +135,7 @@ class AliHFEtrdPIDqa : public TNamed{
     void AnalyseNTracklets(Int_t nTracklets);
     Int_t GetThresholdBin(const TH1 * const input, Double_t efficiency);
     Bool_t CalculateEfficiency(const TH1 * const input, Int_t threshbin, Double_t *params);
-    TF1 *MakeThresholds(TGraph *input);
+    TF1 *MakeThresholds(TGraph *input, Double_t lowerLimit, Double_t upperLimit);
 
     void CreateLikelihoodHistogram();
     void CreateQAHistogram();
@@ -160,6 +161,7 @@ class AliHFEtrdPIDqa : public TNamed{
     TList *fThresholds;           //! List for Threshold Graphs
 
     Bool_t fShowMessage;         // Display debug messages
+    Bool_t fTotalChargeInSlice0;  // Flag for Foreward/Backward compatibility in TRD total charge calculation
   
   ClassDef(AliHFEtrdPIDqa, 3)     // QA class for TRD PID 
 };
index 839e4d4..28b20f2 100644 (file)
@@ -134,25 +134,44 @@ void AliHFEtrdPIDqaV1::Browse(TBrowser *b){
       listTPCnsigma->SetOwner();
 
       TH2 *hptr = NULL; 
-      for(Int_t ispec = 0; ispec < 4; ispec++){
-        for(Int_t istep = 0; istep < 2; istep++){
-          hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          if(hptr){
-            hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
-            listTM->Add(hptr);
+      TList *lctm, *lclike, *lccharge, *lcTPCsig;
+      for(Int_t icent = -1; icent < 11; icent++){
+        lctm = new TList;
+        lctm->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+        lctm->SetOwner();
+        listTM->Add(lctm);
+        lclike = new TList;
+        lclike->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+        lclike->SetOwner();
+        listLike->Add(lclike);
+        lccharge = new TList;
+        lccharge->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+        lccharge->SetOwner();
+        listCharge->Add(lccharge);
+        lcTPCsig = new TList;
+        lcTPCsig->SetName(icent < 0 ? "MinBias" : Form("Centrality%d", icent));
+        lcTPCsig->SetOwner();
+        listTPCnsigma->Add(lcTPCsig);
+        for(Int_t ispec = 0; ispec < 4; ispec++){
+          for(Int_t istep = 0; istep < 2; istep++){
+            hptr = MakeTRDspectrumTM(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+            if(hptr){
+              hptr->SetName(Form("hTRDtm%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+              lctm->Add(hptr);
+            }
+            hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+            hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+            lclike->Add(hptr);
+            hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+            hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+            lccharge->Add(hptr);
+            hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec], icent);
+            hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
+            lcTPCsig->Add(hptr);
           }
-          hptr = MakeTRDlikelihoodDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          hptr->SetName(Form("hTRDlike%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
-          listLike->Add(hptr);
-          hptr = MakeTRDchargeDistribution(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          hptr->SetName(Form("hTRDcharge%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
-          listCharge->Add(hptr);
-          hptr = MakeTPCspectrumNsigma(static_cast<AliHFEdetPIDqa::EStep_t>(istep), specind[ispec]);
-          hptr->SetName(Form("hTPCspectrum%s%s", specnames[ispec].Data(), istep == 0 ? "Before" : "After"));
-          listTPCnsigma->Add(hptr);
         }
       }
-      
+        
       b->Add(listTM, "Projections Truncated Mean");
       b->Add(listLike, "Projections Likelihood distribution");
       b->Add(listCharge, "Projections Tracklet Charge");
@@ -188,27 +207,27 @@ void AliHFEtrdPIDqaV1::Initialize(){
   Int_t nBinsTPCSigma[5] = {kPIDbins, kPbins, tpcSigmaBins, kSteps, kCentralityBins};
   Double_t minTPCSigma[5] = {kMinPID, kMinP, -12., 0., 0.};
   Double_t maxTPCSigma[5] = {kMaxPID, kMaxP, 12., 2., 11.};
-  fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 4, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
+  fHistos->CreateTHnSparse("hTPCsigma", "TPC sigma; species p [GeV/c]; TPC dEdx - <dE/dx>|_{el} [#sigma]; selection step", 5, nBinsTPCSigma, minTPCSigma, maxTPCSigma);
   fHistos->Sumw2("hTPCsigma");
   // Create Monitoring histogram for the Likelihood distribution
   Int_t nBinsTRDlike[5] = {kPIDbins, kPbins, trdLikelihoodBins, kSteps, kCentralityBins};
   Double_t minTRDlike[5] = {kMinPID, kMinP, 0., 0., 0.};
   Double_t maxTRDlike[5] = {kMaxPID, kMaxP, 1., 2., 11.};
-  fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 4, nBinsTRDlike, minTRDlike, maxTRDlike);
+  fHistos->CreateTHnSparse("hTRDlikelihood", "TRD Likelihood Distribution; species; p [GeV/c]; TRD electron Likelihood; selection step", 5, nBinsTRDlike, minTRDlike, maxTRDlike);
   fHistos->Sumw2("hTRDlikelihood");
   // Create Monitoring histogram for the TRD total charge
   const Int_t kTRDchargeBins = 1000;
   Int_t nBinsTRDcharge[5] = {kPIDbins, kPbins, kTRDchargeBins, kSteps, kCentralityBins};
   Double_t minTRDcharge[5] = {kMinPID, kMinP, 0., 0., 0.};
   Double_t maxTRDcharge[5] = {kMaxPID, kMaxP, 100000., 2., 11.};
-  fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 4, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
+  fHistos->CreateTHnSparse("hTRDcharge", "Total TRD charge; species; p [GeV/c]; TRD charge [a.u.]; selection step", 5, nBinsTRDcharge, minTRDcharge, maxTRDcharge);
   fHistos->Sumw2("hTRDcharge");
   // Monitoring of the TRD truncated mean according to version 1
   const Int_t kTRDtmBins = 1000;
   Int_t nBinsTRDtm[5] = {kPIDbins, kPbins, kTRDtmBins, kSteps, kCentralityBins};
   Double_t minTRDtm[5] = {kMinPID, kMinP, 0., 0., 0.};
   Double_t maxTRDtm[5] = {kMaxPID, kMaxP, 20000., 2., 11.};
-  fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 4, nBinsTRDtm, minTRDtm, maxTRDtm);
+  fHistos->CreateTHnSparse("hTRDtruncatedMean", "TRD truncated Mean; species; p [GeV/c]; TRD signal [a.u.]; selection step", 5, nBinsTRDtm, minTRDtm, maxTRDtm);
   fHistos->Sumw2("hTRDtruncatedMean");
 }
 
@@ -242,12 +261,13 @@ void AliHFEtrdPIDqaV1::ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa
   }
   for(UInt_t ily = 0; ily < 6; ily++){
     container[2] = trdpid ? trdpid->GetChargeLayer(track->GetRecTrack(), ily, anatype) : 0;
+    if(container[2] < 1e-3) continue; // Filter out 0 entries
     fHistos->Fill("hTRDcharge", container);
   }
 }
 
 //_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
   //
   // Get the TPC control histogram for the TRD selection step (either before or after PID)
   //
@@ -260,6 +280,25 @@ TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t
     // cut on species (if available)
     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+    if(histo->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo = kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form("Centrality %d", centralityClass);
+    }
+  } else {
+    histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+    centname = centtitle = "MinBias";
+    hasCentralityInfo = kTRUE;
+  }
   histo->GetAxis(3)->SetRange(step + 1, step + 1); 
 
   TH2 *hSpec = histo->Projection(2, 1);
@@ -267,19 +306,20 @@ TH2 *AliHFEtrdPIDqaV1::MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t
   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
-  TString histname = Form("hSigmaTPC%s%s", specID.Data(), stepname.Data());
-  TString histtitle = Form("TPC Sigma for %s %s PID", speciesname.Data(), stepname.Data());
+  TString histname = Form("hSigmaTPC%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+  TString histtitle = Form("TPC Sigma for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
   hSpec->SetName(histname.Data());
   hSpec->SetTitle(histtitle.Data());
 
   // Unset range on the original histogram
   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
-  histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+  histo->GetAxis(3)->SetRange(0, histo->GetAxis(3)->GetNbins());
+  if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
   return hSpec; 
 }
 
 //_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
   //
   // Get the TPC control histogram for the TRD selection step (either before or after PID)
   //
@@ -292,6 +332,25 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t spe
     // cut on species (if available)
     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+     if(histo->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo= kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form("Centrality %d", centralityClass);
+    }
+  } else {
+    histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+    centname = centtitle = "MinBias";
+    hasCentralityInfo= kFALSE;
+  }
   histo->GetAxis(3)->SetRange(step + 1, step + 1); 
 
   TH2 *hSpec = histo->Projection(2, 1);
@@ -299,19 +358,20 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t spe
   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
-  TString histname = Form("hTMTRD%s%s", specID.Data(), stepname.Data());
-  TString histtitle = Form("TRD Truncated Mean for %s %s PID", speciesname.Data(), stepname.Data());
+  TString histname = Form("hTMTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+  TString histtitle = Form("TRD Truncated Mean for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
   hSpec->SetName(histname.Data());
   hSpec->SetTitle(histtitle.Data());
 
   // Unset range on the original histogram
   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+  if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
   return hSpec; 
 }
 
 //_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
   //
   // Make Histogram for TRD Likelihood distribution
   //
@@ -324,6 +384,25 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t ste
     // cut on species (if available)
     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+    if(histo->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo= kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form("Centrality %d", centralityClass);
+    }
+  } else {
+    histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+    centname = centtitle = "MinBias";
+    hasCentralityInfo= kTRUE;
+  }
   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 
   TH2 *hSpec = histo->Projection(2, 1);
@@ -331,19 +410,20 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t ste
   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
-  TString histname = Form("hLikeElTRD%s%s", specID.Data(), stepname.Data());
-  TString histtitle = Form("TRD electron Likelihood for %s %s PID", speciesname.Data(), stepname.Data());
+  TString histname = Form("hLikeElTRD%s%s%s", specID.Data(), stepname.Data(),centname.Data());
+  TString histtitle = Form("TRD electron Likelihood for %s %s PID %s", speciesname.Data(), stepname.Data(),centtitle.Data());
   hSpec->SetName(histname.Data());
   hSpec->SetTitle(histtitle.Data());
 
   // Unset range on the original histogram
   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+  if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
   return hSpec; 
 }
 
 //_________________________________________________________
-TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species){
+TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species, Int_t centralityClass){
   //
   // Make Histogram for TRD Likelihood distribution
   //
@@ -356,6 +436,25 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, I
     // cut on species (if available)
     histo->GetAxis(0)->SetRange(species + 2, species + 2); // undef + underflow
   }
+  TString centname, centtitle;
+  Bool_t hasCentralityInfo = kTRUE;
+  if(centralityClass > -1){
+    if(histo->GetNdimensions() < 5){
+      AliError("Centrality Information not available");
+      centname = centtitle = "MinBias";
+      hasCentralityInfo= kFALSE;
+    } else {
+      // Project centrality classes
+      // -1 is Min. Bias
+      histo->GetAxis(4)->SetRange(centralityClass+1, centralityClass+1);
+      centname = Form("Cent%d", centralityClass);
+      centtitle = Form("Centrality %d", centralityClass);
+    }
+  } else {
+    histo->GetAxis(4)->SetRange(1, histo->GetAxis(4)->GetNbins()-1);
+    centname = centtitle = "MinBias";
+    hasCentralityInfo= kTRUE;
+  }
   histo->GetAxis(3)->SetRange(step + 1, step + 1);
 
   TH2 *hSpec = histo->Projection(2, 1);
@@ -363,14 +462,15 @@ TH2 *AliHFEtrdPIDqaV1::MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, I
   TString stepname = step == AliHFEdetPIDqa::kBeforePID ? "before" : "after";
   TString speciesname = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "all Particles";
   TString specID = species > -1 && species < AliPID::kSPECIES ? AliPID::ParticleName(species) : "unid";
-  TString histname = Form("hChargeTRD%s%s", specID.Data(), stepname.Data());
-  TString histtitle = Form("TRD total charge for %s %s PID", speciesname.Data(), stepname.Data());
+  TString histname = Form("hChargeTRD%s%s%s", specID.Data(), stepname.Data(), centname.Data());
+  TString histtitle = Form("TRD total charge for %s %s PID %s", speciesname.Data(), stepname.Data(), centtitle.Data());
   hSpec->SetName(histname.Data());
   hSpec->SetTitle(histtitle.Data());
 
   // Unset range on the original histogram
   histo->GetAxis(0)->SetRange(0, histo->GetAxis(0)->GetNbins());
   histo->GetAxis(2)->SetRange(0, histo->GetAxis(2)->GetNbins());
+  if(hasCentralityInfo)histo->GetAxis(4)->SetRange(0, histo->GetAxis(4)->GetNbins());
   return hSpec; 
 }
 
index 06f08a0..922c383 100644 (file)
@@ -45,10 +45,12 @@ class AliHFEtrdPIDqaV1 : public AliHFEdetPIDqa{
     virtual void Initialize();
     virtual void ProcessTrack(const AliHFEpidObject *track, AliHFEdetPIDqa::EStep_t step);
 
-    TH2 *MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
-    TH2 *MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
-    TH2 *MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
-    TH2 *MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1);
+    AliHFEcollection *GetListOfHistograms() const { return fHistos; }
+
+    TH2 *MakeTPCspectrumNsigma(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+    TH2 *MakeTRDspectrumTM(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+    TH2 *MakeTRDlikelihoodDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
+    TH2 *MakeTRDchargeDistribution(AliHFEdetPIDqa::EStep_t step, Int_t species = -1, Int_t centralityClass = -1);
   protected:
     void ProcessESDtrack(const AliESDtrack *esdtrack, AliHFEdetPIDqa::EStep_t step, Int_t species);
     void ProcessAODtrack(const AliAODTrack *aodtrack, AliHFEdetPIDqa::EStep_t step, Int_t species);