Modifications to p-Pb analysis from Hongsheng Zhu
authorYuri Kharlov <Yuri.Kharlov@cern.ch>
Thu, 16 Jan 2014 15:25:10 +0000 (19:25 +0400)
committerYuri Kharlov <Yuri.Kharlov@cern.ch>
Thu, 16 Jan 2014 15:25:10 +0000 (19:25 +0400)
PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx
PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.h
PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.cxx
PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.h
PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx
PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.h
PWGGA/PHOSTasks/UserTasks/macros/AddTaskPHOSpPbPi0.C

index 07f5f8c..2c16a20 100644 (file)
 
 ClassImp(AliAnalysisTaskSEPHOSpPbPi0)
 
-const Float_t AliAnalysisTaskSEPHOSpPbPi0::kLogWeight = 4.5;
 Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgRemovePileup    = kFALSE;
 Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgUseFiducialCut  = kFALSE;
-Double_t AliAnalysisTaskSEPHOSpPbPi0::fgDecaliWidth     = 0.055;
-Double_t AliAnalysisTaskSEPHOSpPbPi0::fgCuts[5]         = { 0.3, 2., 0.2, 2.5, 7e-8 };
+Bool_t   AliAnalysisTaskSEPHOSpPbPi0::fgUseTOFCut       = kFALSE;
+Double_t AliAnalysisTaskSEPHOSpPbPi0::fgCuts[5]         = { 0.3, 2., 0.2, 2.5, 1e-7 };
 
 //________________________________________________________________________
 AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0():
   AliAnalysisTaskSE(), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
-  fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
+  fPHOSGeo(0), fOutputListQA(0), fOutputListRD(0), fOutputListMC(0), fHeader(0), fCaloClArr(0)
 {
   //
   // Default constructor
   //
   for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
-
-  // Init bad channel map
-  for (Int_t i=0; i<5; i++)
-    fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
 }
 //________________________________________________________________________
 AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0(const char *name):
   AliAnalysisTaskSE(name), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10), fRunNumber(-1),
-  fPHOSGeo(0), fPHOSCalibData(0), fListQA(0), fListRD(0), fListMC(0), fHeader(0), fCaloClArr(0)
+  fPHOSGeo(0), fOutputListQA(0), fOutputListRD(0), fOutputListMC(0), fHeader(0), fCaloClArr(0)
 {
   // Constructor
   for (Int_t i=0; i<10; i++) { for (Int_t j=0; j<10; j++) fEventList[i][j] = 0; }
 
-  // Init bad channel map
-  for (Int_t i=0; i<5; i++)
-    fPHOSBadMap[i] = new TH2I(Form("PHOS_BadMap_mod%d", i), Form("PHOS_BadMap_mod%d", i), 64, 0., 64., 56, 0., 56.);
-
   DefineOutput(1, TList::Class());
   DefineOutput(2, TList::Class());
   DefineOutput(3, TList::Class());
@@ -100,12 +91,11 @@ AliAnalysisTaskSEPHOSpPbPi0::~AliAnalysisTaskSEPHOSpPbPi0()
   //
   // Default destructor
   //
-  if (fListQA)      { delete fListQA;      fListQA     = NULL; }
-  if (fListRD)      { delete fListRD;      fListRD     = NULL; }
-  if (fListMC)      { delete fListMC;      fListMC     = NULL; }
-  if (fHeader)      { delete fHeader;      fHeader     = NULL; }
-  if (fCaloClArr)   { delete fCaloClArr;   fCaloClArr  = NULL; }
-
+  if (fOutputListQA)  { delete fOutputListQA;  fOutputListQA = NULL; }
+  if (fOutputListRD)  { delete fOutputListRD;  fOutputListRD = NULL; }
+  if (fOutputListMC)  { delete fOutputListMC;  fOutputListMC = NULL; }
+  if (fHeader)        { delete fHeader;        fHeader       = NULL; }
+  if (fCaloClArr)     { delete fCaloClArr;     fCaloClArr    = NULL; }
 }
 //________________________________________________________________________
 void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
@@ -115,19 +105,19 @@ void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
   AliPHOSpPbPi0Header::SetIsMC(fIsMC);
   AliPHOSpPbPi0Header::SetUseFiducialCut(fgUseFiducialCut);
 
-  if (!fHeader)     fHeader     = new AliPHOSpPbPi0Header();
-  if (!fCaloClArr)  fCaloClArr  = new TClonesArray("AliCaloClusterInfo", 0);
-  if (!fListQA)     fListQA     = new TList();
-  if (!fListRD)     fListRD     = new TList();
-  if (!fListMC)     fListMC     = new TList();
+  if (!fHeader)       fHeader       = new AliPHOSpPbPi0Header();
+  if (!fCaloClArr)    fCaloClArr    = new TClonesArray("AliCaloClusterInfo", 0);
+  if (!fOutputListQA) fOutputListQA = new TList();   fOutputListQA->SetOwner(kTRUE);
+  if (!fOutputListRD) fOutputListRD = new TList();   fOutputListRD->SetOwner(kTRUE);
+  if (!fOutputListMC) fOutputListMC = new TList();   fOutputListMC->SetOwner(kTRUE);
 
   fHeader->SetNCent(fCentralityBin.GetSize()-1);
-  fHeader->CreateHistograms(fListQA, fListRD, fListMC);
+  fHeader->CreateHistograms(fOutputListQA, fOutputListRD, fOutputListMC);
 
   // Post output data.
-  PostData(1, fListQA);
-  PostData(2, fListRD);
-  if (fIsMC) PostData(3, fListMC);
+  PostData(1, fOutputListQA);
+  PostData(2, fOutputListRD);
+  if (fIsMC) PostData(3, fOutputListMC);
 
   return;
 }
@@ -160,24 +150,23 @@ void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *)
 
   // Fill Event info
   fHeader->SetEventInfo(fInputHandler);
-  fHeader->FillHistosEvent(fListQA);
+  fHeader->FillHistosEvent(fOutputListQA);
 
   // PHOS Geometry and Misalignment initialization at the first time it runs
   if(fRunNumber != fInputEvent->GetRunNumber()) {
     fRunNumber = fInputEvent->GetRunNumber();
-    fHeader->SetIspARun(fRunNumber>195344 && fRunNumber<197388);  // flag for pA collisions
     PHOSInitialize(esd);
   }
 
   // Event Selection
   if (!fHeader->IsSelected())                     return;
-  if (fgRemovePileup && fHeader->IsPileupSPD())   return;
+  if (fgRemovePileup && fHeader->IsPileup())      return;
 
   // Fill PHOS cells QA histograms
-  fHeader->FillHistosCaloCellsQA(fListQA, fInputEvent->GetPHOSCells(), fPHOSGeo);
+  fHeader->FillHistosCaloCellsQA(fOutputListQA, fInputEvent->GetPHOSCells(), fPHOSGeo);
 
   // Fill PHOS cluster Clones Array
-  FillCaloClusterInfo(aod, esd);
+  FillCaloClusterInfo(/*aod, esd*/);
   aod = 0x0;
 
   if (!fCaloClArr->GetEntriesFast()) return;
@@ -189,25 +178,22 @@ void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *)
   TList *eventList = fEventList[zvtx][cent];
 
   // Fill cluster histograms
-  fHeader->FillHistosCaloCluster(fListQA, fCaloClArr, cent);
+  fHeader->FillHistosCaloCluster(fOutputListQA, fCaloClArr, cent);
 
   // Fill pi0 histograms
-  fHeader->FillHistosPi0(fListRD, fCaloClArr, cent);
+  fHeader->FillHistosPi0(fOutputListRD, fCaloClArr, cent);
 
   // Fill mixed pi0 histograms
-  fHeader->FillHistosMixPi0(fListRD, fCaloClArr, eventList, cent);
+  fHeader->FillHistosMixPi0(fOutputListRD, fCaloClArr, eventList, cent);
 
-  // Fill MC info
+  // Fill MC info JUST FOR ESD
   AliStack *stack = 0x0;
   if (fIsMC) {
-    if (esd) {
-      if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
-        if (static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
-          stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
-      }
-      fHeader->FillHistosMC(fListMC, stack, fPHOSGeo, cent);
-    } else
-      fHeader->FillHistosMC(fListMC, MCEvent(), fPHOSGeo, cent);
+    if (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()) {
+      if (static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent())
+        stack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
+    }
+    fHeader->FillHistosMC(fOutputListMC, stack, fCaloClArr, fPHOSGeo, cent);
   }
   esd = 0x0;
 
@@ -238,7 +224,6 @@ void AliAnalysisTaskSEPHOSpPbPi0::Terminate(Option_t *)
 void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
 {
   // Initialize PHOS Geometry ,misalignment and calibration 
-//TGeoManager::Import("$ALICE_ROOT/test/QA/geometry.root");
 
   fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP");
   AliOADBContainer geomContainer("phosGeo");
@@ -257,21 +242,10 @@ void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
       else fPHOSGeo->SetMisalMatrix(modMatrix, mod);
     }
   }
-
-  TRandom random;   random.SetSeed(0); // the seed is set to the current  machine clock
-  fPHOSCalibData = new AliPHOSCalibData();   fPHOSCalibData->SetName("PHOSCalibData");
-  for(Int_t iMod=1; iMod<6; iMod++) {
-    for(Int_t iCol=1; iCol<57; iCol++) {
-      for(Int_t iRow=1; iRow<65; iRow++) {
-        Float_t value = fIsMC ? random.Gaus(1., fgDecaliWidth) : 0.;  // ADC channel Emc
-        fPHOSCalibData->SetADCchannelEmc(iMod, iCol, iRow, value);
-      }
-    }
-  }
 }
 
 //________________________________________________________________________
-void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, AliESDEvent* const esd)
+void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(/*AliAODEvent* const aod, AliESDEvent* const esd*/)
 {
   // Fill calo cluster info
   if (fCaloClArr) fCaloClArr->Clear();
@@ -282,8 +256,6 @@ void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, Al
   Int_t    relID[4]     = {0,0,0,0 }; // module = relID[0]; cellX = relID[2]; cellZ = relID[3];
   Float_t  position[3]  = {0.,0.,0.};
   Double_t vtx[3]       = {0,0,0};   fHeader->GetXYZ(vtx);   TVector3 vtxVector(vtx); 
-  Double_t magfield     = fInputEvent->GetMagneticField();
-  TF1 *nonLinCorr = new TF1("Non-linear", "0.0241+1.0504*x+0.000249*x*x", 0., 50.); // GCB's non-linear correction function
 
   TClonesArray &caloRef = *fCaloClArr;
   TLorentzVector momentum;
@@ -294,37 +266,16 @@ void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, Al
     if (!(clust && clust->IsPHOS() && clust->E()>fgCuts[0]))                 { clust=0; continue; }
     if (!(clust->GetNCells()>(Int_t)fgCuts[1] && clust->GetM02()>fgCuts[2])) { clust=0; continue; } // To remove exotic clusters
     if (!(clust->GetDistanceToBadChannel()>fgCuts[3]))                       { clust=0; continue; }
-//  if (!(TMath::Abs(clust->GetTOF()<fgCuts[4])))                            { clust=0; continue; } // remove clusters from pileup or same benches
+    if (fgUseTOFCut && !(TMath::Abs(clust->GetTOF())<fgCuts[4]))             { clust=0; continue; } // remove clusters from pileup or same benches
     clust->GetPosition(position); TVector3 global(position); fPHOSGeo->GlobalPos2RelId(global,relID);
     if (relID[0] == 2)                                                       { clust=0; continue; } // !remove module 2
-    if (!IsGoodCaloCluster(relID[0], relID[2], relID[3]))                    { clust=0; continue; } // calo cluster selection 
-
-    caloCluster = new AliCaloClusterInfo(clust, esd, relID, magfield);
-    if (fgUseFiducialCut && !caloCluster->IsInFiducialRegion(relID[2], relID[3])) { delete caloCluster; caloCluster=0; clust=0; continue; } // Fiducial cut 
-
-    if (aod) { // TODO recalibration for AOD
-      if (fIsMC) {
-        AliPHOSAodCluster aodClust(*(AliAODCaloCluster*) (clust));
-        aodClust.Recalibrate(fPHOSCalibData, aod->GetPHOSCells());
-        aodClust.EvalAll(kLogWeight, vtxVector);        // recalculate all cluster parameters
-        aodClust.SetE(nonLinCorr->Eval(aodClust.E()));  // non-linear correction
-        aodClust.GetMomentum(momentum, vtx);
-      } else clust->GetMomentum(momentum, vtx);
-      caloCluster->SetLorentzVector(momentum);
-    } else { // for ESD
-      AliPHOSEsdCluster esdClust(*(AliESDCaloCluster*) (clust));
-      esdClust.Recalibrate(fPHOSCalibData, esd->GetPHOSCells());
-      esdClust.EvalAll(kLogWeight, vtxVector);       // recalculate all cluster parameters
-      esdClust.SetE(nonLinCorr->Eval(esdClust.E()));  // non-linear correction
-      esdClust.GetMomentum(momentum, vtx);
-      caloCluster->SetLorentzVector(momentum);
-    }
 
-    if (fIsMC && !(momentum.E()>fgCuts[0])) { delete caloCluster; caloCluster=0; clust=0; continue; } // check for MC
+    caloCluster = new AliCaloClusterInfo(clust, relID);
+//  if (fgUseFiducialCut && !caloCluster->IsInFiducialRegion(relID[2], relID[3])) { delete caloCluster; caloCluster=0; clust=0; continue; } // Fiducial cut 
 
-    Double_t disp = caloCluster->TestDisp();
-    if (disp < 2.5*2.5) caloCluster->SetPIDBit(BIT(1)); // set Disp1 bit
-    if (disp < 1.5*1.5) caloCluster->SetPIDBit(BIT(3)); // set Disp2 bit
+    clust->GetMomentum(momentum, vtx);
+    caloCluster->SetLorentzVector(momentum);
+    if (fIsMC) caloCluster->SetLabels(clust->GetNLabels(), clust->GetLabels());
 
     clust = 0;
 
@@ -332,14 +283,3 @@ void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(AliAODEvent* const aod, Al
     delete caloCluster; caloCluster=0;
   }  // end loop of all clusters
 }
-
-//________________________________________________________________________
-Bool_t AliAnalysisTaskSEPHOSpPbPi0::IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ)
-{
-  // !Check whether this cluster is not in bad channel
-
-  if (!(iMod>0 && iMod<5 && fPHOSBadMap[iMod]))          return kTRUE;   //No bad maps for this Module
-  if (fPHOSBadMap[iMod]->GetBinContent(cellX,cellZ)>0)   return kFALSE;
-
-  return kTRUE;
-}
index 7a50f08..c4b9935 100644 (file)
@@ -1,5 +1,5 @@
 #ifndef ALIANALYSISTASKSEPHOSPPBPI0_H
-#define ALIANALYSISTASKSEPHOSPPBPI0_H
+#define ALIANALYSISTAKSSEPHOSPPBPI0_H
 
 /* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
  * See cxx source for full Copyright notice                               */
@@ -21,7 +21,6 @@ class TClonesArray;
 class AliAODEvent;
 class AliESDEvent;
 class AliPHOSGeoUtils;
-class AliPHOSCalibData;
 class AliPHOSpPbPi0Header;
 
 class AliAnalysisTaskSEPHOSpPbPi0 : public AliAnalysisTaskSE {
@@ -35,20 +34,13 @@ class AliAnalysisTaskSEPHOSpPbPi0 : public AliAnalysisTaskSE {
   virtual void UserExec(Option_t *option);
   virtual void Terminate(Option_t *option);
 
-  void SetPHOSBadMap(Int_t mod,TH2I *hMap)
-  {
-    if(fPHOSBadMap[mod]) delete fPHOSBadMap[mod];
-    fPHOSBadMap[mod]=new TH2I(*hMap);
-    printf("Set %s \n",fPHOSBadMap[mod]->GetName());
-  }
-
   void SetUseMC(Bool_t isMC=kFALSE)                           { fIsMC          = isMC;                         }
   void SetXBins(const TArrayF& tCent, const TArrayI& tBuffer) { fCentralityBin = tCent; fBufferSize = tBuffer; }
-  void SetEventCuts(Double_t cuts[4])                   const { AliPHOSpPbPi0Header::SetSelectionCuts(cuts);   }
+  void SetEventCuts(Double_t cuts[3])                   const { AliPHOSpPbPi0Header::SetSelectionCuts(cuts);   }
 
   static void SetRemovePileup(Bool_t rm=kFALSE)               { fgRemovePileup                   = rm;         }
   static void SetUseFiducialCut(Bool_t fc=kFALSE)             { fgUseFiducialCut                 = fc;         }
-  static void SetDecaliWidth(Double_t width)                  { fgDecaliWidth                    = width;      }
+  static void SetUseTOFCut(Bool_t tof=kFALSE)                 { fgUseTOFCut                      = tof;        }
   static void SetCaloClCuts(Double_t cuts[5])                 { for (Int_t i=5; i--; ) fgCuts[i] = cuts[i];    }          
 
  private:
@@ -57,15 +49,11 @@ class AliAnalysisTaskSEPHOSpPbPi0 : public AliAnalysisTaskSE {
   AliAnalysisTaskSEPHOSpPbPi0& operator=(const AliAnalysisTaskSEPHOSpPbPi0&); // not implemented
 
   void PHOSInitialize(AliESDEvent* const esd);
-  void FillCaloClusterInfo(AliAODEvent* const aod, AliESDEvent* const esd);
-
-  Bool_t IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ);
-
-  static const Float_t kLogWeight;
+  void FillCaloClusterInfo(/*AliAODEvent* const aod, AliESDEvent* const esd*/);
 
   static Bool_t        fgRemovePileup;      // flag of remove pileup events
   static Bool_t        fgUseFiducialCut;    // flag of use fiducial cut
-  static Double_t      fgDecaliWidth;       // decalibration width
+  static Bool_t        fgUseTOFCut;         // flag of use cluster TOF cut
   static Double_t      fgCuts[5];           // 0, min of cluster Energy
                                             // 1, min of NCells
                                             // 2, min of M02
@@ -77,18 +65,16 @@ class AliAnalysisTaskSEPHOSpPbPi0 : public AliAnalysisTaskSE {
   TArrayI              fBufferSize;         // Buffer size for event mixing
 
   Int_t                fRunNumber;          // Run Number
-  TH2I                *fPHOSBadMap[5];      // Container of PHOS bad channels map
   AliPHOSGeoUtils     *fPHOSGeo;            // PHOS geometry
-  AliPHOSCalibData    *fPHOSCalibData;      // PHOS calibration
 
   TList               *fEventList[10][10];  // Event list for mixing
-  TList               *fListQA;             // output list of QA histograms
-  TList               *fListRD;             // output list of Real Data histograms
-  TList               *fListMC;             // output list of MC histograms
+  TList               *fOutputListQA;       // output list of QA histograms
+  TList               *fOutputListRD;       // output list of RD histograms
+  TList               *fOutputListMC;       // output list of MC histograms
   AliPHOSpPbPi0Header *fHeader;             // info at event level
   TClonesArray        *fCaloClArr;          // Container of Calo clusters Info
    
-  ClassDef(AliAnalysisTaskSEPHOSpPbPi0, 1);
+  ClassDef(AliAnalysisTaskSEPHOSpPbPi0, 2);
 };
 
 #endif
index 556d8e4..e77fad4 100644 (file)
@@ -22,7 +22,9 @@
 /////////////////////////////////////////////////////////////
 
 #include <iostream>
+#include <TParticle.h>
 
+#include "AliStack.h"
 #include "AliESDEvent.h"
 #include "AliAODTrack.h"
 #include "AliESDtrack.h"
@@ -38,7 +40,7 @@ ClassImp(AliCaloClusterInfo)
 //-------------------------------------------------------------------------------------------
 AliCaloClusterInfo::AliCaloClusterInfo():
   TObject(), fLorentzVector(), fModule(0), fTRUNumber(0), fNCells(0), fPIDBit(0x0),
-  fDistToBad(0.), fEmcCpvDistance(0.), fM02(0.), fM20(0.), fTOF(0.)
+  fLabels(0), fDistToBad(0.), fM02(0.), fM20(0.), fTOF(0.)
 {
   //
   // Default constructor
@@ -46,20 +48,20 @@ AliCaloClusterInfo::AliCaloClusterInfo():
 } 
 
 //-------------------------------------------------------------------------------------------
-AliCaloClusterInfo::AliCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, Int_t relID[4], Double_t mf):
+AliCaloClusterInfo::AliCaloClusterInfo(AliVCluster* const clust, Int_t relID[4]):
   TObject(), fLorentzVector(), fModule(0), fTRUNumber(0), fNCells(0), fPIDBit(0x0),
-  fDistToBad(0.), fEmcCpvDistance(0.), fM02(0.), fM20(0.), fTOF(0.)
+  fLabels(0), fDistToBad(0.), fM02(0.), fM20(0.), fTOF(0.)
 {
   //
   // constructor
   //
-  this->FillCaloClusterInfo(clust, esd, relID, mf);
+  this->FillCaloClusterInfo(clust, relID);
 } 
 
 //-------------------------------------------------------------------------------------------
 AliCaloClusterInfo::AliCaloClusterInfo(const AliCaloClusterInfo &src):
   TObject(src), fLorentzVector(src.fLorentzVector), fModule(src.fModule), fTRUNumber(src.fTRUNumber),fNCells(src.fNCells), fPIDBit(src.fPIDBit),
-  fDistToBad(src.fDistToBad), fEmcCpvDistance(src.fEmcCpvDistance), fM02(src.fM02), fM20(src.fM20), fTOF(src.fTOF)
+  fLabels(src.fLabels), fDistToBad(src.fDistToBad), fM02(src.fM02), fM20(src.fM20), fTOF(src.fTOF)
 {
   //
   // copy constructor
@@ -79,8 +81,8 @@ AliCaloClusterInfo& AliCaloClusterInfo::operator=(const AliCaloClusterInfo &src)
   fTRUNumber       = src.fTRUNumber;
   fNCells          = src.fNCells;
   fPIDBit          = src.fPIDBit;
+  fLabels          = src.fLabels;
   fDistToBad       = src.fDistToBad;
-  fEmcCpvDistance  = src.fEmcCpvDistance;
   fM02             = src.fM02;
   fM20             = src.fM20;
   fTOF             = src.fTOF;
@@ -97,82 +99,50 @@ AliCaloClusterInfo::~AliCaloClusterInfo()
 }
 
 //-----------------------------------------------------------------------------
-void AliCaloClusterInfo::FillCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, Int_t relID[4], Double_t mf)
+void AliCaloClusterInfo::FillCaloClusterInfo(AliVCluster* const clust, Int_t relID[4])
 {
   // extract information of calo clusters
 
-  Short_t  trkCharge = 0;
-  Double_t trkDz = 0., trkDx = 0., trkPt = -1.;
-
-  if (!esd) { // for AOD
-    AliAODTrack *trkAOD = 0x0;
-    if (clust->GetNTracksMatched()>0) {
-      trkAOD = dynamic_cast <AliAODTrack*> (clust->GetTrackMatched(0));
-      if (trkAOD) {
-        trkPt     = trkAOD->Pt();
-        trkCharge = trkAOD->Charge();
-      }
-    }
-  } else { // for ESD
-    Int_t iESDtrack = clust->GetTrackMatchedIndex();
-    AliESDtrack *trkESD = 0x0;
-    if (iESDtrack>-1) { 
-      trkESD = esd->GetTrack(iESDtrack);
-      if (trkESD) {
-        trkPt     = trkESD->Pt();
-        trkCharge = trkESD->Charge();
-      }
-    }
-  }
-  trkDz = clust->GetTrackDz();
-  trkDx = clust->GetTrackDx();
-
   fModule         = relID[0];
   fTRUNumber      = GetTRUNumber(relID[2], relID[3]);
   fNCells         = clust->GetNCells();
   fDistToBad      = clust->GetDistanceToBadChannel();
-  fEmcCpvDistance = clust->GetEmcCpvDistance();
   fM02            = clust->GetM02();
   fM20            = clust->GetM20();
   fTOF            = clust->GetTOF();
 
-  Double_t cpv = TestCpv(trkPt, trkCharge, trkDz, trkDx, mf);
-  if (trkPt == -1. || cpv>2.) this->SetPIDBit(BIT(0));
-  if (trkPt == -1. || cpv>4.) this->SetPIDBit(BIT(2));
+  Double_t cpv  = clust->GetEmcCpvDistance(); // Distance in sigmas filled by tender
+  Double_t disp = TestDisp();                 // Dispersion in sigmas filled by tender
+  if (cpv  > 2.)      this->SetPIDBit(BIT(0));
+  if (cpv  > 4.)      this->SetPIDBit(BIT(2));
+  if (disp < 2.5*2.5) this->SetPIDBit(BIT(1));
+  if (disp > 1.5*1.5) this->SetPIDBit(BIT(3));
 
   return;
 }
 
 //-----------------------------------------------------------------------------
-Double_t AliCaloClusterInfo::TestCpv(Double_t trkPt, Short_t trkCharge, Double_t trkDz, Double_t trkDx, Double_t mf)
+Bool_t AliCaloClusterInfo::IsInFiducialRegion(Int_t cellX, Int_t cellZ)
 {
-  // Parameterization of LHC10h period
-
-  Double_t meanX=0;
-  Double_t meanZ=0.;
-  Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-trkPt/1.02053e-01)+
-                           6.58365e-01*5.91917e-01*5.91917e-01/((trkPt-9.61306e-01)*(trkPt-9.61306e-01)+5.91917e-01*5.91917e-01)+1.59219);
-  Double_t sz = TMath::Min(2.75,4.90341e+02*1.91456e-02*1.91456e-02/(trkPt*trkPt+1.91456e-02*1.91456e-02)+1.60);
-
-  if(mf<0.) { // Field --
-    meanZ = -0.468318;
-    if (trkCharge>0)
-      meanX = TMath::Min(7.3, 3.89994*1.20679*1.20679/(trkPt*trkPt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-trkPt*3.33650e+01));
-    else
-      meanX =-TMath::Min(7.7,3.86040*0.912499*0.912499/(trkPt*trkPt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-trkPt*2.57070e+01));
-  }
-  else {      // Field ++
-    meanZ = -0.468318;
-    if (trkCharge>0)
-      meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(trkPt*trkPt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-trkPt*3.08451e+01));
-    else
-      meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(trkPt*trkPt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-trkPt*2.40913e+01));
-  }
+  const Int_t edgeX = 2;
+  const Int_t edgeZ = 2;
+  if (cellX >edgeX && cellX<(65-edgeX) && cellZ>edgeZ && cellZ <(57-edgeZ)) return kTRUE;
 
-  Double_t rx = (trkDx-meanX)/sx;
-  Double_t rz = (trkDz-meanZ)/sz;
+  return kFALSE;
+}
 
-  return TMath::Sqrt(rx*rx+rz*rz);
+//-----------------------------------------------------------------------------
+Bool_t AliCaloClusterInfo::CheckIsClusterFromPi0(AliStack* const stack, Int_t &pi0Indx)
+{
+  if (GetLabel()>-1) {
+    TParticle* gamma = stack->Particle(GetLabel());
+    if (gamma->GetPdgCode() == 22) {
+      pi0Indx = gamma->GetFirstMother();
+      if (pi0Indx>-1 && ((TParticle*)stack->Particle(pi0Indx))->GetPdgCode() == 111) 
+        return kTRUE;
+      else return kFALSE;
+    } else return kFALSE; 
+  } else return kFALSE;
 }
 
 //-----------------------------------------------------------------------------
@@ -195,16 +165,6 @@ Double_t AliCaloClusterInfo::TestDisp()
 }
 
 //-----------------------------------------------------------------------------
-Bool_t AliCaloClusterInfo::IsInFiducialRegion(Int_t cellX, Int_t cellZ)
-{
-  const Int_t edgeX = 2;
-  const Int_t edgeZ = 2;
-  if (cellX >edgeX && cellX<(65-edgeX) && cellZ>edgeZ && cellZ <(57-edgeZ)) return kTRUE;
-
-  return kFALSE;
-}
-
-//-----------------------------------------------------------------------------
 Int_t AliCaloClusterInfo::GetTRUNumber(Int_t cellX, Int_t cellZ)
 {
   // Return TRU region number for given cell.
index a41fa17..4f58127 100644 (file)
 
 #include <TObject.h>
 #include <TString.h>
+#include <TArrayI.h>
 #include <TLorentzVector.h>
 
+class AliStack;
 class AliESDEvent;
 class AliVCluster;
 
@@ -22,35 +24,39 @@ class AliCaloClusterInfo : public TObject{
  public:
  
   AliCaloClusterInfo();
-  AliCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, Int_t relID[4], Double_t mf);
+  AliCaloClusterInfo(AliVCluster* const clust, Int_t relID[4]);
   AliCaloClusterInfo(const AliCaloClusterInfo &src);
   AliCaloClusterInfo& operator=(const AliCaloClusterInfo &src);
   virtual ~AliCaloClusterInfo();
 
   TLorentzVector LorentzVector() const { return fLorentzVector; }
 
-  Int_t    GetModule()         const { return fModule;          }
-  Int_t    GetTRUNumber()      const { return fTRUNumber;       }
-  Int_t    GetNCells()         const { return fNCells;          }
-  UInt_t   GetPIDBit()         const { return fPIDBit;          }
-  Double_t GetDistToBad()      const { return fDistToBad;       }
-  Double_t GetEmcCpvDistance() const { return fEmcCpvDistance;  }
-  Double_t GetM02()            const { return fM02;             }
-  Double_t GetM20()            const { return fM20;             }
-  Double_t GetTOF()            const { return fTOF;             }
+  Int_t    GetModule()          const { return fModule;          }
+  Int_t    GetTRUNumber()       const { return fTRUNumber;       }
+  Int_t    GetNCells()          const { return fNCells;          }
+  TArrayI* GetLabelsArray()     const { return fLabels;          }
+  Int_t    GetLabel()           const { if(fLabels && fLabels->GetSize() >0)        return fLabels->At(0); else return -1;   }
+  Int_t    GetLabelAt(UInt_t i) const { if(fLabels && i<(UInt_t)fLabels->GetSize()) return fLabels->At(i); else return -999; }
+
+  UInt_t   GetPIDBit()          const { return fPIDBit;          }
+  Double_t GetDistToBad()       const { return fDistToBad;       }
+  Double_t GetM02()             const { return fM02;             }
+  Double_t GetM20()             const { return fM20;             }
+  Double_t GetTOF()             const { return fTOF;             }
 
   void SetLorentzVector(TLorentzVector momentum) { fLorentzVector = momentum; }
-  void SetPIDBit(UInt_t bit)                     { fPIDBit       |= bit;      }
+  void SetLabels(UInt_t size, Int_t* labels)     { if(fLabels) delete fLabels; fLabels = new TArrayI(size, labels); }
 
-  Bool_t   IsInFiducialRegion(Int_t cellX, Int_t cellZ);
-  Double_t TestDisp();
+  Bool_t IsInFiducialRegion(Int_t cellX, Int_t cellZ);
+  Bool_t CheckIsClusterFromPi0(AliStack* const stack, Int_t &pi0Indx);
 
  private:
 
-  void FillCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, Int_t relID[4], Double_t mf);
+  void FillCaloClusterInfo(AliVCluster* const clust, Int_t relID[4]);
+  void SetPIDBit(UInt_t bit)   { fPIDBit |= bit; }
+  Double_t TestDisp();
 
-  Int_t    GetTRUNumber(Int_t cellX, Int_t cellZ);
-  Double_t TestCpv(Double_t trkPt, Short_t trkCharge, Double_t trkDz, Double_t trkDx, Double_t mf);
+  Int_t  GetTRUNumber(Int_t cellX, Int_t cellZ);
 
   TLorentzVector fLorentzVector;
 
@@ -58,13 +64,13 @@ class AliCaloClusterInfo : public TObject{
   Int_t    fTRUNumber;          // TRU Number
   Int_t    fNCells;             // Number of cells in cluster
   UInt_t   fPIDBit;             // PID Bit
-  Double_t fDistToBad ;         // Distance to nearest bad channel
-  Double_t fEmcCpvDistance;     // Distance from PHOS EMC rec.point to the closest CPV rec.point
+  TArrayI* fLabels;             // list of primaries that generated the cluster, ordered in deposited energy
+  Double_t fDistToBad;          // Distance to nearest bad channel
   Double_t fM02;                // lambda0
   Double_t fM20;                // lambda1
   Double_t fTOF;
 
-  ClassDef(AliCaloClusterInfo,1);
+  ClassDef(AliCaloClusterInfo, 2);
 };
 
 #endif // #ifdef ALICALOCLUSTERINFO_H
index f249e0e..abf45b9 100644 (file)
 #include "AliAODMCParticle.h"
 #include "AliVEvent.h"
 #include "AliVVertex.h"
+#include "AliAnalysisUtils.h"
 #include "AliVCaloCells.h"
 #include "AliPHOSGeoUtils.h"
-#include "AliAODEvent.h"
-#include "AliESDEvent.h"
 #include "AliCentrality.h"
 
 #include "AliCaloClusterInfo.h"
@@ -49,20 +48,13 @@ class TNamed;
 ClassImp(AliPHOSpPbPi0Header)
 
 Bool_t   AliPHOSpPbPi0Header::fgIsMC           = kFALSE;
-Bool_t   AliPHOSpPbPi0Header::fgIspARun        = kFALSE;
 Bool_t   AliPHOSpPbPi0Header::fgUseFiducialCut = kFALSE;
 Int_t    AliPHOSpPbPi0Header::fgNCent          = 10;
-Double_t AliPHOSpPbPi0Header::fgCuts[4]        = { 1., 10., 0., 100. };
+Double_t AliPHOSpPbPi0Header::fgCuts[3]        = { 10., 0., 100. };
 
 //_____________________________________________________________________________
 AliPHOSpPbPi0Header::AliPHOSpPbPi0Header() :
-TNamed(),
-fFiredTriggerClass(),
-fSelMask(0),
-fVtxContrsN(0),
-fIsVertexOK(kFALSE),
-fIsPileupSPD(kFALSE),
-fCentrality(0.)
+  TNamed(), fFiredTriggerClass(), fSelMask(0), fIsVertexOK(kFALSE), fIsPileup(kFALSE), fCentrality(0.)
 {
   //
   // default constructor
@@ -72,13 +64,8 @@ fCentrality(0.)
 
 //_____________________________________________________________________________
 AliPHOSpPbPi0Header::AliPHOSpPbPi0Header(const AliPHOSpPbPi0Header &src) :
-TNamed(),
-fFiredTriggerClass(src.fFiredTriggerClass),
-fSelMask(src.fSelMask),
-fVtxContrsN(src.fVtxContrsN),
-fIsVertexOK(src.fIsVertexOK),
-fIsPileupSPD(src.fIsPileupSPD),
-fCentrality(src.fCentrality)
+  TNamed(), fFiredTriggerClass(src.fFiredTriggerClass), fSelMask(src.fSelMask), fIsVertexOK(src.fIsVertexOK),
+  fIsPileup(src.fIsPileup), fCentrality(src.fCentrality)
 {
   //
   // copy constructor
@@ -97,9 +84,8 @@ AliPHOSpPbPi0Header& AliPHOSpPbPi0Header::operator=(const AliPHOSpPbPi0Header &s
 
   fFiredTriggerClass            = src.fFiredTriggerClass;
   fSelMask                      = src.fSelMask;
-  fVtxContrsN                   = src.fVtxContrsN;
   fIsVertexOK                   = src.fIsVertexOK;
-  fIsPileupSPD                  = src.fIsPileupSPD;
+  fIsPileup                     = src.fIsPileup;
   fCentrality                   = src.fCentrality;
   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
 
@@ -120,16 +106,13 @@ void AliPHOSpPbPi0Header::SetEventInfo(AliInputEventHandler* const handler)
   // fill info at event level
 
   AliVEvent *event = handler->GetEvent();
-  AliAODEvent *aod = 0x0;   aod = dynamic_cast<AliAODEvent*>(event);
-  AliESDEvent *esd = 0x0;   esd = dynamic_cast<AliESDEvent*>(event);
 
-  fIsVertexOK   = CheckEventVertex(aod, esd);
   fSelMask      = handler->IsEventSelected();
-  fIsPileupSPD  = (aod && !aod->GetTracklets()) ? event->IsPileupFromSPD(3,0.8,3.,2.,5.) :
-                                                  event->IsPileupFromSPDInMultBins(); //TODO not sure!!!
+  fIsVertexOK   = CheckEventVertex(event);
 
   AliCentrality *cent = event->GetCentrality();
-  if (cent) fCentrality = cent->GetCentralityPercentile("V0M");
+  if (cent) fCentrality = cent->GetCentralityPercentile("V0A");
+  else      fCentrality = -1.;
 
 //this->SetTitle(vertex->GetTitle());
   return;
@@ -139,40 +122,27 @@ void AliPHOSpPbPi0Header::SetEventInfo(AliInputEventHandler* const handler)
 Bool_t AliPHOSpPbPi0Header::IsSelected()
 {
   // select event according to the event selection cuts
-  if (fVtxContrsN<(Int_t)fgCuts[0])                    return kFALSE; // num. of vtx contributors cut
-  if (TMath::Abs(fVtx[2])>fgCuts[1])                   return kFALSE; // Vtxz cut
-  if (fCentrality<fgCuts[2] || fCentrality>fgCuts[3])  return kFALSE; // centrality selection
-  if (fgIspARun && !fIsVertexOK)                       return kFALSE; // pA vertex cut
+  if (TMath::Abs(fVtx[2])>fgCuts[0])                   return kFALSE; // Vtxz cut
+  if (fCentrality<fgCuts[1] || fCentrality>fgCuts[2])  return kFALSE; // centrality selection
+  if (!fIsVertexOK)                                    return kFALSE; // pA vertex cut
 
   return kTRUE;
 }
 
 //_____________________________________________________________________________
-Bool_t AliPHOSpPbPi0Header::CheckEventVertex(AliAODEvent* const aod, AliESDEvent* const esd)
+Bool_t AliPHOSpPbPi0Header::CheckEventVertex(AliVEvent* const event)
 {
   // check event vertex
 
-  Bool_t   isAOD = kFALSE;
-  if (aod) isAOD = kTRUE;
-  if (esd) isAOD = kFALSE;
-
   // set event basic info
-  fFiredTriggerClass       = isAOD ? aod->GetFiredTriggerClasses() : esd->GetFiredTriggerClasses();
-  const AliVVertex* trkVtx = isAOD ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertex()) :
-                                     dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertex()) ;
-  trkVtx->GetXYZ(fVtx);   fVtxContrsN  = trkVtx->GetNContributors();
-
-  if (!trkVtx || fVtxContrsN<1)                                                           return kFALSE;
-  if (!fgIspARun)                                                                         return kTRUE;  // Following cuts are for pA vertex
-  if (!((TString)trkVtx->GetTitle()).Contains("VertexerTracks"))                          return kFALSE;
-
-  const AliVVertex* spdVtx = isAOD ? dynamic_cast<const AliVVertex*>(aod->GetPrimaryVertexSPD()) :
-                                     dynamic_cast<const AliVVertex*>(esd->GetPrimaryVertexSPD()) ;
-  Double_t cov[6]={0.};  spdVtx->GetCovarianceMatrix(cov);
+  fFiredTriggerClass       = event->GetFiredTriggerClasses();
+  const AliVVertex* trkVtx = event->GetPrimaryVertex();
+  if (!trkVtx)                               return kFALSE;
+  trkVtx->GetXYZ(fVtx);
 
-  if (spdVtx->GetNContributors()<1)                                                       return kFALSE;
-  if (((TString)spdVtx->GetTitle()).Contains("vertexer:Z") && (TMath::Sqrt(cov[5])>0.25)) return kFALSE; // Double_t zRes = TMath::Sqrt(cov[5]);
-  if (TMath::Abs(spdVtx->GetZ() - fVtx[2])>0.5)                                           return kFALSE;
+  AliAnalysisUtils *utils = new AliAnalysisUtils();
+  if (!utils->IsVertexSelected2013pA(event)) return kFALSE;
+  fIsPileup = utils->IsPileUpEvent(event);
 
   return kTRUE;
 }
@@ -193,58 +163,47 @@ void AliPHOSpPbPi0Header::CreateHistograms(TList *listQA, TList *listRD, TList *
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosEvent(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosEvent(TList *listQA)
 {
   // create histograms at event level
 
-  if (!list) list = new TList();
-  list->SetOwner();
-  Bool_t oldStatus = TH1::AddDirectoryStatus();
-  TH1::AddDirectory(kFALSE);
-
-  const Int_t nhs    = 4;
-  TString tName[nhs] = { "VtxNcontr", "Centrality", "Vz", "Pileup" };
-  Int_t   nbins[nhs] = {       202  ,         110 , 500 ,     500  };
-  Double_t xlow[nhs] = {        -2.5,         -10., -25.,     -25. };
-  Double_t  xup[nhs] = {       199.5,         100.,  25.,      25. };
+  const Int_t nhs    = 3;
+  TString tName[nhs] = { "Centrality", "Vz", "Pileup" };
+  Int_t   nbins[nhs] = {         110 , 500 ,     500  };
+  Double_t xlow[nhs] = {         -10., -25.,     -25. };
+  Double_t  xup[nhs] = {         100.,  25.,      25. };
 
   TString hName;
   TH1D   *histo = 0;
   for (Int_t i=0; i<nhs; i++) {
     hName = Form("hEvent_%s", tName[i].Data());
     histo = new TH1D(hName.Data(), hName.Data(), nbins[i], xlow[i], xup[i]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    histo->Sumw2(); listQA->Add(histo); histo = 0;
 
     hName = Form("hEventSel_%s", tName[i].Data());
     histo = new TH1D(hName.Data(), hName.Data(), nbins[i], xlow[i], xup[i]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    histo->Sumw2(); listQA->Add(histo); histo = 0;
   }
-  TH1::AddDirectory(oldStatus);
 
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosCaloCellsQA(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosCaloCellsQA(TList *listQA)
 {
   // create QA histograms for CaloCells
 
-  if (!list) list = new TList();
-  list->SetOwner();
-  Bool_t oldStatus1 = TH1::AddDirectoryStatus(); TH1::AddDirectory(kFALSE);
-  Bool_t oldStatus2 = TH2::AddDirectoryStatus(); TH2::AddDirectory(kFALSE);
-
   TH1D *histo1       = 0;
   TH2D *histo2       = 0;
   const Int_t nhs    = 3;
   const Int_t nMod   = 3;
   TString hName = "hCaloCellsQA_E";
   histo1 = new TH1D(hName, hName, 300, 0., 30.);
-  histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+  histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
   for (Int_t i=0; i<nMod; i++) {
-    hName = Form("hCaloCellsQA_E_Mod%d", i+1);
+    hName  = Form("hCaloCellsQA_E_Mod%d", i+1);
     histo1 = new TH1D(hName.Data(), hName.Data(), 300, 0., 30.);
-    histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+    histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
   }
 
   TString tName[nhs] = { "CentNCells", "Nxz", "Exz" };
@@ -257,41 +216,33 @@ void AliPHOSpPbPi0Header::CreateHistosCaloCellsQA(TList *list)
 
   for (Int_t i=0; i<nhs; i++) {
     if (i == 0) {
-      hName = Form("hCaloCellsQA_%s", tName[i].Data());
+      hName  = Form("hCaloCellsQA_%s", tName[i].Data());
       histo2 = new TH2D(hName.Data(), hName.Data(), nbinsx[i], xlow[i], xup[i], nbinsy[i], ylow[i], yup[i]);
-      histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+      histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
     } else {
       for (Int_t j=0; j<nMod; j++) {
-        hName = Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), j+1);
+        hName  = Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), j+1);
         histo2 = new TH2D(hName.Data(), hName.Data(), nbinsx[i], xlow[i], xup[i], nbinsy[i], ylow[i], yup[i]);
-        histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+        histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
       }
     }
   }
 
-  TH1::AddDirectory(oldStatus1);
-  TH2::AddDirectory(oldStatus2);
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosCaloCluster(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosCaloCluster(TList *listQA)
 {
   // create histograms for CaloCluster
 
-  if (!list) list = new TList();
-  list->SetOwner();
-
-  Bool_t oldStatus1 = TH1::AddDirectoryStatus();   TH1::AddDirectory(kFALSE);
-  Bool_t oldStatus2 = TH2::AddDirectoryStatus();   TH2::AddDirectory(kFALSE);
-
   const Int_t nMods = 3;
   const Int_t nTRUs = 8;
                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs };   // PID
   TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"           };
                          // { kPtClu, kEtaClu, kPhiClu, kM02Clu, kM20Clu, kTOFClu, kNCellsClu, kNClustersClu, kVarsClu };   // clusters
   TString tName[kVarsClu] = {   "Pt",   "Eta",   "Phi",   "M02",   "M20",   "TOF",   "NCells",   "NClusters"           };
-  Int_t    bins[kVarsClu] = {   500 ,  300   ,   120  ,    200 ,    200 ,    400 ,      1000 ,          20             };
+  Int_t    bins[kVarsClu] = {   500 ,  300   ,   120  ,    200 ,    200 ,    400 ,       200 ,          20             };
   Double_t xMin[kVarsClu] = {     0.,   -0.15,     4.4,      0.,      0.,   -2e-6,         0.,           0.            };
   Double_t xMax[kVarsClu] = {    50.,    0.15,     5.6,      2.,      2.,    2e-6,      1000.,          20.            };
 
@@ -302,79 +253,72 @@ void AliPHOSpPbPi0Header::CreateHistosCaloCluster(TList *list)
 
   hName   = Form("hCaloCluster_%s", tName[kNCellsClu].Data());
   hncells = new TH1I(hName.Data(), hName.Data(), bins[kNCellsClu], (Int_t)xMin[kNCellsClu], (Int_t)xMax[kNCellsClu]);
-  hncells->Sumw2(); list->Add(hncells); hncells = 0;
+  hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 
   hName   = Form("hCaloCluster_%s", tName[kNClustersClu].Data());
   hncells = new TH1I(hName.Data(), hName.Data(), bins[kNClustersClu], (Int_t)xMin[kNClustersClu], (Int_t)xMax[kNClustersClu]);
-  hncells->Sumw2(); list->Add(hncells); hncells = 0;
+  hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 
   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
     hName  = Form("hCaloCluster_%s_%s", tName[kPtClu].Data(), namePID[iPID].Data());
     histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
-    histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+    histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 
     for (Int_t icent=0; icent<fgNCent; icent++) {
       hName  = Form("hCaloCluster_%s_%s_cent%d", tName[kPtClu].Data(), namePID[iPID].Data(), icent);
       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
-      histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+      histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
     }
 
-    hName = Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(), namePID[iPID].Data());
+    hName  = Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(), namePID[iPID].Data());
     histo2 = new TH2D(hName.Data(), hName.Data(), bins[kEtaClu], xMin[kEtaClu], xMax[kEtaClu], bins[kPhiClu], xMin[kPhiClu], xMax[kPhiClu]);
-    histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+    histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
     for (Int_t i=0; i<kVarsClu-3; i++) {
-      hName = Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
+      hName  = Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
       histo2 = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
-      histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+      histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
     }
   }
 
   for (Int_t iMod=0; iMod<nMods; iMod++) {
     hName   = Form("hCaloCluster_%s_Mod%d", tName[kNCellsClu].Data(), iMod+1);
     hncells = new TH1I(hName.Data(), hName.Data(), bins[kNCellsClu], (Int_t)xMin[kNCellsClu], (Int_t)xMax[kNCellsClu]);
-    hncells->Sumw2(); list->Add(hncells); hncells = 0;
+    hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 
     hName   = Form("hCaloCluster_%s_Mod%d", tName[kNClustersClu].Data(), iMod+1);
     hncells = new TH1I(hName.Data(), hName.Data(), bins[kNClustersClu], (Int_t)xMin[kNClustersClu], (Int_t)xMax[kNClustersClu]);
-    hncells->Sumw2(); list->Add(hncells); hncells = 0;
+    hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 
     for (Int_t iTRU=0; iTRU<nTRUs; iTRU++) {
       hName  = Form("hCaloCluster_%s_Mod%d_TRU%d", tName[kPtClu].Data(), iMod+1, iTRU+1);
       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
-      histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+      histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
     }
 
     for (Int_t iPID=0; iPID<kPIDs; iPID++) {
       hName  = Form("hCaloCluster_%s_Mod%d_%s", tName[kPtClu].Data(), iMod+1, namePID[iPID].Data());
       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
-      histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+      histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
     }
   }
 
-  hName = Form("hCaloCluster_%s%s", tName[kEtaClu].Data(), tName[kPhiClu].Data());
+  hName  = Form("hCaloCluster_%s%s", tName[kEtaClu].Data(), tName[kPhiClu].Data());
   histo2 = new TH2D(hName.Data(), hName.Data(), bins[kEtaClu], xMin[kEtaClu], xMax[kEtaClu], bins[kPhiClu], xMin[kPhiClu], xMax[kPhiClu]);
-  histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+  histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
   for (Int_t i=0; i<kVarsClu-3; i++) {
-    hName = Form("hCaloCluster_%s%s", tName[0].Data(), tName[i+1].Data());
+    hName  = Form("hCaloCluster_%s%s", tName[0].Data(), tName[i+1].Data());
     histo2 = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
-    histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+    histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
   }
 
-  TH1::AddDirectory(oldStatus1);
-  TH2::AddDirectory(oldStatus2);
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosPi0(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosPi0(TList *listRD)
 {
   // create histograms for Pi0
 
-  if (!list) list = new TList();
-  list->SetOwner();
-
-  Bool_t oldStatus = TH2::AddDirectoryStatus();   TH2::AddDirectory(kFALSE);
-
   const Int_t nComb = 2;
   TString srcMod[nComb]  =  { "Mod11", "Mod33" }; 
                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs      };  // PID
@@ -385,214 +329,227 @@ void AliPHOSpPbPi0Header::CreateHistosPi0(TList *list)
   Double_t xMin[kVarsPi0] = {     0.,   -0.15,     4.4,      0.,          0.,          0.           };
   Double_t xMax[kVarsPi0] = {    50.,    0.15,     5.6,      1., TMath::Pi(),          1.           };
 
+  const Int_t    nPtbins     = 32;
+  Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
+                                  2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
+                                  6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
+                                 30.0, 36.0, 44.0 };
+
   TH2D   *histo = 0;
   TString hName;
 
-  hName  = Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data());
+  hName = Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data());
   histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaPi0], xMin[kEtaPi0], xMax[kEtaPi0], bins[kPhiPi0], xMin[kPhiPi0], xMax[kPhiPi0]);
-  histo->Sumw2(); list->Add(histo); histo = 0;
+  histo->Sumw2(); listRD->Add(histo); histo = 0;
   for (Int_t i=0; i<kVarsPi0-2; i++) {
-    hName  = Form("hPi0_%s%s", tName[0].Data(), tName[i+1].Data());
+    hName = Form("hPi0_%s%s", tName[0].Data(), tName[i+1].Data());
     histo = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    histo->Sumw2(); listRD->Add(histo); histo = 0;
   }
 
   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
-    hName  = Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data());
-    histo = new TH2D(hName.Data(), hName.Data(), bins[kPtPi0], xMin[kPtPi0], xMax[kPtPi0], bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    hName = Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data());
+    histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
+    histo->Sumw2(); listRD->Add(histo); histo = 0;
 
-    hName  = Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(), namePID[iPID].Data());
+    hName = Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(), namePID[iPID].Data());
     histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaPi0], xMin[kEtaPi0], xMax[kEtaPi0], bins[kPhiPi0], xMin[kPhiPi0], xMax[kPhiPi0]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    histo->Sumw2(); listRD->Add(histo); histo = 0;
     for (Int_t icent=0; icent<fgNCent; icent++) {
-      hName  = Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), icent);
-      histo = new TH2D(hName.Data(), hName.Data(), bins[kPtPi0], xMin[kPtPi0], xMax[kPtPi0], bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
-      histo->Sumw2(); list->Add(histo); histo = 0;
+      hName = Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), icent);
+      histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
+      histo->Sumw2(); listRD->Add(histo); histo = 0;
     }
 
     for (Int_t i=0; i<kVarsPi0-2; i++) {
-      hName  = Form("hPi0_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
+      hName = Form("hPi0_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
       histo = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
-      histo->Sumw2(); list->Add(histo); histo = 0;
+      histo->Sumw2(); listRD->Add(histo); histo = 0;
     }
     for (Int_t iComb=0; iComb<nComb; iComb++) {
-      hName  = Form("hPi0_%s%s_%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
-      histo = new TH2D(hName.Data(), hName.Data(), bins[kPtPi0], xMin[kPtPi0], xMax[kPtPi0], bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
-      histo->Sumw2(); list->Add(histo); histo = 0;
+      hName = Form("hPi0_%s%s_%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
+      histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
+      histo->Sumw2(); listRD->Add(histo); histo = 0;
     }
   }
 
-  TH2::AddDirectory(oldStatus);
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosMixPi0(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosMixPi0(TList *listRD)
 {
   // create Histograms for Mixed Pi0
 
-  if (!list) list = new TList();
-  list->SetOwner();
-
-  Bool_t oldStatus = TH2::AddDirectoryStatus(); TH2::AddDirectory(kFALSE);
   const Int_t nComb = 2;
-  TString  srcMod[nComb] =  { "Mod11", "Mod33" };
-                         // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs       };  // PID
-  TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"                 };
+  TString  srcMod[nComb]     = { "Mod11", "Mod33" };
+                             //{ kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs    };  // PID
+  TString namePID[kPIDs]     = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"              };
                              //{ kPtMixPi0, kEtaMixPi0, kPhiMixPi0, kInvMassMixPi0, kVarsMixPi0      };  // Mix Pi0
   TString tName[kVarsMixPi0] = {      "Pt",      "Eta",      "Phi",      "InvMass"                   };
   Int_t    bins[kVarsMixPi0] = {      500 ,     300   ,      120  ,          1000                    };
   Double_t xMin[kVarsMixPi0] = {        0.,      -0.15,        4.4,             0.                   };
   Double_t xMax[kVarsMixPi0] = {       50.,       0.15,        5.6,             1.                   };
 
+  const Int_t    nPtbins     = 32;
+  Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
+                                  2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
+                                  6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
+                                 30.0, 36.0, 44.0 };
+
   TH2D   *histo = 0;
   TString hName;
 
-  hName  = Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data());
+  hName = Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data());
   histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaMixPi0], xMin[kEtaMixPi0], xMax[kEtaMixPi0],
                                                bins[kPhiMixPi0], xMin[kPhiMixPi0], xMax[kPhiMixPi0]);
-  histo->Sumw2(); list->Add(histo); histo = 0;
+  histo->Sumw2(); listRD->Add(histo); histo = 0;
   
   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
-    hName  = Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data());
-    histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMixPi0],      xMin[kPtMixPi0],      xMax[kPtMixPi0],
-                                                 bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
-    hName  = Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(), namePID[iPID].Data());
+    hName = Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data());
+    histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
+    histo->Sumw2(); listRD->Add(histo); histo = 0;
+    hName = Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(), namePID[iPID].Data());
     histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaMixPi0], xMin[kEtaMixPi0], xMax[kEtaMixPi0], 
                                                  bins[kPhiMixPi0], xMin[kPhiMixPi0], xMax[kPhiMixPi0]);
-    histo->Sumw2(); list->Add(histo); histo = 0;
+    histo->Sumw2(); listRD->Add(histo); histo = 0;
     for (Int_t icent=0; icent<fgNCent; icent++) {
-      hName  = Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), icent);
-      histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMixPi0],      xMin[kPtMixPi0],      xMax[kPtMixPi0],
-                                                   bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
-      histo->Sumw2(); list->Add(histo); histo = 0;
+      hName = Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), icent);
+      histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
+      histo->Sumw2(); listRD->Add(histo); histo = 0;
     } 
     for (Int_t iComb=0; iComb<nComb; iComb++) {
-      hName  = Form("hMixPi0_%s%s_%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
-      histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMixPi0],      xMin[kPtMixPi0],      xMax[kPtMixPi0],
-                                                   bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
-      histo->Sumw2(); list->Add(histo); histo = 0;
+      hName = Form("hMixPi0_%s%s_%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
+      histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
+      histo->Sumw2(); listRD->Add(histo); histo = 0;
     }
   }
 
-  TH2::AddDirectory(oldStatus);
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::CreateHistosMC(TList *list)
+void AliPHOSpPbPi0Header::CreateHistosMC(TList *listMC)
 {
   // create Histograms for MC
-  if (!list) list = new TList();
-  list->SetOwner();
-
-  Bool_t oldStatus1 = TH2::AddDirectoryStatus(); TH1::AddDirectory(kFALSE);
-  Bool_t oldStatus2 = TH2::AddDirectoryStatus(); TH2::AddDirectory(kFALSE);
-  const Int_t method = 2;
-  const char *mName[method] = {"IsPrimary", "RCut"};
-  const Int_t cut    = 3;
-  const char *cName[cut]    = {"WideY", "NarrY", "Acc"};
-  const Int_t type   = 3;
-  const char *pName[type]   = {"Pi0", "Eta", "Gamma"};
-                       //  { kVertexMC, kPtMC, kRapidityMC,          kPhiMC, kWeightMC, kVarsMC };   // MC
-  TString tName[kVarsMC] = {  "Vertex",  "Pt",  "Rapidity",           "Phi",  "Weight"          };
-  Int_t    bins[kVarsMC] = {     1000 ,  500 ,        200 ,            360 ,     100            };
-  Double_t xMin[kVarsMC] = {        0.,    0.,         -1.,              0.,       0.           };
-  Double_t xMax[kVarsMC] = {      500.,   50.,          1.,  TMath::TwoPi(),      10.           };
+
+  const Int_t acc    = 2;
+  const char *aName[acc]   = {"Pi0InAcc", "GammaInAcc"};
+  const Int_t types  = 5;
+  const char *pName[types] = {"InclusivePi0", "PrimaryPi0", "2ndPi0FromK0s", "2ndPi0FromMaterials", "2ndPi0FromOtherDecays"};
+                         //  { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
+  TString namePID[kPIDs]   = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
+                         //  { kPtMC, kRadiusMC, kRapidityMC,          kPhiMC, kInvMassMC, kVarsMC  };  // MC
+  TString tName[kVarsMC]   = {  "Pt",  "Radius",  "Rapidity",           "Phi",  "InvMass"           };
+  Int_t    bins[kVarsMC]   = {  500 ,      1000,        200 ,            360 ,       500            };
+  Double_t xMin[kVarsMC]   = {    0.,        0.,         -1.,              0.,         0.           };
+  Double_t xMax[kVarsMC]   = {   50.,      500.,          1.,  TMath::TwoPi(),         0.5          };
+
+  const Int_t    nPtbins     = 32;
+  Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
+                                  2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
+                                  6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
+                                 30.0, 36.0, 44.0 };
 
   TH1D   *histo1 = 0;
   TH2D   *histo2 = 0;
   TString hName;
 
-  for (Int_t iType=0; iType<type; iType++) {
-    hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kVertexMC].Data());
-    histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kVertexMC], xMin[kVertexMC], xMax[kVertexMC]);
-    histo2->Sumw2(); list->Add(histo2); histo2 = 0;
+  for (Int_t iType=0; iType<types; iType++) {
+    hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data());
+    histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
+    histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
     hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kRapidityMC].Data());
     histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kRapidityMC], xMin[kRapidityMC], xMax[kRapidityMC]);
-    histo2->Sumw2(); list->Add(histo2); histo2 = 0;
-    hName  = Form("hMC%s_%s%s_%s", pName[iType], tName[kPtMC].Data(), tName[kVertexMC].Data(), mName[0]);
-    histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kVertexMC], xMin[kVertexMC], xMax[kVertexMC]);
-    histo2->Sumw2(); list->Add(histo2); histo2 = 0;
-    for (Int_t iMethod=0; iMethod<method; iMethod++) {
-      for (Int_t iCent=0; iCent<fgNCent; iCent++) {
-        for (Int_t i=0; i<kVarsMC-2; i++) {
-          hName  = Form("hMC%s_%s_cent%d_%s", pName[iType], tName[i+1].Data(), iCent, mName[iMethod]);
-          histo1 = new TH1D(hName.Data(), hName.Data(), bins[i+1], xMin[i+1], xMax[i+1] );
-          histo1->Sumw2(); list->Add(histo1); histo1 = 0;
+    histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
+    for (Int_t iAcc=0; iAcc<acc; iAcc++) {
+      hName  = Form("hMC%s_%s_%s", pName[iType], tName[kPtMC].Data(), aName[iAcc]);
+      histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
+      histo1->Sumw2(); listMC->Add(histo1); histo1 = 0;
+    }
+    for (Int_t iPID=0; iPID<kPIDs; iPID++) {
+      hName  = Form("hMC%s_%s%s_%s_Reco", pName[iType], tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data());
+      histo2 = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
+      histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
+    }
 
-        }
-        for (Int_t iCut=0; iCut<cut; iCut++) {
-          hName  = Form("hMC%s_%s_cent%d_%s_%s", pName[iType], tName[kPtMC].Data(), iCent, cName[iCut], mName[iMethod]);
-          histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
-          histo1->Sumw2(); list->Add(histo1); histo1 = 0;
-        }
+    for (Int_t iCent=0; iCent<fgNCent; iCent++) {
+      hName  = Form("hMC%s_%s%s_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), iCent);
+      histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
+      histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
+      hName  = Form("hMC%s_%s%s_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRapidityMC].Data(), iCent);
+      histo2 = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kRapidityMC], xMin[kRapidityMC], xMax[kRapidityMC]);
+      histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
+      for (Int_t iAcc=0; iAcc<acc; iAcc++) {
+        hName  = Form("hMC%s_%s_%s_cent%d", pName[iType], tName[kPtMC].Data(), aName[iAcc], iCent);
+        histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
+        histo1->Sumw2(); listMC->Add(histo1); histo1 = 0;
+      }
+      for (Int_t iPID=0; iPID<kPIDs; iPID++) {
+        hName  = Form("hMC%s_%s%s_%s_Reco_cent%d", pName[iType], tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data(), iCent);
+        histo2 = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
+        histo2->Sumw2(); listMC->Add(histo2); histo2 = 0;
       }
     }
   }
 
-
-  TH1::AddDirectory(oldStatus1);
-  TH2::AddDirectory(oldStatus2);
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosEvent(TList *list)
+void AliPHOSpPbPi0Header::FillHistosEvent(TList *listQA)
 {
   // fill histograms at event level according to event selection cuts
 
-  if (!list) return;
+  if (!listQA) return;
 
-  const Int_t nhs      = 3;
-  TString tName[nhs+1] = { "VtxNcontr", "Centrality",    "Vz", "Pileup" };
-  Double_t dist[nhs]   = { fVtxContrsN,  fCentrality, fVtx[2]           };
+  const Int_t nhs      = 2;
+  TString tName[nhs+1] = { "Centrality",    "Vz", "Pileup" };
+  Double_t dist[nhs]   = {  fCentrality, fVtx[2]           };
   for (Int_t i=nhs; i--;) {
-    ((TH1D*)list->FindObject(Form("hEvent_%s",tName[i].Data())))->Fill(dist[i]);
-    if (this->IsSelected()) ((TH1D*)list->FindObject(Form("hEventSel_%s",tName[i].Data())))->Fill(dist[i]);
+    ((TH1D*)listQA->FindObject(Form("hEvent_%s",tName[i].Data())))->Fill(dist[i]);
+    if (this->IsSelected()) ((TH1D*)listQA->FindObject(Form("hEventSel_%s",tName[i].Data())))->Fill(dist[i]);
   }
-  if (fIsPileupSPD) {
-    ((TH1D*)list->FindObject(Form("hEvent_%s",tName[nhs].Data())))->Fill(dist[nhs-1]);
-    if (this->IsSelected())  ((TH1D*)list->FindObject(Form("hEventSel_%s",tName[nhs].Data())))->Fill(dist[nhs-1]);
+  if (fIsPileup) {
+    ((TH1D*)listQA->FindObject(Form("hEvent_%s",tName[nhs].Data())))->Fill(dist[nhs-1]);
+    if (this->IsSelected())  ((TH1D*)listQA->FindObject(Form("hEventSel_%s",tName[nhs].Data())))->Fill(dist[nhs-1]);
   }
 
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosCaloCellsQA(TList *list, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo)
+void AliPHOSpPbPi0Header::FillHistosCaloCellsQA(TList *listQA, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo)
 {
   // fill QA histograms for calocells according to evnet selection cuts
 
-  if (!list) return;
+  if (!listQA) return;
 
   const Int_t nhs    = 3;
   TString tName[nhs] = { "CentNCells", "Nxz", "Exz" };
 
-  ((TH2D*)list->FindObject(Form("hCaloCellsQA_%s", tName[0].Data())))->Fill(fCentrality, cells->GetNumberOfCells());
+  ((TH2D*)listQA->FindObject(Form("hCaloCellsQA_%s", tName[0].Data())))->Fill(fCentrality, cells->GetNumberOfCells());
   for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
     Int_t relID[4] = {0,0,0,0};
     phosGeo->AbsToRelNumbering(cells->GetCellNumber(iCell),relID); // Int_t mod   = relID[0]; Int_t cellX = relID[2]; Int_t cellZ = relID[3];
     Double_t wight[2] = { 1., cells->GetAmplitude(iCell)};
-    ((TH1D*)list->FindObject(Form("hCaloCellsQA_%s", "E")))->Fill(wight[1]);
-    ((TH1D*)list->FindObject(Form("hCaloCellsQA_%s_Mod%d", "E", relID[0])))->Fill(wight[1]);
+    ((TH1D*)listQA->FindObject(Form("hCaloCellsQA_%s", "E")))->Fill(wight[1]);
+    ((TH1D*)listQA->FindObject(Form("hCaloCellsQA_%s_Mod%d", "E", relID[0])))->Fill(wight[1]);
     for (Int_t i=1; i<nhs; i++)
-      ((TH2D*)list->FindObject(Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), relID[0])))->Fill(relID[2], relID[3], wight[i-1]);
+      ((TH2D*)listQA->FindObject(Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), relID[0])))->Fill(relID[2], relID[3], wight[i-1]);
   }
 
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosCaloCluster(TList *list, TClonesArray* const caloClArr, Int_t cent)
+void AliPHOSpPbPi0Header::FillHistosCaloCluster(TList *listQA, TClonesArray* const caloClArr, Int_t cent)
 {
   // fill histograms for calocluster according to evnet && calocluster candidates selection cuts
 
-  if (!list)  return;
+  if (!listQA)  return;
                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs       };  // PID
-  TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"                 };
+  TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"                 };
   UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)     };
                          // { kPtClu, kEtaClu, kPhiClu, kM02Clu, kM20Clu, kTOFClu, kNCellsClu, kNClustersClu, kVarsClu };   // clusters
   TString tName[kVarsClu] = {   "Pt",   "Eta",   "Phi",   "M02",   "M20",   "TOF",   "NCells",   "NClusters"           };
@@ -616,45 +573,45 @@ void AliPHOSpPbPi0Header::FillHistosCaloCluster(TList *list, TClonesArray* const
                                   caloCluster->GetTOF(),
                                   (Double_t)caloCluster->GetNCells() };
  
-    ((TH1I*)list->FindObject(Form("hCaloCluster_%s", tName[kNCellsClu].Data())))->Fill((Int_t)vars[kNCellsClu]);
-    ((TH1I*)list->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNCellsClu].Data(), iMod)))->Fill((Int_t)vars[kNCellsClu]);
-    ((TH1D*)list->FindObject(Form("hCaloCluster_%s_Mod%d_TRU%d", tName[kPtClu].Data(), iMod, iTRU)))->Fill(vars[kPtClu]);
+    ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s", tName[kNCellsClu].Data())))->Fill((Int_t)vars[kNCellsClu]);
+    ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNCellsClu].Data(), iMod)))->Fill((Int_t)vars[kNCellsClu]);
+    ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d_TRU%d", tName[kPtClu].Data(), iMod, iTRU)))->Fill(vars[kPtClu]);
 
     for (Int_t iPID=0; iPID<kPIDs; iPID++) {
       if ((pidBit&PIDBIT[iPID])==PIDBIT[iPID]) {
-        ((TH1D*)list->FindObject(Form("hCaloCluster_%s_%s", tName[kPtClu].Data(), namePID[iPID].Data())))->Fill(vars[kPtClu]);
-        ((TH1D*)list->FindObject(Form("hCaloCluster_%s_Mod%d_%s", tName[kPtClu].Data(), iMod, namePID[iPID].Data())))->Fill(vars[kPtClu]);
-        ((TH1D*)list->FindObject(Form("hCaloCluster_%s_%s_cent%d", tName[kPtClu].Data(), namePID[iPID].Data(), cent)))->Fill(vars[kPtClu]);
-        ((TH2D*)list->FindObject(Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(),
+        ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_%s", tName[kPtClu].Data(), namePID[iPID].Data())))->Fill(vars[kPtClu]);
+        ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d_%s", tName[kPtClu].Data(), iMod, namePID[iPID].Data())))->Fill(vars[kPtClu]);
+        ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_%s_cent%d", tName[kPtClu].Data(), namePID[iPID].Data(), cent)))->Fill(vars[kPtClu]);
+        ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(),
                                                               namePID[iPID].Data())))->Fill(vars[kEtaClu], vars[kPhiClu]);
         for (Int_t j=0; j<kVarsClu-3; j++)
-          ((TH2D*)list->FindObject(Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[j+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[j+1]);
+          ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[j+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[j+1]);
       }
     }
 
-    ((TH2D*)list->FindObject(Form("hCaloCluster_%s%s", tName[kEtaClu].Data(), tName[kPhiClu].Data())))->Fill(vars[kEtaClu], vars[kPhiClu]);
+    ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s", tName[kEtaClu].Data(), tName[kPhiClu].Data())))->Fill(vars[kEtaClu], vars[kPhiClu]);
     for (Int_t j=0; j<kVarsClu-3; j++) {
-      ((TH2D*)list->FindObject(Form("hCaloCluster_%s%s", tName[0].Data(), tName[j+1].Data())))->Fill(vars[0], vars[j+1]);
+      ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s", tName[0].Data(), tName[j+1].Data())))->Fill(vars[0], vars[j+1]);
     }
   } // End loop calo cluster
 
-  ((TH1I*)list->FindObject(Form("hCaloCluster_%s", tName[kNClustersClu].Data())))->Fill(entries);
+  ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s", tName[kNClustersClu].Data())))->Fill(entries);
   for (Int_t jMod=1; jMod<4; jMod++) {
     if (jMod==2) continue;
-    ((TH1I*)list->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNClustersClu].Data(), jMod)))->Fill(nclusters[jMod-1]);
+    ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNClustersClu].Data(), jMod)))->Fill(nclusters[jMod-1]);
   }
 
   return;
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosPi0(TList *list, TClonesArray* const caloClArr, Int_t cent)
+void AliPHOSpPbPi0Header::FillHistosPi0(TList *listRD, TClonesArray* const caloClArr, Int_t cent)
 {
   // fill histograms for pi0 according to evnet && pi0 candidates selection cuts
 
-  if (!list) return;
+  if (!listRD) return;
                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
-  TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
+  TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
   UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
                          // { kPtPi0, kEtaPi0, kPhiPi0, kAsyPi0, kAnglePi0, kInvMassPi0, kVarsPi0  };  // Pi0
   TString tName[kVarsPi0] = {   "Pt",   "Eta",   "Phi",   "Asy",   "Angle",   "InvMass"            };
@@ -678,7 +635,7 @@ void AliPHOSpPbPi0Header::FillHistosPi0(TList *list, TClonesArray* const caloClA
       jGamma  = jCaloCluster->LorentzVector();
       if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
 
-      if (iMod==1 && jMod==1 )      srcMod = 11; // Mod11
+      if (iMod==1 && jMod==1)       srcMod = 11; // Mod11
       else if (iMod==3 && jMod==3)  srcMod = 33; // Mod33
       else if (iMod==2 && jMod==2)  srcMod = 22; // Mod22
       else if (iMod==1 && jMod==2)  srcMod = 12; // Mod12
@@ -691,23 +648,23 @@ void AliPHOSpPbPi0Header::FillHistosPi0(TList *list, TClonesArray* const caloClA
       vars[kAsyPi0]     = TMath::Abs((iGamma-jGamma).E())/(iGamma+jGamma).E();
       vars[kAnglePi0]   = jGamma.Angle(iGamma.Vect());
 
-      ((TH2D*)list->FindObject(Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data())))->Fill(vars[kEtaPi0], vars[kPhiPi0]);
+      ((TH2D*)listRD->FindObject(Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data())))->Fill(vars[kEtaPi0], vars[kPhiPi0]);
       for (Int_t k=0; k<kVarsPi0-2; k++) {
-        ((TH2D*)list->FindObject(Form("hPi0_%s%s", tName[0].Data(), tName[k+1].Data())))->Fill(vars[0], vars[k+1]);
+        ((TH2D*)listRD->FindObject(Form("hPi0_%s%s", tName[0].Data(), tName[k+1].Data())))->Fill(vars[0], vars[k+1]);
       }
 
       for (Int_t iPID=0; iPID<kPIDs; iPID++) {
         if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
-          ((TH2D*)list->FindObject(Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(),
+          ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(),
                                                         namePID[iPID].Data())))->Fill(vars[kEtaPi0], vars[kPhiPi0]);
-          ((TH2D*)list->FindObject(Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
+          ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
                                                         namePID[iPID].Data())))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
-          ((TH2D*)list->FindObject(Form("hPi0_%s%s_%s_Mod%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
+          ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s_Mod%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
                                                               namePID[iPID].Data(), srcMod)))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
-          ((TH2D*)list->FindObject(Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), 
+          ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), 
                                                                namePID[iPID].Data(), cent)))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
           for (Int_t k=0; k<kVarsPi0-2; k++) {
-            ((TH2D*)list->FindObject(Form("hPi0_%s%s_%s", tName[0].Data(), tName[k+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[k+1]);
+            ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[0].Data(), tName[k+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[k+1]);
           }
         }
       }
@@ -720,11 +677,11 @@ void AliPHOSpPbPi0Header::FillHistosPi0(TList *list, TClonesArray* const caloClA
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *list, TClonesArray* const iCaloClArr, TList* const eventList, Int_t cent)
+void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *listRD, TClonesArray* const iCaloClArr, TList* const eventList, Int_t cent)
 {
   // fill histograms for mixed pi0 according to evnet && pi0 candidates selection cuts
 
-  if (!list) return;
+  if (!listRD) return;
                         // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
   TString namePID[kPIDs] = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
   UInt_t  PIDBIT[kPIDs]  = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
@@ -765,17 +722,17 @@ void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *list, TClonesArray* const iCal
         vars[kPhiMixPi0]     = (iGamma+jGamma).Phi() + TMath::TwoPi();
         vars[kInvMassMixPi0] = (iGamma+jGamma).M();
 
-        ((TH2D*)list->FindObject(Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
+        ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
         for (Int_t iPID=0; iPID<kPIDs; iPID++) {
           if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
-            ((TH2D*)list->FindObject(Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(),
-                                                             namePID[iPID].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
-            ((TH2D*)list->FindObject(Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
-                                                             namePID[iPID].Data())))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
-            ((TH2D*)list->FindObject(Form("hMixPi0_%s%s_%s_Mod%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
-                                                                   namePID[iPID].Data(), srcMod)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
-            ((TH2D*)list->FindObject(Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
-                                                                    namePID[iPID].Data(), cent)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
+            ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(),
+                                                               namePID[iPID].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
+            ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
+                                                               namePID[iPID].Data())))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
+            ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s_Mod%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
+                                                                     namePID[iPID].Data(), srcMod)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
+            ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
+                                                                      namePID[iPID].Data(), cent)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
           }
         }
         jCaloCluster = 0;
@@ -789,213 +746,127 @@ void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *list, TClonesArray* const iCal
 }
 
 //_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosMC(TList *list, AliMCEvent* const mcEvent, AliPHOSGeoUtils* const phosGeo, Int_t cent)
+void AliPHOSpPbPi0Header::FillHistosMC(TList *listMC, AliStack* const stack, TClonesArray* const caloClArr, AliPHOSGeoUtils* const phosGeo, Int_t cent)
 {
-  // fill histograms for AOD MC
-  if (!list) return;
+  // fill histograms for ESD MC
+  if (!listMC)  return;
+  if (!stack)   return;
+
+                         // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
+  TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
+  UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
 
+  // MC truth info
   TString pName;
-  const Int_t ntrks = mcEvent->GetNumberOfTracks();
-  for(Int_t i=0; i<ntrks; i++) {
-    AliAODMCParticle *pMC = (AliAODMCParticle*)mcEvent->GetTrack(i);
-    Int_t pdg = pMC->GetPdgCode();
-    if (pdg == 111)      pName = "Pi0";   // Pi0
-    else if (pdg == 221) pName = "Eta";   // Eta
-    else if (pdg == 22)  pName = "Gamma"; // Gamma
-    else continue;
-
-    //Primary particle
+  const Int_t nParticles = stack->GetNtrack();
+  Int_t iMod = 0, jMod = 0;
+  for(Int_t iParticle=0; iParticle<nParticles; iParticle++) {
+    TParticle* pMC  = stack->Particle(iParticle);
+    if (pMC->GetPdgCode() != 111) continue;   // Pi0
+    pName = ClassifyMCPi0(iParticle, stack);
+
+    //Pi0 Information
     Double_t pt     = pMC->Pt();
-    Double_t r      = TMath::Sqrt(TMath::Power(pMC->Xv(), 2) + TMath::Power(pMC->Yv(), 2));
+    Double_t r      = pMC->R();
     Double_t y      = pMC->Y();
     Double_t phi    = pMC->Phi();   while (phi<0.) phi += TMath::TwoPi(); while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
-//  Double_t weight = PrimaryParticleWeight(pdg, pt);
-
-    ((TH2D*)list->FindObject(Form("hMC%s_PtVertex", pName.Data())))->Fill(pt, r);
-    ((TH2D*)list->FindObject(Form("hMC%s_PtRapidity", pName.Data())))->Fill(pt, y);
-
-    Int_t mod1 = 0, mod2 = 0;
-    AliAODMCParticle *gamma1=0x0, *gamma2=0x0;
-    if (r<1.) {
-      ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_RCut", pName.Data(), cent)))->Fill(pt);
-      ((TH1D*)list->FindObject(Form("hMC%s_Rapidity_cent%d_RCut", pName.Data(), cent)))->Fill(y);
-      ((TH1D*)list->FindObject(Form("hMC%s_Phi_cent%d_RCut", pName.Data(), cent)))->Fill(phi);
-      if (y<1.)   ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_WideY_RCut", pName.Data(), cent)))->Fill(pt); 
-      if (y<0.15) ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_NarrY_RCut", pName.Data(), cent)))->Fill(pt); 
-      if (pMC->GetNDaughters()==2) { // Do not account Dalitz decays
-        gamma1 = (AliAODMCParticle*)mcEvent->GetTrack(pMC->GetDaughter(0));
-        gamma2 = (AliAODMCParticle*)mcEvent->GetTrack(pMC->GetDaughter(1));
-        //Number of pi0s decayed into acceptance
-        mod1 = HitPHOSModule(gamma1, phosGeo);
-        mod2 = HitPHOSModule(gamma2, phosGeo);
-        if (mod1!=-1 && mod1!=2 && mod2!=-1 && mod2!=2 && TMath::Abs(mod1-mod2)<2) // !remove module 2
-          ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_Acc_RCut", pName.Data(), cent)))->Fill(pt);
-      }
+
+    ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius"))->Fill(pt, r);
+    ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_cent%d", cent)))->Fill(pt, r);
+    ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRapidity"))->Fill(pt, y);
+    ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRapidity_cent%d", cent)))->Fill(pt, y);
+
+    ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius", pName.Data())))->Fill(pt, r);
+    ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_cent%d", pName.Data(), cent)))->Fill(pt, r);
+    ((TH2D*)listMC->FindObject(Form("hMC%s_PtRapidity", pName.Data())))->Fill(pt, y);
+    ((TH2D*)listMC->FindObject(Form("hMC%s_PtRapidity_cent%d", pName.Data(), cent)))->Fill(pt, y);
+
+    if ((((phi>260.*TMath::DegToRad()) && (phi<280.*TMath::DegToRad())) || ((phi>300.*TMath::DegToRad()) && (phi<320.*TMath::DegToRad()))) && (TMath::Abs(y)<0.135)) {
+      ((TH1D*)listMC->FindObject("hMCInclusivePi0_Pt_Pi0InAcc"))->Fill(pt);
+      ((TH1D*)listMC->FindObject(Form("hMCInclusivePi0_Pt_Pi0InAcc_cent%d", cent)))->Fill(pt);
+
+      ((TH1D*)listMC->FindObject(Form("hMC%s_Pt_Pi0InAcc", pName.Data())))->Fill(pt);
+      ((TH1D*)listMC->FindObject(Form("hMC%s_Pt_Pi0InAcc_cent%d", pName.Data(), cent)))->Fill(pt);
     }
-    if (pMC->IsPhysicalPrimary()) {
-      ((TH2D*)list->FindObject(Form("hMC%s_PtVertex_IsPrimary", pName.Data())))->Fill(pt, r);
-      ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_IsPrimary", pName.Data(), cent)))->Fill(pt);
-      ((TH1D*)list->FindObject(Form("hMC%s_Rapidity_cent%d_IsPrimary", pName.Data(), cent)))->Fill(y);
-      ((TH1D*)list->FindObject(Form("hMC%s_Phi_cent%d_IsPrimary", pName.Data(), cent)))->Fill(phi);
-      if (y<1.)   ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_WideY_IsPrimary", pName.Data(), cent)))->Fill(pt); 
-      if (y<0.15) ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_NarrY_IsPrimary", pName.Data(), cent)))->Fill(pt); 
-      if (pMC->GetNDaughters()==2) { // Do not account Dalitz decays
-        gamma1 = (AliAODMCParticle*)mcEvent->GetTrack(pMC->GetDaughter(0));
-        gamma2 = (AliAODMCParticle*)mcEvent->GetTrack(pMC->GetDaughter(1));
-        //Number of pi0s decayed into acceptance
-        mod1 = HitPHOSModule(gamma1, phosGeo);
-        mod2 = HitPHOSModule(gamma2, phosGeo);
-        if (mod1!=-1 && mod1!=2 && mod2!=-1 && mod2!=2 && TMath::Abs(mod1-mod2)<2) // !remove module 2
-          ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_Acc_IsPrimary", pName.Data(), cent)))->Fill(pt);
+
+    TParticle *iGamma = 0x0, *jGamma = 0x0;
+    if (pMC->GetNDaughters()==2) { // Do not account Dalitz decays
+      iGamma = stack->Particle(pMC->GetFirstDaughter());
+      jGamma = stack->Particle(pMC->GetLastDaughter());
+      //Pi0's daughter particles in PHOS acceptance
+      iMod = HitPHOSModule(iGamma, phosGeo);
+      jMod = HitPHOSModule(jGamma, phosGeo);
+      if (iMod!=-1 && iMod!=2 && jMod!=-1 && jMod!=2 && TMath::Abs(iMod-jMod)<2) { // !remove module 2
+        ((TH1D*)listMC->FindObject("hMCInclusivePi0_Pt_GammaInAcc"))->Fill(pt);
+        ((TH1D*)listMC->FindObject(Form("hMCInclusivePi0_Pt_GammaInAcc_cent%d", cent)))->Fill(pt);
+
+        ((TH1D*)listMC->FindObject(Form("hMC%s_Pt_GammaInAcc", pName.Data())))->Fill(pt);
+        ((TH1D*)listMC->FindObject(Form("hMC%s_Pt_GammaInAcc_cent%d", pName.Data(), cent)))->Fill(pt);
       }
     }
   }
 
-  return;
-}
+  // Reconstruction info
+  Int_t  entries  = caloClArr->GetEntriesFast();
+  Int_t  iPi0Indx =-1, jPi0Indx =-1;
+  UInt_t iPIDBit  = 0, jPIDBit  = 0;
+  AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
+  TLorentzVector iLorentz, jLorentz;
+  for(Int_t i=0; i<entries-1; i++) { // Loop calo cluster i
+    iCaloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
+    if (!iCaloCluster->CheckIsClusterFromPi0(stack, iPi0Indx)) { iCaloCluster = 0; continue; }
+    iMod     = iCaloCluster->GetModule();
+    iPIDBit  = iCaloCluster->GetPIDBit();
+    iLorentz = iCaloCluster->LorentzVector();
 
-//_____________________________________________________________________________
-void AliPHOSpPbPi0Header::FillHistosMC(TList *list, AliStack* const stack, AliPHOSGeoUtils* const phosGeo, Int_t cent)
-{
-  // fill histograms for ESD MC
-  if (!list)  return;
-  if (!stack) return;
+    for (Int_t j=i+1; j<entries; j++) { // Loop calo cluster j
+      jCaloCluster = (AliCaloClusterInfo*)caloClArr->At(j);
+      if (!jCaloCluster->CheckIsClusterFromPi0(stack, jPi0Indx)) { jCaloCluster = 0; continue; }
+      if (jPi0Indx != iPi0Indx) { jCaloCluster = 0; continue; }  // coming from the same pi0
+      pName    = ClassifyMCPi0(iPi0Indx, stack);
+      jMod     = jCaloCluster->GetModule();
+      jPIDBit  = jCaloCluster->GetPIDBit();
+      jLorentz = jCaloCluster->LorentzVector();
+      if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
 
-  TString pName;
-  const Int_t ntrks = stack->GetNtrack();
-  for(Int_t i=0; i<ntrks; i++) {
-    TParticle* pMC = stack->Particle(i);
-    Int_t pdg = pMC->GetPdgCode();
-    if (pdg == 111)      pName = "Pi0";   // Pi0
-    else if (pdg == 221) pName = "Eta";   // Eta
-    else if (pdg == 22)  pName = "Gamma"; // Gamma
-    else continue;
-
-    //Primary particle
-    Double_t pt     = pMC->Pt();
-    Double_t r      = pMC->R();
-    Double_t y      = pMC->Y();
-    Double_t phi    = pMC->Phi();   while (phi<0.) phi += TMath::TwoPi(); while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
-//  Double_t weight = PrimaryParticleWeight(pdg, pt);
-
-    ((TH2D*)list->FindObject(Form("hMC%s_PtVertex", pName.Data())))->Fill(pt, r);
-    ((TH2D*)list->FindObject(Form("hMC%s_PtRapidity", pName.Data())))->Fill(pt, y);
-
-    Int_t mod1 = 0, mod2 = 0;
-    TParticle *gamma1 = 0x0, *gamma2 = 0x0;
-    if (r<1.) { // R_Cut
-      ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_RCut", pName.Data(), cent)))->Fill(pt);
-      ((TH1D*)list->FindObject(Form("hMC%s_Rapidity_cent%d_RCut", pName.Data(), cent)))->Fill(y);
-      ((TH1D*)list->FindObject(Form("hMC%s_Phi_cent%d_RCut", pName.Data(), cent)))->Fill(phi);
-      if (y<0.5)  ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_WideY_RCut", pName.Data(), cent)))->Fill(pt);
-      if (y<0.15) ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_NarrY_RCut", pName.Data(), cent)))->Fill(pt);
-      if (pMC->GetNDaughters()==2) { // Do not account Dalitz decays
-        gamma1 = stack->Particle(pMC->GetFirstDaughter());
-        gamma2 = stack->Particle(pMC->GetLastDaughter());
-        //Number of pi0s decayed into acceptance
-        mod1 = HitPHOSModule(gamma1, phosGeo);
-        mod2 = HitPHOSModule(gamma2, phosGeo);
-        if (mod1!=-1 && mod1!=2 && mod2!=-1 && mod2!=2 && TMath::Abs(mod1-mod2)<2) // !remove module 2
-          ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_Acc_RCut", pName.Data(), cent)))->Fill(pt);
-      }
-    }
-    if (stack->IsPhysicalPrimary(i)) { // IsPrimary
-      ((TH2D*)list->FindObject(Form("hMC%s_PtVertex_IsPrimary", pName.Data())))->Fill(pt, r);
-      ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_IsPrimary", pName.Data(), cent)))->Fill(pt);
-      ((TH1D*)list->FindObject(Form("hMC%s_Rapidity_cent%d_IsPrimary", pName.Data(), cent)))->Fill(y);
-      ((TH1D*)list->FindObject(Form("hMC%s_Phi_cent%d_IsPrimary", pName.Data(), cent)))->Fill(phi);
-      if (y<0.5)  ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_WideY_IsPrimary", pName.Data(), cent)))->Fill(pt);
-      if (y<0.15) ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_NarrY_IsPrimary", pName.Data(), cent)))->Fill(pt);
-      if (pMC->GetNDaughters()==2) { // Do not account Dalitz decays
-        gamma1 = stack->Particle(pMC->GetFirstDaughter());
-        gamma2 = stack->Particle(pMC->GetLastDaughter());
-        //Number of pi0s decayed into acceptance
-        mod1 = HitPHOSModule(gamma1, phosGeo);
-        mod2 = HitPHOSModule(gamma2, phosGeo);
-        if (mod1!=-1 && mod1!=2 && mod2!=-1 && mod2!=2 && TMath::Abs(mod1-mod2)<2) // !remove module 2
-          ((TH1D*)list->FindObject(Form("hMC%s_Pt_cent%d_Acc_IsPrimary", pName.Data(), cent)))->Fill(pt);
+      Double_t pi0Pt      = (iLorentz+jLorentz).E();
+      Double_t pi0InvMass = (iLorentz+jLorentz).M();
+
+      for (Int_t iPID=0; iPID<kPIDs; iPID++) {
+        if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
+          ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtInvMass_%s_Reco", namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
+          ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtInvMass_%s_Reco_cent%d", namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
+
+          ((TH2D*)listMC->FindObject(Form("hMC%s_PtInvMass_%s_Reco", pName.Data(), namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
+          ((TH2D*)listMC->FindObject(Form("hMC%s_PtInvMass_%s_Reco_cent%d", pName.Data(), namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
+        }
       }
-    }
-  }
+    } // End loop calo cluster j 
+  } // End loop calo cluster i
 
   return;
 }
 
 //________________________________________________________________________
-Double_t AliPHOSpPbPi0Header::PrimaryParticleWeight(Int_t pdg, Double_t pt)
+TString AliPHOSpPbPi0Header::ClassifyMCPi0(Int_t index, AliStack* const stack)
 {
 
-  Int_t type=0 ;
-  if (pdg == 111 || TMath::Abs(pdg)==211) type = 1; // Pi0/Eta 
-  else if (TMath::Abs(pdg)<1000)          type = 2; // Kaon-like
-  else                                    type = 3; // baryons
-
-  if (type==1) {
-    if (fCentrality<5.)       // 0-5
-      return (1.662990+1.140890*pt-0.192088*pt*pt)/(1.-0.806630*pt+0.304771*pt*pt)+0.141690*pt;
-    else if (fCentrality<10.) // 5-10
-      return (1.474351+0.791492*pt-0.066369*pt*pt)/(1.-0.839338*pt+0.317312*pt*pt)+0.093289*pt;
-    else if (fCentrality<20.) // 10-20
-      return (1.174728+0.959681*pt-0.137695*pt*pt)/(1.-0.788873*pt+0.299538*pt*pt)+0.128759*pt; 
-    else if (fCentrality<40.) // 20-40
-      return (0.927335+0.475349*pt+0.004364*pt*pt)/(1.-0.817966*pt+0.309787*pt*pt)+0.086899*pt; 
-    else if (fCentrality<60.) // 40-60
-      return (0.676878+0.190680*pt+0.077031*pt*pt)/(1.-0.790623*pt+0.305183*pt*pt)+0.064510*pt; 
-    else if (fCentrality<80.) // 60-80
-      return (0.684726-0.606262*pt+0.409963*pt*pt)/(1.-1.080061*pt+0.456933*pt*pt)+0.005151*pt; 
-  }
-  if (type==2) {
-    if (fCentrality<5.)       // 0-5
-      return (-0.417131+2.253936*pt-0.337731*pt*pt)/(1.-0.909892*pt+0.316820*pt*pt)+0.157312*pt;
-    else if (fCentrality<10.) // 5-10
-      return (-0.352275+1.844466*pt-0.248598*pt*pt)/(1.-0.897048*pt+0.316462*pt*pt)+0.132461*pt; 
-    else if (fCentrality<20.) // 10-20
-      return (-0.475481+1.975108*pt-0.336013*pt*pt)/(1.-0.801028*pt+0.276705*pt*pt)+0.188164*pt; 
-    else if (fCentrality<40.) // 20-40
-      return (-0.198954+1.068789*pt-0.103540*pt*pt)/(1.-0.848354*pt+0.299209*pt*pt)+0.112939*pt; 
-    else if (fCentrality<60.) // 40-60
-      return (-0.111052+0.664041*pt-0.019717*pt*pt)/(1.-0.804916*pt+0.300779*pt*pt)+0.095784*pt;
-    else if (fCentrality<80.) // 60-80
-      return (0.202788-0.439832*pt+0.564585*pt*pt)/(1.-1.254029*pt+0.679444*pt*pt)+0.016235*pt;
-  }
-  if (type==3) {
-    if (fCentrality<5.)       // 0-5
-      return (-1.312732+2.743568*pt-0.375775*pt*pt)/(1.-0.717533*pt+0.164694*pt*pt)+0.164445*pt;
-    else if (fCentrality<10.) // 5-10
-      return (-1.229425+2.585889*pt-0.330164*pt*pt)/(1.-0.715892*pt+0.167386*pt*pt)+0.133085*pt;
-    else if (fCentrality<20.) // 10-20
-      return (-1.135677+2.397489*pt-0.320355*pt*pt)/(1.-0.709312*pt+0.164350*pt*pt)+0.146095*pt;
-    else if (fCentrality<40.) // 20-40
-      return (-0.889993+1.928263*pt-0.220785*pt*pt)/(1.-0.715991*pt+0.174729*pt*pt)+0.095098*pt;
-    else if (fCentrality<60.) // 40-60
-      return (-0.539237+1.329118*pt-0.115439*pt*pt)/(1.-0.722906*pt+0.186832*pt*pt)+0.059267*pt;
-    else if (fCentrality<80.) // 60-80
-      return (-0.518126+1.327628*pt-0.130881*pt*pt)/(1.-0.665649*pt+0.184300*pt*pt)+0.081701*pt;
-  }
-
-  return 1.;
-}
+    TString  pName;
+    TParticle* pMC  = stack->Particle(index);
+    Int_t      mpdg = 0;
 
-//________________________________________________________________________
-Int_t AliPHOSpPbPi0Header::HitPHOSModule(AliAODMCParticle* const pMC, AliPHOSGeoUtils* const phosGeo)
-{ 
+    if (pMC->GetFirstMother()>-1) {
+      TParticle* pMother = stack->Particle(pMC->GetFirstMother());
+      mpdg = TMath::Abs(pMother->GetPdgCode());
+    }
 
-  Int_t mod=0, relID[4]={0,0,0,0}; 
-  Double_t x=0., z=0.;
-  Double_t vtx[3]; pMC->XvYvZv(vtx);
-  Double_t theta = pMC->Theta();
-  Double_t phi   = pMC->Phi();
-    
-  if (!phosGeo->ImpactOnEmc(vtx, theta, phi, mod, z, x)) return -1;
-  phosGeo->RelPosToRelId(mod, x, z, relID);
-  if (fgUseFiducialCut) {
-    const Int_t edgeX = 2;
-    const Int_t edgeZ = 2;
-    if (relID[2] >edgeX && relID[2]<(65-edgeX) && relID[3]>edgeZ && relID[3] <(57-edgeZ)) return relID[0];
-    else return -1;
-  }
+    // Classify Pi0 sources
+    if (index<stack->GetNprimary())                              pName = "PrimaryPi0";
+    else if (mpdg==310)                                          pName = "2ndPi0FromK0s";
+    else if (mpdg==321 || mpdg==211 || mpdg==2212 || mpdg==2112) pName = "2ndPi0FromMaterials";     // K+-, pi+-, p(pbar), n(nbar)
+    else                                                         pName = "2ndPi0FromOtherDecays";   // all other secondary pi0s
 
-  return relID[0];
+    return pName;
 }
 
 //________________________________________________________________________
index 0ef68f1..261767e 100644 (file)
@@ -19,11 +19,8 @@ class TParticle;
 class TClonesArray;
 
 class AliInputEventHandler;
-class AliMCEvent;
 class AliStack;
-class AliAODMCParticle;
-class AliAODEvent;
-class AliESDEvent;
+class AliVEvent;
 class AliVCaloCells;
 class AliPHOSGeoUtils;
 
@@ -43,67 +40,62 @@ class AliPHOSpPbPi0Header : public TNamed {
   Double_t Vz()                    const { return fVtx[2];                         }
   TString  FiredTriggerClass()     const { return fFiredTriggerClass;              }
   UInt_t   SelectionMask()         const { return fSelMask;                        }
-  Int_t    VtxContrsN()            const { return fVtxContrsN;                     }
   Bool_t   IsVertexOK()            const { return fIsVertexOK;                     }
-  Bool_t   IsPileupSPD()           const { return fIsPileupSPD;                    }
+  Bool_t   IsPileup()              const { return fIsPileup;                       }
   Float_t  Centrality()            const { return fCentrality;                     }
+
   Bool_t   IsSelected();
 
   void SetEventInfo(AliInputEventHandler* const handler);
 
   void CreateHistograms(TList *listQA, TList *listRD, TList *listMC);
-  void FillHistosEvent(TList *list);
-  void FillHistosCaloCellsQA(TList *list, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo);
-  void FillHistosCaloCluster(TList *list, TClonesArray* const caloClArr, Int_t cent);
-  void FillHistosPi0(TList *list, TClonesArray* const caloClArr, Int_t cent);
-  void FillHistosMixPi0(TList *list, TClonesArray* const caloClArr, TList* const eventlist, Int_t cent);
-  void FillHistosMC(TList *list, AliMCEvent* const mcEvent, AliPHOSGeoUtils* const phosGeo, Int_t cent);
-  void FillHistosMC(TList *list, AliStack*   const stack,   AliPHOSGeoUtils* const phosGeo, Int_t cent);
-
-  static void SetIsMC(Bool_t isMC=kFALSE)         { fgIsMC                         = isMC;    }
-  static void SetIspARun(Bool_t ispARun=kFALSE)   { fgIspARun                      = ispARun; }
-  static void SetUseFiducialCut(Bool_t fc=kFALSE) { fgUseFiducialCut               = fc;      }
-  static void SetNCent(Int_t ncent=10)            { fgNCent                        = ncent;   }
-  static void SetSelectionCuts(Double_t cuts[4])  { for (Int_t i=4; i--;) fgCuts[i]=cuts[i];  }
+  void FillHistosEvent(TList *listQA);
+  void FillHistosCaloCellsQA(TList *listQA, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo);
+  void FillHistosCaloCluster(TList *listQA, TClonesArray* const caloClArr, Int_t cent);
+  void FillHistosPi0(TList *listRD, TClonesArray* const caloClArr, Int_t cent);
+  void FillHistosMixPi0(TList *listRD, TClonesArray* const caloClArr, TList* const eventlist, Int_t cent);
+  void FillHistosMC(TList *listMC, AliStack* const stack, TClonesArray* const caloClArr, AliPHOSGeoUtils* const phosGeo, Int_t cent);
+
+  static void SetIsMC(Bool_t isMC=kFALSE)         { fgIsMC                          = isMC;    }
+  static void SetUseFiducialCut(Bool_t fc=kFALSE) { fgUseFiducialCut                = fc;      }
+  static void SetNCent(Int_t ncent=10)            { fgNCent                         = ncent;   }
+  static void SetSelectionCuts(Double_t cuts[3])  { for (Int_t i=3; i--;) fgCuts[i] = cuts[i]; }
 
  private :
 
-  void CreateHistosEvent(TList *list);
-  void CreateHistosCaloCellsQA(TList *list);
-  void CreateHistosCaloCluster(TList *list);
-  void CreateHistosPi0(TList *list);
-  void CreateHistosMixPi0(TList *list);
-  void CreateHistosMC(TList *list);
+  void CreateHistosEvent(TList *listQA);
+  void CreateHistosCaloCellsQA(TList *listQA);
+  void CreateHistosCaloCluster(TList *listQA);
+  void CreateHistosPi0(TList *listRD);
+  void CreateHistosMixPi0(TList *listRD);
+  void CreateHistosMC(TList *listMC);
 
-  Bool_t    CheckEventVertex(AliAODEvent* const aod, AliESDEvent* const esd);
-  Int_t     HitPHOSModule(AliAODMCParticle* const pMC, AliPHOSGeoUtils* const phosGeo);
+  Bool_t    CheckEventVertex(AliVEvent* const event);
+  TString   ClassifyMCPi0(Int_t index, AliStack* const stack);
   Int_t     HitPHOSModule(TParticle* const pMC, AliPHOSGeoUtils* const phosGeo);
-  Double_t  PrimaryParticleWeight(Int_t pdg, Double_t pt);
 
   static Bool_t   fgIsMC;           // flag to use MC
   static Bool_t   fgIspARun;        // flag to use pA vertex cut
   static Bool_t   fgUseFiducialCut; // flag to use fiducial cut
   static Int_t    fgNCent;          // # of centrality bins
-  static Double_t fgCuts[4];        // 0, low limit of num. of vtx contributors
-                                    // 1, up limit of vz
-                                    // 2, centrality max
-                                    // 3, centrality min
+  static Double_t fgCuts[3];        // 0, up limit of vz
+                                    // 1, centrality max
+                                    // 2, centrality min
 
   enum { kAll,      kCpv,       kDisp,       kBoth,          kCpv2,      kDisp2,      kBoth2,     kPIDs                   };   // PID
   enum { kPtClu,    kEtaClu,    kPhiClu,     kM02Clu,        kM20Clu,    kTOFClu,     kNCellsClu, kNClustersClu, kVarsClu };   // clusters
   enum { kPtPi0,    kEtaPi0,    kPhiPi0,     kAsyPi0,        kAnglePi0,  kInvMassPi0, kVarsPi0                            };   // pi0
   enum { kPtMixPi0, kEtaMixPi0, kPhiMixPi0,  kInvMassMixPi0, kVarsMixPi0                                                  };   // Mixed pi0
-  enum { kVertexMC, kPtMC,      kRapidityMC, kPhiMC,         kWeightMC,  kVarsMC                                          };   // MC
+  enum { kPtMC,     kRadiusMC,  kRapidityMC, kPhiMC,         kInvMassMC, kVarsMC                                          };   // MC
 
-  Double_t fVtx[3];                       // position of vtx
-  TString  fFiredTriggerClass;            // trigger class
-  UInt_t   fSelMask;                      // mask of physics selection
-  Int_t    fVtxContrsN;                   // num. of contributors of vtx rec
-  Bool_t   fIsVertexOK;                   // is vertex OK
-  Bool_t   fIsPileupSPD;                  // is Pileup from SPD
-  Float_t  fCentrality;                   // event certrality
+  Double_t fVtx[3];                 // position of vtx
+  TString  fFiredTriggerClass;      // trigger class
+  UInt_t   fSelMask;                // mask of physics selection
+  Bool_t   fIsVertexOK;             // is vertex OK
+  Bool_t   fIsPileup;               // is Pileup from SPD
+  Float_t  fCentrality;             // event certrality
 
-  ClassDef(AliPHOSpPbPi0Header, 1)
+  ClassDef(AliPHOSpPbPi0Header, 2)
 };
 
 #endif
index 92269b3..300b521 100644 (file)
@@ -1,7 +1,7 @@
-AliAnalysisTaskSEPHOSpPbPi0* AddTaskPHOSpPbPi0(UInt_t triggerTag = 0, Bool_t isMCtruth=kFALSE, Bool_t isBadMap=kFALSE, Double_t width = 0.)
+AliAnalysisTaskSEPHOSpPbPi0* AddTaskPHOSpPbPi0(UInt_t triggerTag = 0, Bool_t isMCtruth=kFALSE)
 {
 // Creates a task to analysis pi0 in p-Pb collisions with PHOS
-// Author: H. Zhu - 05/14/2013
+// Author: H. Zhu - 05/28/2013
 
   AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
   if (!mgr) {
@@ -25,42 +25,29 @@ AliAnalysisTaskSEPHOSpPbPi0* AddTaskPHOSpPbPi0(UInt_t triggerTag = 0, Bool_t isM
   Float_t cBin[nBins+1] = {0., 20., 40., 60., 80., 90., 100.}; TArrayF tCent(nBins+1, cBin);
   Int_t   buffer[nBins] = { 40,  80,  80,  100, 100, 100    }; TArrayI tBuffer(nBins, buffer);
 
-  Double_t cutsEvent[4] = {  1.,    // 0, min of VtxNCont
-                            10.,    // 1, max of VtxZ
-                             0.,    // 2, min of Centrality
-                           100.  }; // 3, max of Centrality
+  Double_t cutsEvent[3] = { 10.,    // 0, max of VtxZ
+                             0.,    // 1, min of Centrality
+                           100.  }; // 2, max of Centrality
 
   Double_t cutsCaloCl[5] = { 0.3,    // 0, min of cluster Energy
                              2.,     // 1, min of NCells
                              0.2,    // 2, min of M02
                              2.5,    // 3, min of DistToBad
-                             7e-8 }; // 4, max of TOF
+                             1e-7 }; // 4, max of TOF
 
   AliAnalysisTaskSEPHOSpPbPi0 *task = new AliAnalysisTaskSEPHOSpPbPi0("TaskPHOSpPbPi0");
   task->SetUseMC(isMCtruth);
   task->SetXBins(tCent, tBuffer);
   task->SetEventCuts(cutsEvent);
   task->SetCaloClCuts(cutsCaloCl);
-  task->SetDecaliWidth(width);    // for decalibration
-  task->SetRemovePileup(kFALSE);  // not remove pileup events
-  task->SetUseFiducialCut(kTRUE); // use fiducial cut
+  task->SetRemovePileup(kTRUE);    // remove pileup events
+  task->SetUseFiducialCut(kTRUE);  // use fiducial cut
+  task->SetUseTOFCut(kTRUE);       // use TOF cut
   task->SetDebugLevel(-1);
 
   if (triggerTag==0) task->SelectCollisionCandidates(AliVEvent::kINT7);
-  if (triggerTag==1) task->SelectCollisionCandidates(AliVEvent::kPHI7);
-  if (triggerTag==2) task->SelectCollisionCandidates(AliVEvent::kMB);
-
-  // Bad maps for PHOS
-  if (isBadMap) {
-    AliOADBContainer badmapContainer("phosBadMap");
-    badmapContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSBadMaps.root","phosBadMap");
-    TObjArray *maps = (TObjArray*)badmapContainer.GetObject(144871,"phosBadMap"); //run number LHC11a
-
-    for (Int_t mod=1;mod<4; mod++) {
-      TH2I * h = (TH2I*)maps->At(mod);
-      if(h) task->SetPHOSBadMap(mod,h);
-    }
-  }
+  else if (triggerTag==1) task->SelectCollisionCandidates(AliVEvent::kPHI7);
+  else if (triggerTag==2) task->SelectCollisionCandidates(AliVEvent::kMB);
 
   mgr->AddTask(task);