user task from hongshen
authorloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 27 Apr 2013 13:33:28 +0000 (13:33 +0000)
committerloizides <loizides@f7af4fe6-9843-0410-8265-dc069ae4e863>
Sat, 27 Apr 2013 13:33:28 +0000 (13:33 +0000)
PWGGA/CMakelibPWGGAPHOSTasks.pkg
PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.h [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.cxx [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.h [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.h [new file with mode: 0644]
PWGGA/PHOSTasks/UserTasks/macros/AddTaskPHOSpPbPi0.C [new file with mode: 0644]
PWGGA/PWGGAPHOSTasksLinkDef.h

index fcd0df3..c1e0381 100644 (file)
@@ -41,7 +41,10 @@ set ( SRCS
  PHOSTasks/CaloCellQA/AliCaloCellsQA.cxx
  PHOSTasks/CaloCellQA/AliAnalysisTaskCaloCellsQA.cxx
  PHOSTasks/omega3pi/AliAnalysisTaskOmegaPi0PiPi.cxx
- )
+ PHOSTasks/UserTasks/AliCaloClusterInfo.cxx
+ PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx
+ PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx
+)
 
 string ( REPLACE ".cxx" ".h" HDRS "${SRCS}" )
 
diff --git a/PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx b/PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.cxx
new file mode 100644 (file)
index 0000000..0204292
--- /dev/null
@@ -0,0 +1,289 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+///////////////////////////////////////////////////////////////////////////
+//
+// AliAnalysisTaskSE for the gamma and pi0 from pPb collision analysis
+//
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch
+//                   hszhu@iopp.ccnu.edu.cn
+///////////////////////////////////////////////////////////////////////////
+
+#include <TH2I.h>
+#include <TList.h>
+#include <TMath.h>
+#include <TArray.h>
+#include <TClonesArray.h>
+
+#include "AliAnalysisManager.h"
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliVEvent.h"
+#include "AliESDEvent.h"
+#include "AliAODEvent.h"
+#include "AliESDHeader.h"
+#include "AliAODHeader.h"
+#include "AliVCluster.h"
+#include "AliVCaloCells.h"
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSGeometry.h"
+#include "AliOADBContainer.h"
+#include "AliInputEventHandler.h"
+#include "AliAnalysisManager.h"
+#include "AliPHOSpPbPi0Header.h"
+#include "AliCaloClusterInfo.h"
+#include "AliAnalysisTaskSEPHOSpPbPi0.h"
+
+ClassImp(AliAnalysisTaskSEPHOSpPbPi0)
+
+Int_t    AliAnalysisTaskSEPHOSpPbPi0::fgMinNCells        = 2;
+Double_t AliAnalysisTaskSEPHOSpPbPi0::fgMinClusterEnergy = 0.3;
+Double_t AliAnalysisTaskSEPHOSpPbPi0::fgMinM02           = 0.2;
+Double_t AliAnalysisTaskSEPHOSpPbPi0::fgMinDistToBad     = 2.5;
+
+//________________________________________________________________________
+AliAnalysisTaskSEPHOSpPbPi0::AliAnalysisTaskSEPHOSpPbPi0():
+  AliAnalysisTaskSE(), fIsMC(kFALSE), fCentralityBin(10), fBufferSize(10),
+  fRunNumber(-1), fPHOSGeo(0), fList(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; }
+
+  // Set 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), fList(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; }
+
+  // Set 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());
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskSEPHOSpPbPi0::~AliAnalysisTaskSEPHOSpPbPi0()
+{
+  //
+  // Default destructor
+  //
+  if (fList)        { delete fList;        fList      = NULL; }
+  if (fHeader)      { delete fHeader;      fHeader    = NULL; }
+  if (fCaloClArr)   { delete fCaloClArr;   fCaloClArr = NULL; }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEPHOSpPbPi0::UserCreateOutputObjects()
+{
+  // Create the output container
+  // Initialize the PHOS geometry
+//fPHOSGeo = new AliPHOSGeoUtils("PHOSGeo");
+//fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP");
+
+  if (!fHeader)    fHeader    = new AliPHOSpPbPi0Header();
+  if (!fCaloClArr) fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);
+  if (!fList)      fList      = new TList();
+
+  fHeader->SetIsMC(fIsMC);
+  fHeader->SetNCent(fCentralityBin.GetSize());
+  fHeader->CreateHistograms(fList);
+
+  // Post output data.
+  PostData(1, fList);
+  return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEPHOSpPbPi0::UserExec(Option_t *) 
+{
+  // Main loop
+  // Called for each event
+
+  if (fIsMC) {
+    if (MCEvent()) {
+      if (MCEvent()->GetNumberOfTracks()<=0)
+           { AliError("MC event not found. Nothing done!"); return; }
+    } else { AliError("MC event not found. Nothing done!"); return; }
+  }
+
+  Int_t       nclsts = 0;
+  AliAODEvent *aod   = 0;
+  AliESDEvent *esd   = 0;
+
+  if (((TString)fInputEvent->IsA()->GetName())=="AliAODEvent") {
+    aod = dynamic_cast<AliAODEvent*>(fInputEvent);
+    if (!aod) { AliError("AOD event not found. Nothing done!");   return; }
+    if (!fIsMC && (aod->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
+    if (!fHeader->IspAVertexOK(aod))                              return; // check p-A collision vertex
+    nclsts = aod->GetNumberOfCaloClusters();   if (nclsts<1)      return;
+  } else {
+    esd = dynamic_cast<AliESDEvent*>(fInputEvent);
+    if (!esd) { AliError("ESD event not found. Nothing done!");   return; }
+    if (!fIsMC && (esd->GetHeader()->GetEventType()!=7))          return; // check event type; should be PHYSICS = 7 for data and 0 for MC
+    if (!fHeader->IspAVertexOK(esd))                              return; // check p-A collision vertex
+    nclsts = esd->GetNumberOfCaloClusters();   if (nclsts<1)      return;
+  }
+
+  // Fill Event info
+  fHeader->SetEventInfo(fInputHandler);
+  fHeader->FillHistosEvnH(fList);
+
+  // PHOS Geometry and Misalignment initialization at the first time it runs
+  if(fRunNumber != fInputEvent->GetRunNumber()) {
+    fRunNumber = fInputEvent->GetRunNumber();
+    PHOSInitialize(esd);
+  }
+
+  // Fill PHOS cells QA histograms
+  fHeader->FillHistosCaloCellsQA(fList, fInputEvent->GetPHOSCells(), fPHOSGeo);
+//aod = 0;
+
+  // Fill PHOS cluster Clones Array
+  FillCaloClusterInfo(nclsts, esd);
+
+  if (!fCaloClArr->GetEntriesFast()) return;
+  Int_t zvtx = (Int_t)((fHeader->Vz() + 10.)/2.);    if (zvtx<0) zvtx = 0; if (zvtx>9) zvtx = 9;
+  Int_t cent = TMath::BinarySearch<Double_t>(fCentralityBin.GetSize()-1, fCentralityBin.GetArray(), fHeader->Centrality());
+  if (!fEventList[zvtx][cent]) fEventList[zvtx][cent] = new TList(); 
+  TList *eventList = fEventList[zvtx][cent];
+
+  // Fill cluster histograms
+  fHeader->FillHistosCaloCluster(fList, fCaloClArr, cent);
+
+  // Fill pi0 histograms
+  fHeader->FillHistosPi0(fList, fCaloClArr, cent);
+
+  // Fill mixed pi0 histograms
+  fHeader->FillHistosMixPi0(fList, fCaloClArr, eventList, cent);
+
+  // Fill MC info
+  AliStack *stack = 0;
+  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(fList, stack, fPHOSGeo, cent);
+    } else
+      fHeader->FillHistosMC(fList, MCEvent(), fPHOSGeo, cent);
+  }
+//esd = 0;
+
+  // Fill event list for mixing
+  if(fCaloClArr->GetEntriesFast()>0) {
+    eventList->AddFirst(fCaloClArr);   fCaloClArr = 0;
+    //fCaloClArr->Clear();
+    if(eventList->GetSize()>fBufferSize[cent]) { // Remove redundant events
+      TClonesArray * tmp = static_cast<TClonesArray*>(eventList->Last()) ;
+      eventList->RemoveLast();
+      delete tmp ;
+    }
+  }
+
+  return;
+}      
+
+//________________________________________________________________________
+void AliAnalysisTaskSEPHOSpPbPi0::Terminate(Option_t *) 
+{
+  // Terminate analysis
+
+  // add the correction matrix
+
+  return;
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEPHOSpPbPi0::PHOSInitialize(AliESDEvent* const esd)
+{
+  // Initialize PHOS Geometry and misalignment
+  // PHOS Geometry
+  AliOADBContainer geomContainer("phosGeo");
+  geomContainer.InitFromFile("$ALICE_ROOT/OADB/PHOS/PHOSGeometry.root","PHOSRotationMatrixes");
+  TObjArray *matrixes = (TObjArray*)geomContainer.GetObject(fRunNumber,"PHOSRotationMatrixes");
+  fPHOSGeo =  AliPHOSGeometry::GetInstance("IHEP");
+  for (Int_t mod=0; mod<5; mod++) {
+    if (!matrixes->At(mod)) continue;
+    else fPHOSGeo->SetMisalMatrix(((TGeoHMatrix*)matrixes->At(mod)),mod);
+  }
+
+  // sets the PHOS Misalignment vertex if ESD
+  if (esd) {
+    for (Int_t mod=0; mod<5; mod++) {
+      const TGeoHMatrix* modMatrix = fInputEvent->GetPHOSMatrix(mod);
+      if (!modMatrix) continue;
+      else fPHOSGeo->SetMisalMatrix(modMatrix, mod);
+    }
+  }
+}
+
+//________________________________________________________________________
+void AliAnalysisTaskSEPHOSpPbPi0::FillCaloClusterInfo(Int_t nclsts, AliESDEvent* const esd)
+{
+  // Fill calo cluster info
+  if (fCaloClArr) fCaloClArr->Clear();
+  else fCaloClArr = new TClonesArray("AliCaloClusterInfo", 0);
+
+  TClonesArray       &caloRef     = *fCaloClArr;
+  Int_t              countN       = 0;
+  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);
+  TLorentzVector     momentum;
+  AliVCluster        *clust       = 0;
+  AliCaloClusterInfo *caloCluster = 0;
+  for (Int_t iclst=0; iclst<nclsts; iclst++) {  // loop over all clusters
+    clust = fInputEvent->GetCaloCluster(iclst);
+    if (!(clust && clust->IsPHOS() && clust->E()>fgMinClusterEnergy))     { clust=0; continue; }
+    if (!(clust->GetNCells()>fgMinNCells && clust->GetM02()>fgMinM02))    { clust=0; continue; } // To remove exotic clusters
+    if (!(clust->GetDistanceToBadChannel()>fgMinDistToBad))               { clust=0; continue; }
+    clust->GetPosition(position); TVector3 global(position);
+    fPHOSGeo->GlobalPos2RelId(global,relId);
+    if (!IsGoodCaloCluster(relId[0], relId[2], relId[3]))                 { clust=0; continue; } // calo cluster selection 
+    if (relId[0] == 2)                                                    { clust=0; continue; } // !remove module 2
+
+    caloCluster = new AliCaloClusterInfo(clust, esd, fPHOSGeo, vtx);
+//  if (esd && !(caloCluster->LorentzVector().E()>fgMinClusterEnergy)) { delete caloCluster; caloCluster=0; clust=0; continue; } // check again for ESD
+    if (caloCluster->TestCPV(fHeader->MagneticField())) caloCluster->SetPIDBit(BIT(0));          // set CPV bit
+
+    clust = 0;
+
+    new(caloRef[countN++]) AliCaloClusterInfo(*caloCluster);
+
+    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;
+}
diff --git a/PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.h b/PWGGA/PHOSTasks/UserTasks/AliAnalysisTaskSEPHOSpPbPi0.h
new file mode 100644 (file)
index 0000000..68f5367
--- /dev/null
@@ -0,0 +1,83 @@
+#ifndef ALIANALYSISTASKSEPHOSPPBPI0_cxx\r
+#define ALIANALYSISTAKSSEPHOSPPBPI0_cxx\r
+\r
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *\r
+ * See cxx source for full Copyright notice                               */\r
+\r
+//*************************************************************************\r
+// Class AliAnalysisTaskSEPHOSpPbPi0\r
+// AliAnalysisTaskSE for the gamma and pi0 from pPb collision analysis\r
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch\r
+//                   hszhu@iopp.ccnu.edu.cn\r
+//*************************************************************************\r
+\r
+#include "AliAnalysisTaskSE.h"\r
+\r
+class TH2I;\r
+class TList;\r
+class TString;\r
+class TArray;\r
+class TClonesArray;\r
+class AliESDEvent;\r
+class AliPHOSGeoUtils;\r
+class AliPHOSpPbPi0Header;\r
+\r
+class AliAnalysisTaskSEPHOSpPbPi0 : public AliAnalysisTaskSE {\r
+ public:\r
+\r
+  AliAnalysisTaskSEPHOSpPbPi0();\r
+  AliAnalysisTaskSEPHOSpPbPi0(const char *name);\r
+  virtual ~AliAnalysisTaskSEPHOSpPbPi0(); \r
+  \r
+  virtual void UserCreateOutputObjects();\r
+  virtual void UserExec(Option_t *option);\r
+  virtual void Terminate(Option_t *option);\r
+\r
+  void SetPHOSBadMap(Int_t mod,TH2I *hMap)\r
+  {\r
+    if(fPHOSBadMap[mod]) delete fPHOSBadMap[mod];\r
+    fPHOSBadMap[mod]=new TH2I(*hMap);\r
+    printf("Set %s \n",fPHOSBadMap[mod]->GetName());\r
+  }\r
+\r
+  void SetUseMC(Bool_t isMC)                                  { fIsMC          = isMC;                         }\r
+  void SetXBins(const TArrayD& tCent, const TArrayI& tBuffer) { fCentralityBin = tCent; fBufferSize = tBuffer; }\r
+  void SetLogWeight(Float_t logWeight)                  const { AliCaloClusterInfo::SetLogWeight(logWeight);   }\r
+\r
+  static void SetMinNCells(Int_t ncells=2)                    { fgMinNCells                         = ncells;  }\r
+  static void SetMinClusterEnergy(Double_t energy=0.3)        { fgMinClusterEnergy                  = energy;  }\r
+  static void SetMinM02(Double_t m02=0.2)                     { fgMinM02                            = m02;     }\r
+  static void SetMinDistToBad(Double_t dist=2.5)              { fgMinDistToBad                      = dist;    }\r
+\r
+ private:\r
+\r
+  AliAnalysisTaskSEPHOSpPbPi0(const AliAnalysisTaskSEPHOSpPbPi0&);            // not implemented\r
+  AliAnalysisTaskSEPHOSpPbPi0& operator=(const AliAnalysisTaskSEPHOSpPbPi0&); // not implemented\r
+\r
+  void PHOSInitialize(AliESDEvent* const esd);\r
+  void FillCaloClusterInfo(Int_t nclsts, AliESDEvent* const esd);\r
+\r
+  Bool_t IsGoodCaloCluster(Int_t iMod, Int_t cellX, Int_t cellZ);\r
+\r
+  static Int_t    fgMinNCells;\r
+  static Double_t fgMinClusterEnergy;\r
+  static Double_t fgMinM02;\r
+  static Double_t fgMinDistToBad;\r
+\r
+  Bool_t               fIsMC;               // flag of whether the input is MC\r
+  TArrayD              fCentralityBin;      // Centrality bin\r
+  TArrayI              fBufferSize;         // Buffer size for event mixing\r
+\r
+  Int_t                fRunNumber;          // Run Number\r
+  TH2I                *fPHOSBadMap[5];      // Container for PHOS bad channels map\r
+  AliPHOSGeoUtils     *fPHOSGeo;            // PHOS geometry\r
+\r
+  TList               *fEventList[10][10];  // Event list for mixing\r
+  TList               *fList;               // output list of histograms\r
+  AliPHOSpPbPi0Header *fHeader;             // output for info at ev level\r
+  TClonesArray        *fCaloClArr;          // output clones array for Calo clusters\r
+   \r
+  ClassDef(AliAnalysisTaskSEPHOSpPbPi0, 1);\r
+};\r
+\r
+#endif\r
diff --git a/PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.cxx b/PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.cxx
new file mode 100644 (file)
index 0000000..6c30f16
--- /dev/null
@@ -0,0 +1,376 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/////////////////////////////////////////////////////////////
+//
+// class used to extract and store reco info of calo cluster
+//
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch
+//                   hszhu@iopp.ccnu.edu.cn
+/////////////////////////////////////////////////////////////
+
+#include <iostream>
+
+#include <TF1.h>
+#include <TVector3.h>
+
+#include "AliESDEvent.h"
+#include "AliAODTrack.h"
+#include "AliESDtrack.h"
+#include "AliVCluster.h"
+#include "AliESDCaloCluster.h"
+#include "AliESDCaloCells.h"
+#include "AliPHOSGeoUtils.h"
+#include "AliPHOSEsdCluster.h"
+#include "AliPHOSCalibData.h"
+#include "AliCaloClusterInfo.h"
+
+class TObject;
+
+ClassImp(AliCaloClusterInfo) 
+
+Float_t AliCaloClusterInfo::fgLogWeight = 0.07;
+
+//-----------------------------------------------------------------------------
+Double_t NonLinear(Double_t *x, Double_t * /*par*/)
+{ 
+  // a = par[0], b = par[1].
+  // 1+a*exp(-e/b)
+  
+// return 0.0241+1.0504*x[0]+0.000249*x[0]*x[0];
+ return 1.015*(0.0241+1.0504*x[0]+0.000249*x[0]*x[0]);
+  
+}
+
+//-------------------------------------------------------------------------------------------
+AliCaloClusterInfo::AliCaloClusterInfo():
+  TObject(), fLorentzVector(), fModule(0), fNCells(0), fTRUNumber(0), fNTracksMatched(0), fTrackCharge(0), fPIDBit(0),
+  fDistToBad(0.), fEmcCpvDistance(0.), fM02(0.), fM20(0.), fTOF(0.), fTrackDz(0.), fTrackDx(0.), fTrackPt(0.)
+{
+  //
+  // Default constructor
+  //
+} 
+
+//-------------------------------------------------------------------------------------------
+AliCaloClusterInfo::AliCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, AliPHOSGeoUtils* const phosGeo, Double_t vtx[3]):
+  TObject(), fLorentzVector(), fModule(0), fNCells(0), fTRUNumber(0), fNTracksMatched(0), fTrackCharge(0), fPIDBit(0),
+  fDistToBad(0.), fEmcCpvDistance(0.), fM02(0.), fM20(0.), fTOF(0.), fTrackDz(0.), fTrackDx(0.), fTrackPt(0.)
+{
+  //
+  // constructor
+  //
+  this->FillCaloClusterInfo(clust, esd, phosGeo, vtx);
+} 
+
+//-------------------------------------------------------------------------------------------
+AliCaloClusterInfo::AliCaloClusterInfo(const AliCaloClusterInfo &src):
+  TObject(src), fLorentzVector(src.fLorentzVector), fModule(src.fModule), fNCells(src.fNCells), fTRUNumber(src.fTRUNumber),
+  fNTracksMatched(src.fNTracksMatched), fTrackCharge(src.fTrackCharge), fPIDBit(src.fPIDBit), fDistToBad(src.fDistToBad),
+  fEmcCpvDistance(src.fEmcCpvDistance), fM02(src.fM02), fM20(src.fM20), fTOF(src.fTOF), fTrackDz(src.fTrackDz),
+  fTrackDx(src.fTrackDx), fTrackPt(src.fTrackPt)
+{
+  //
+  // copy constructor
+  //
+} 
+
+//-------------------------------------------------------------------------------------------
+AliCaloClusterInfo& AliCaloClusterInfo::operator=(const AliCaloClusterInfo &src)
+{
+  //
+  // assignment constructor
+  //
+  if(&src==this) return *this;
+
+  fLorentzVector   = src.fLorentzVector;
+  fModule          = src.fModule;
+  fNCells          = src.fNCells;
+  fTRUNumber       = src.fTRUNumber;
+  fNTracksMatched  = src.fNTracksMatched;
+  fTrackCharge     = src.fTrackCharge;
+  fPIDBit          = src.fPIDBit;
+  fDistToBad       = src.fDistToBad;
+  fEmcCpvDistance  = src.fEmcCpvDistance;
+  fM02             = src.fM02;
+  fM20             = src.fM20;
+  fTOF             = src.fTOF;
+  fTrackDz         = src.fTrackDz;
+  fTrackDx         = src.fTrackDx;
+  fTrackPt         = src.fTrackPt;
+
+  return *this;
+} 
+
+//-------------------------------------------------------------------------------------------
+AliCaloClusterInfo::~AliCaloClusterInfo()
+{
+  //
+  // destructor
+  //
+}
+
+//-----------------------------------------------------------------------------
+void AliCaloClusterInfo::FillCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, AliPHOSGeoUtils* const phosGeo, Double_t vtx[3])
+{
+  // extract information of calo clusters
+  AliAODTrack *trkAOD = 0x0;
+  AliESDtrack *trkESD = 0x0;
+
+  Int_t   relId[4]    = {0,0,0,0}; // module = relId[0]; cellX = relId[2]; cellZ = relId[3];
+  Float_t position[3] = {0,0,0};
+  clust->GetPosition(position); TVector3 global(position);
+  phosGeo->GlobalPos2RelId(global,relId);
+  fModule         = relId[0];
+  fTRUNumber      = GetTRUNumber(relId[2], relId[3]);
+
+  AliPHOSCalibData *calibData = 0x0;
+  calibData = new AliPHOSCalibData();
+  if (esd) { // TODO recalibration for ESD
+    TVector3 vtxVector(vtx);
+    TF1 *nonLinCorr = new TF1("Non-linear", NonLinear, 0., 40., 0);
+
+    AliPHOSEsdCluster phosClust( *(AliESDCaloCluster*) (clust) );
+    phosClust.Recalibrate(calibData, esd->GetPHOSCells());
+    Reclusterize(&phosClust, phosGeo);               // reclusterize 
+    phosClust.EvalAll(fgLogWeight, vtxVector);       // recalculate all cluster parameters
+    phosClust.SetE(nonLinCorr->Eval(phosClust.E())); // non-linear correction
+
+    TVector3 localPos;
+    const Float_t shiftX[6] = { 0.,-2.3,-2.11,-1.53, 0., 0. };
+    const Float_t shiftZ[6] = { 0.,-0.4, 0.52, 0.8,  0., 0. };
+    phosGeo->Global2Local(localPos, global, fModule);
+    phosGeo->Local2Global(fModule, localPos.X()+shiftX[fModule], localPos.Z()+shiftZ[fModule], global);
+    position[0] = global.X();
+    position[1] = global.Y();
+    position[2] = global.Z();
+    phosClust.SetPosition(position);
+
+    phosClust.GetMomentum(fLorentzVector, vtx);
+
+    //TODO: Check, this may be LHC10h specific:
+    if (fModule==2) fLorentzVector *= 135.5/134.0;
+    if (fModule==3) fLorentzVector *= 135.5/137.2;
+
+    Int_t iESDtrack = clust->GetTrackMatchedIndex();
+    if (iESDtrack>-1) trkESD = esd->GetTrack(iESDtrack);
+    if (!trkESD) { fTrackPt = 0.; fTrackCharge = 0; }
+    else {
+      fTrackPt     = trkESD->Pt();
+      fTrackCharge = trkESD->Charge();
+    }
+  }else {
+    clust->GetMomentum(fLorentzVector, vtx);
+
+    trkAOD = dynamic_cast <AliAODTrack*> (clust->GetTrackMatched(0));
+    if (!trkAOD) { fTrackPt = 0.; fTrackCharge = 0; }
+    else {
+      fTrackPt     = trkAOD->Pt();
+      fTrackCharge = trkAOD->Charge();
+    }
+  }
+
+  fNCells         = clust->GetNCells();
+  fNTracksMatched = clust->GetNTracksMatched();
+  fDistToBad      = clust->GetDistanceToBadChannel();
+  fEmcCpvDistance = clust->GetEmcCpvDistance();
+  fM02            = clust->GetM02();
+  fM20            = clust->GetM20();
+  fTOF            = clust->GetTOF();
+  fTrackDz        = clust->GetTrackDz();
+  fTrackDx        = clust->GetTrackDx();
+  if (TestDisp())                             fPIDBit |= BIT(1); // Disp
+  if (IsInFiducialRegion(/*relId[2], relId[3])*/)) fPIDBit |= BIT(2); // Fiducial
+
+  return;
+}
+
+//-----------------------------------------------------------------------------
+void AliCaloClusterInfo::Reclusterize(AliVCluster *clust, AliPHOSGeoUtils* const phosGeo)
+{
+  // Reclusterize to have continuous cluster
+
+  const Int_t oldMulDigit = clust->GetNCells();
+  Double32_t       *elist = clust->GetCellsAmplitudeFraction();
+  UShort_t         *dlist = clust->GetCellsAbsId();
+
+  Int_t index[oldMulDigit];
+  Bool_t used[oldMulDigit];
+  for (Int_t i=0; i<oldMulDigit; i++) used[i] = 0;
+  Int_t   inClu = 0;
+  Double_t eMax = 0.;
+  for (Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) { // find maximum
+    if (eMax<elist[iDigit]) {
+      eMax     = elist[iDigit];
+      index[0] = iDigit;
+      inClu    = 1;
+    }
+  }
+  if (inClu==0) return;    // empty cluster
+  used[index[0]] = kTRUE ; // mark as used
+  for (Int_t i=0; i<inClu; i++) {
+    for (Int_t iDigit=0; iDigit<oldMulDigit; iDigit++) {
+      if (used[iDigit]) continue; // already used
+      if (AreNeighbors(dlist[index[i]], dlist[iDigit], phosGeo)) {
+        index[inClu] = iDigit;
+        inClu++;
+        used[iDigit] = kTRUE;
+      }
+    }
+  }
+
+  if (inClu==oldMulDigit) return; // no need to modify
+
+  // copy
+  clust->SetNCells(inClu);
+  UShort_t tmpD[oldMulDigit];
+  Double_t tmpE[oldMulDigit];
+  for (Int_t i=0; i<oldMulDigit; i++) {
+    tmpD[i]=dlist[i];
+    tmpE[i]=elist[i];
+  }
+  // change order of digits in list so that
+  // first inClu cells were true ones
+  for (Int_t i=0; i<inClu; i++) {
+    dlist[i]=tmpD[index[i]];
+    elist[i]=tmpE[index[i]];
+  }
+
+  return;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliCaloClusterInfo::TestCPV(Double_t mf)
+{
+  // Parameterization of LHC10h period
+
+  Double_t meanX=0;
+  Double_t meanZ=0.;
+  Double_t sx = TMath::Min(5.4,2.59719e+02*TMath::Exp(-fTrackPt/1.02053e-01)+
+                           6.58365e-01*5.91917e-01*5.91917e-01/((fTrackPt-9.61306e-01)*(fTrackPt-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/(fTrackPt*fTrackPt+1.91456e-02*1.91456e-02)+1.60);
+
+  if(mf<0.) { // Field --
+    meanZ = -0.468318;
+    if (fTrackCharge>0)
+      meanX = TMath::Min(7.3, 3.89994*1.20679*1.20679/(fTrackPt*fTrackPt+1.20679*1.20679)+0.249029+2.49088e+07*TMath::Exp(-fTrackPt*3.33650e+01));
+    else
+      meanX =-TMath::Min(7.7,3.86040*0.912499*0.912499/(fTrackPt*fTrackPt+0.912499*0.912499)+1.23114+4.48277e+05*TMath::Exp(-fTrackPt*2.57070e+01));
+  }
+  else {      // Field ++
+    meanZ = -0.468318;
+    if (fTrackCharge>0)
+      meanX =-TMath::Min(8.0,3.86040*1.31357*1.31357/(fTrackPt*fTrackPt+1.31357*1.31357)+0.880579+7.56199e+06*TMath::Exp(-fTrackPt*3.08451e+01));
+    else
+      meanX = TMath::Min(6.85, 3.89994*1.16240*1.16240/(fTrackPt*fTrackPt+1.16240*1.16240)-0.120787+2.20275e+05*TMath::Exp(-fTrackPt*2.40913e+01));
+  }
+
+  Double_t rx = (fTrackDx-meanX)/sx;
+  Double_t rz = (fTrackDz-meanZ)/sz;
+
+  return TMath::Sqrt(rx*rx+rz*rz) > 2.;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliCaloClusterInfo::TestDisp()
+{
+  Double_t pt = fLorentzVector.E();
+  Double_t l0 = fM02;
+  Double_t l1 = fM20;
+
+  Double_t l0Mean  = 1.53126+9.50835e+06/(1.+1.08728e+07*pt+1.73420e+06*pt*pt);
+  Double_t l1Mean  = 1.12365+0.123770*TMath::Exp(-pt*0.246551)+5.30000e-03*pt;
+  Double_t l0Sigma = 6.48260e-02+7.60261e+10/(1.+1.53012e+11*pt+5.01265e+05*pt*pt)+9.00000e-03*pt;
+  Double_t l1Sigma = 4.44719e-04+6.99839e-01/(1.+1.22497e+00*pt+6.78604e-07*pt*pt)+9.00000e-03*pt;
+  Double_t c       =-0.35-0.550*TMath::Exp(-0.390730*pt);
+  Double_t R2      = 0.5*(l1-l1Mean)*(l1-l1Mean)/l1Sigma/l1Sigma +
+                     0.5*(l0-l0Mean)*(l0-l0Mean)/l0Sigma/l0Sigma +
+                     0.5*c*(l1-l1Mean)*(l0-l0Mean)/l1Sigma/l0Sigma;
+
+  return R2 < 2.5*2.5;
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliCaloClusterInfo::IsInFiducialRegion(/*Int_t cellX, Int_t cellZ*/)
+{
+/*const Int_t edgeX = 2;
+  const Int_t edgeZ = 8;
+  if (cellX >edgeX && cellX<(65-edgeX) && cellZ>edgeZ && cellZ <(57-edgeZ)) return kTRUE; // remove 2 cells in x direction and 8 cells in z direction
+  return kFALSE;*/
+
+  Double_t eta = TMath::Abs(fLorentzVector.Eta());                                // abs eta
+//Double_t phi = (fLorentzVector.Phi())*TMath::RadToDeg(); if (phi<0.) phi+=360.; // phi in degree
+
+  return (Bool_t)(eta<0.13 /*&& ((phi>260.25 && phi<279.75) || (phi>300.25 && phi<319.75))*/);
+}
+
+//-----------------------------------------------------------------------------
+Bool_t AliCaloClusterInfo::AreNeighbors(Int_t id1, Int_t id2, AliPHOSGeoUtils* const phosGeo)
+{
+  // return true if absId are "Neighbors" (adjacent, including diagornaly,)
+  // false if not.
+
+  Int_t relid1[4] ;
+  phosGeo->AbsToRelNumbering(id1, relid1) ;
+
+  Int_t relid2[4] ;
+  phosGeo->AbsToRelNumbering(id2, relid2) ;
+
+  // if inside the same PHOS module
+  if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) {
+    const Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
+    const Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
+
+    // and if diff in both direction is 1 or less
+    if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
+      return kTRUE; // are neighbors
+  }
+
+  // else false
+  return kFALSE;
+}
+
+//-----------------------------------------------------------------------------
+Int_t AliCaloClusterInfo::GetTRUNumber(Int_t cellX, Int_t cellZ)
+{
+  // Return TRU region number for given cell.
+  // cellX: [1-64], cellZ: [1-56]
+  
+  // RCU0: TRU 1,2
+  if (0<cellX && cellX<17) {
+    if (0<cellZ && cellZ<29) return 2;
+    else                     return 1;
+  }
+
+  // RCU1: TRU 3,4
+  if (16<cellX && cellX<33) {
+    if (0<cellZ && cellZ<29) return 4;
+    else                     return 3;
+  }
+
+  // RCU2: TRU 5,6
+  if (32<cellX && cellX<49) {
+    if (0<cellZ && cellZ<29) return 6;
+    else                     return 5;
+  }
+
+  // RCU3: TRU 7,8
+  if (48<cellX && cellX<65) {
+    if (0<cellZ && cellZ<29) return 8;
+    else                     return 7;
+  }
+
+  return -1;
+}
diff --git a/PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.h b/PWGGA/PHOSTasks/UserTasks/AliCaloClusterInfo.h
new file mode 100644 (file)
index 0000000..0b9c067
--- /dev/null
@@ -0,0 +1,86 @@
+#ifndef ALICALOCLUSTERINFO_H
+#define ALICALOCLUSTERINFO_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//***********************************************************
+// Class AliCaloClusterInfo
+// class used to extract and store reco info of calo cluster
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch
+//                   hszhu@iopp.ccnu.edu.cn
+//***********************************************************
+
+#include <TObject.h>
+#include <TString.h>
+#include <TLorentzVector.h>
+
+class AliESDEvent;
+class AliVCluster;
+class AliPHOSGeoUtils;
+
+class AliCaloClusterInfo : public TObject{
+ public:
+  
+  AliCaloClusterInfo();
+  AliCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, AliPHOSGeoUtils* const phosGeo, Double_t vtx[3]);
+  AliCaloClusterInfo(const AliCaloClusterInfo &src);
+  AliCaloClusterInfo& operator=(const AliCaloClusterInfo &src);
+  virtual ~AliCaloClusterInfo();
+
+  TLorentzVector LorentzVector() const{ return fLorentzVector; }
+
+  Int_t    GetModule()         const { return fModule;          }
+  Int_t    GetNCells()         const { return fNCells;          }
+  Int_t    GetTRUNumber()      const { return fTRUNumber;       }
+  Int_t    GetNTracksMatched() const { return fNTracksMatched;  }
+  Short_t  GetTrackCharge()    const { return fTrackCharge;     }
+  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;             }
+  Double_t GetTrackDz()        const { return fTrackDz;         }
+  Double_t GetTrackDx()        const { return fTrackDx;         }
+  Double_t GetTrackPt()        const { return fTrackPt;         }
+
+  void SetPIDBit(UInt_t bit)         { fPIDBit    |= bit;       }
+
+  Bool_t TestCPV(Double_t mf);
+
+  static void SetLogWeight(Float_t logWeight=0.05) { fgLogWeight = logWeight; }
+
+ private:
+
+  void FillCaloClusterInfo(AliVCluster* const clust, AliESDEvent* const esd, AliPHOSGeoUtils* const phosGeo, Double_t vtx[3]);
+  void Reclusterize(AliVCluster *clust, AliPHOSGeoUtils* const phosGeo);
+
+  Bool_t TestDisp();
+  Bool_t AreNeighbors(Int_t id1,Int_t id2, AliPHOSGeoUtils* const phosGeo);
+  Bool_t IsInFiducialRegion(/*Int_t cellX, Int_t cellZ*/);
+  Int_t  GetTRUNumber(Int_t cellX, Int_t cellZ);
+
+  static Float_t fgLogWeight; 
+
+  TLorentzVector fLorentzVector;
+
+  Int_t    fModule;
+  Int_t    fNCells;                      // Number of cells in cluster
+  Int_t    fTRUNumber;                   // TRU Number
+  Int_t    fNTracksMatched;
+  Short_t  fTrackCharge;
+  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
+  Double_t fM02;                         // lambda0
+  Double_t fM20;                         // lambda1
+  Double_t fTOF;
+  Double_t fTrackDz;
+  Double_t fTrackDx;
+  Double_t fTrackPt;
+
+  ClassDef(AliCaloClusterInfo,1);
+};
+
+#endif // #ifdef ALICALOCLUSTERINFO_H
diff --git a/PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx b/PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx
new file mode 100644 (file)
index 0000000..8742acb
--- /dev/null
@@ -0,0 +1,996 @@
+/**************************************************************************
+ * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+///////////////////////////////////////////////////////////////////////////
+//
+// class used to extract ,store info and fill histograms at event level
+//
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch
+//                   hszhu@iopp.ccnu.edu.cn
+///////////////////////////////////////////////////////////////////////////
+
+#include <iostream>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TMath.h>
+#include <TList.h>
+#include <TClonesArray.h>
+
+#include "AliInputEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliAODMCParticle.h"
+#include "AliVEvent.h"
+#include "AliVVertex.h"
+#include "AliVCaloCells.h"
+#include "AliPHOSGeoUtils.h"
+#include "AliAODEvent.h"
+#include "AliESDEvent.h"
+#include "AliCentrality.h"
+
+#include "AliCaloClusterInfo.h"
+#include "AliPHOSpPbPi0Header.h"
+
+class TNamed;
+
+ClassImp(AliPHOSpPbPi0Header)
+
+Bool_t AliPHOSpPbPi0Header::fgIsMC  = kFALSE;
+Int_t  AliPHOSpPbPi0Header::fgNCent = 10;
+
+//_____________________________________________________________________________
+AliPHOSpPbPi0Header::AliPHOSpPbPi0Header() :
+TNamed(),
+fFiredTriggerClass(),
+fSelMask(0),
+fVtxContrsN(0),
+fIsPileupSPD(kFALSE),
+fCentrality(0.),
+fMagneticField(0.)
+{
+  //
+  // default constructor
+  //
+  for (Int_t i=3; i--;) fVtx[i] = 0.;
+}
+
+//_____________________________________________________________________________
+AliPHOSpPbPi0Header::AliPHOSpPbPi0Header(const AliPHOSpPbPi0Header &src) :
+TNamed(),
+fFiredTriggerClass(src.fFiredTriggerClass),
+fSelMask(src.fSelMask),
+fVtxContrsN(src.fVtxContrsN),
+fIsPileupSPD(src.fIsPileupSPD),
+fCentrality(src.fCentrality),
+fMagneticField(src.fMagneticField)
+{
+  //
+  // copy constructor
+  //
+  for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
+}
+
+//_____________________________________________________________________________
+AliPHOSpPbPi0Header& AliPHOSpPbPi0Header::operator=(const AliPHOSpPbPi0Header &src)
+{
+  //
+  // assignment constructor
+  //
+
+  if(&src==this) return *this;
+
+  fFiredTriggerClass            = src.fFiredTriggerClass;
+  fSelMask                      = src.fSelMask;
+  fVtxContrsN                   = src.fVtxContrsN;
+  fIsPileupSPD                  = src.fIsPileupSPD;
+  fCentrality                   = src.fCentrality;
+  fMagneticField                = src.fMagneticField;
+  for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
+
+  return *this;
+}
+
+//_____________________________________________________________________________
+AliPHOSpPbPi0Header::~AliPHOSpPbPi0Header()
+{
+  //
+  // default destructor
+  //
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::SetEventInfo(AliInputEventHandler* const handler)
+{
+  // fill info at event level
+
+  AliVEvent *event = handler->GetEvent();
+  AliAODEvent *aod = dynamic_cast<AliAODEvent*>(event);
+  AliESDEvent *esd = dynamic_cast<AliESDEvent*>(event);
+
+  if (aod) fFiredTriggerClass = aod->GetFiredTriggerClasses();
+  if (esd) fFiredTriggerClass = esd->GetFiredTriggerClasses();
+  fSelMask = handler->IsEventSelected();
+  
+  const AliVVertex *vertex = event->GetPrimaryVertex(); vertex->GetXYZ(fVtx);
+  fVtxContrsN  = vertex->GetNContributors();
+  fIsPileupSPD = (aod && !aod->GetTracklets()) ? event->IsPileupFromSPD(3,0.8,3.,2.,5.) : event->IsPileupFromSPDInMultBins(); //TODO not sure!!!
+//fIsPileupSPD = event->IsPileupFromSPD(3,0.8,3.,2.,5.);
+  
+  AliCentrality *cent = event->GetCentrality();
+  if (cent) fCentrality = cent->GetCentralityPercentileUnchecked("V0M");
+
+  fMagneticField = event->GetMagneticField();
+
+//this->SetTitle(vertex->GetTitle());
+  return;
+}
+
+//_____________________________________________________________________________
+Bool_t AliPHOSpPbPi0Header::IspAVertexOK(AliAODEvent* const event)
+{
+  // check p-A collection vertex
+  const AliAODVertex* trkVtx = event->GetPrimaryVertex();   Double_t zvtx = trkVtx->GetZ();
+  if (!trkVtx || trkVtx->GetNContributors()<1)                                            return kFALSE;
+  if (!((TString)trkVtx->GetTitle()).Contains("VertexerTracks"))                          return kFALSE;
+  const AliAODVertex* spdVtx = event->GetPrimaryVertexSPD();
+  if (spdVtx->GetNContributors()<1)                                                       return kFALSE;
+  Double_t cov[6]={0}; spdVtx->GetCovarianceMatrix(cov);
+  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() - zvtx)>0.5)                                              return kFALSE;
+  if (TMath::Abs(zvtx>10.))                                                               return kFALSE;
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+Bool_t AliPHOSpPbPi0Header::IspAVertexOK(AliESDEvent* const event)
+{
+  // check p-A collection vertex
+  const AliESDVertex* trkVtx = event->GetPrimaryVertex();   Double_t zvtx = trkVtx->GetZ();
+  if (!trkVtx || trkVtx->GetNContributors()<1)                                            return kFALSE;
+  if (!((TString)trkVtx->GetTitle()).Contains("VertexerTracks"))                          return kFALSE;
+  const AliESDVertex* spdVtx = event->GetPrimaryVertexSPD();
+  if (spdVtx->GetNContributors()<1)                                                       return kFALSE;
+  Double_t cov[6]={0}; spdVtx->GetCovarianceMatrix(cov);
+  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() - zvtx)>0.5)                                              return kFALSE;
+  if (TMath::Abs(zvtx>10.))                                                               return kFALSE;
+
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistograms(TList *list)
+{
+  // create output histos of partcorr analysis according to the MC flag
+
+  this->CreateHistosEvnH(list);
+  this->CreateHistosCaloCellsQA(list);
+  this->CreateHistosCaloCluster(list);
+  this->CreateHistosPi0(list);
+  this->CreateHistosMixPi0(list);
+  if (fgIsMC) this->CreateHistosMC(list);
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosEvnH(TList *list)
+{
+  // 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  ,         100 , 200 ,     200  };
+  Double_t xlow[nhs] = {        -2.5,           0., -10.,     -10. };
+  Double_t  xup[nhs] = {       199.5,         100.,  10.,      10. };
+
+  TH1D *histo = 0;
+  for (Int_t i=0; i<nhs; i++) {
+    char *hName = Form("hEvnH_%s", tName[i].Data());
+    histo = new TH1D(hName, hName, nbins[i], xlow[i], xup[i]);
+    histo->Sumw2(); list->Add(histo); histo = 0;
+  }
+  TH1::AddDirectory(oldStatus);
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosCaloCellsQA(TList *list)
+{
+  // 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;
+  for (Int_t i=0; i<nMod; i++) {
+    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;
+  }
+
+  TString tName[nhs] = { "CentNCells", "Nxz", "Exz" };
+  Int_t  nbinsx[nhs] = {          10 ,   64 ,  64   };
+  Double_t xlow[nhs] = {           0.,   0.5,   0.5 };
+  Double_t  xup[nhs] = {         100.,  64.5,  64.5 };
+  Int_t  nbinsy[nhs] = {        1000 ,   56 ,  56   };
+  Double_t ylow[nhs] = {           0.,   0.5,   0.5 };
+  Double_t  yup[nhs] = {        1000.,  56.5,  56.5 };
+
+  for (Int_t i=0; i<nhs; i++) {
+    if (i == 0) {
+      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;
+    } else {
+      for (Int_t j=0; j<nMod; j++) {
+        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;
+      }
+    }
+  }
+
+  TH1::AddDirectory(oldStatus1);
+  TH2::AddDirectory(oldStatus2);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosCaloCluster(TList *list)
+{
+  // 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,  kFiducial,  kBoth, kBothPlusFiducial, kPIDs                       };   // PID
+  TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Fiducial", "Both", "BothPlusFiducial"                             };
+                         // { 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             };
+  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.            };
+
+  TH1I *hncells = 0;
+  TH1D *histo1  = 0;
+  TH2D *histo2  = 0;
+  TString hName;
+
+  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;
+
+  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;
+
+  for (Int_t icent=0; icent<fgNCent; icent++) {
+    hName  = Form("hCaloCluster_%s_cent%d", tName[kPtClu].Data(), icent);
+    histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
+    histo1->Sumw2(); list->Add(histo1); histo1 = 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;
+
+    if (iPID>kFiducial) {
+      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;
+      }
+
+      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;
+      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());
+        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;
+      }
+    }
+  }
+
+  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;
+
+    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;
+
+    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;
+    }
+
+    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;
+    }
+  }
+
+  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;
+  for (Int_t i=0; i<kVarsClu-3; i++) {
+    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;
+  }
+
+  TH1::AddDirectory(oldStatus1);
+  TH2::AddDirectory(oldStatus2);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosPi0(TList *list)
+{
+  // 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,  kFiducial,  kBoth, kBothPlusFiducial, kPIDs   };  // PID
+  TString namePID[kPIDs]  = {  "All", "Cpv", "Disp", "Fiducial", "Both", "BothPlusFiducial"         };
+                         // { kPtPi0, kEtaPi0, kPhiPi0, kAsyPi0,   kAnglePi0, kInvMassPi0, kVarsPi0 };  // pi0
+  TString tName[kVarsPi0] = {   "Pt",   "Eta",   "Phi",   "Asy",     "Angle",   "InvMass"           };
+  Int_t    bins[kVarsPi0] = {   500 ,  300   ,   120  ,    100 ,        180 ,       1000            };
+  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.           };
+
+  TH2D   *histo = 0;
+  TString hName;
+
+  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;
+  for (Int_t i=0; i<kVarsPi0-2; i++) {
+    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;
+  }
+
+  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;
+    if (iPID>kFiducial) {
+      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;
+      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;
+      }
+
+      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());
+        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;
+      }
+    }
+    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;
+    }
+  }
+
+  TH2::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosMixPi0(TList *list)
+{
+  // 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,  kFiducial,  kBoth, kBothPlusFiducial, kPIDs  };  // PID
+  TString  namePID[kPIDs]    = { "All", "Cpv", "Disp", "Fiducial", "Both", "BothPlusFiducial"        };
+                             //{ 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.                   };
+
+  TH2D   *histo = 0;
+  TString hName;
+
+  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;
+  
+  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;
+    if (iPID>kFiducial) {
+      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;
+      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;
+      } 
+    }
+    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;
+    }
+  }
+
+  TH2::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::CreateHistosMC(TList *list)
+{
+  // create Histograms for MC
+  if (!list) list = new TList();
+  list->SetOwner();
+
+  Bool_t oldStatus = TH2::AddDirectoryStatus(); TH2::AddDirectory(kFALSE);
+  const Int_t method = 2;
+  char *mName[method]= {"IsPrimary", "RCut"};
+  const Int_t cut    = 3;
+  char *cName[cut]   = {"WideY", "NarrY", "Acc"};
+  const Int_t type   = 3;
+  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.           };
+
+  TH2D   *histo = 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());
+    histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kVertexMC], xMin[kVertexMC], xMax[kVertexMC]);
+    histo->Sumw2(); list->Add(histo); histo = 0;
+    hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kRapidityMC].Data());
+    histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kRapidityMC], xMin[kRapidityMC], xMax[kRapidityMC]);
+    histo->Sumw2(); list->Add(histo); histo = 0;
+    hName  = Form("hMC%s_%s%s_%s", pName[iType], tName[kPtMC].Data(), tName[kVertexMC].Data(), mName[0]);
+    histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kVertexMC], xMin[kVertexMC], xMax[kVertexMC]);
+    histo->Sumw2(); list->Add(histo); histo = 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%s_cent%d_%s", pName[iType], tName[i+1].Data(), tName[kWeightMC].Data(), iCent, mName[iMethod]);
+          histo = new TH2D(hName.Data(), hName.Data(), bins[i+1], xMin[i+1], xMax[i+1], bins[kWeightMC], xMin[kWeightMC], xMax[kWeightMC]);
+          histo->Sumw2(); list->Add(histo); histo = 0;
+
+        }
+        for (Int_t iCut=0; iCut<cut; iCut++) {
+          hName  = Form("hMC%s_%s%s_cent%d_%s_%s", pName[iType], tName[kPtMC].Data(), tName[kWeightMC].Data(), iCent, cName[iCut], mName[iMethod]);
+          histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kWeightMC], xMin[kWeightMC], xMax[kWeightMC]);
+          histo->Sumw2(); list->Add(histo); histo = 0;
+        }
+      }
+    }
+  }
+
+
+  TH2::AddDirectory(oldStatus);
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosEvnH(TList *list)
+{
+  // fill histograms at event level according to event selection cuts
+
+  if (!list) return;
+
+  const Int_t nhs      = 3;
+  TString tName[nhs+1] = { "VtxNcontr", "Centrality",    "Vz", "Pileup" };
+  Double_t dist[nhs]   = { fVtxContrsN,  fCentrality, fVtx[2]           };
+  for (Int_t i=nhs; i--;) ((TH1D*)list->FindObject(Form("hEvnH_%s",tName[i].Data())))->Fill(dist[i]);
+  if (fIsPileupSPD)       ((TH1D*)list->FindObject(Form("hEvnH_%s",tName[nhs].Data())))->Fill(dist[nhs-1]);
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosCaloCellsQA(TList *list, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo)
+{
+  // fill QA histograms for calocells according to evnet selection cuts
+
+  if (!list) 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());
+  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]);
+    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]);
+  }
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosCaloCluster(TList *list, TClonesArray* const caloClArr, Int_t cent)
+{
+  // fill histograms for calocluster according to evnet && calocluster candidates selection cuts
+
+  if (!list)  return;
+                         // {       kAll,   kCpv,  kDisp,  kFiducial,           kBoth,   kBothPlusFiducial,   kPIDs    };   // PID
+  TString namePID[kPIDs]  = {      "All",  "Cpv", "Disp", "Fiducial",          "Both",   "BothPlusFiducial"            };
+  UInt_t  PIDBIT[kPIDs]   = { 0x00000000, BIT(0), BIT(1),     BIT(2), (BIT(0)|BIT(1)), (BIT(0)|BIT(1)|BIT(2))          };
+                         // { kPtClu, kEtaClu, kPhiClu, kM02Clu, kM20Clu, kTOFClu, kNCellsClu, kNClustersClu, kVarsClu };   // clusters
+  TString tName[kVarsClu] = {   "Pt",   "Eta",   "Phi",   "M02",   "M20",   "TOF",   "NCells",   "NClusters"           };
+
+  Int_t entries = caloClArr->GetEntries();
+  Int_t iMod = 0, iTRU = 0, nclusters[3] = { 0, 0, 0 };
+  UInt_t pidBit = 0;
+  AliCaloClusterInfo *caloCluster = 0;
+  TLorentzVector gamma;
+  for(Int_t i=0; i<entries; i++) { // Loop calo cluster
+    caloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
+    iMod   = caloCluster->GetModule();   nclusters[iMod-1]++;
+    iTRU   = caloCluster->GetTRUNumber();
+    pidBit = caloCluster->GetPIDBit();
+    gamma  = caloCluster->LorentzVector();
+    Double_t vars[kVarsClu-1] = { gamma.E(),
+                                  gamma.Eta(),
+                                  (gamma.Phi()+TMath::TwoPi()),
+                                  caloCluster->GetM02(),
+                                  caloCluster->GetM20(),
+                                  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]);
+    ((TH1D*)list->FindObject(Form("hCaloCluster_%s_cent%d", tName[kPtClu].Data(), cent)))->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]);
+        if (iPID>kFiducial) {
+          ((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(),
+                                                                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*)list->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]);
+    }
+  } // End loop calo cluster
+
+  ((TH1I*)list->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]);
+  }
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosPi0(TList *list, TClonesArray* const caloClArr, Int_t cent)
+{
+  // fill histograms for pi0 according to evnet && pi0 candidates selection cuts
+
+  if (!list) return;
+                         // {       kAll,   kCpv,  kDisp,  kFiducial,         kBoth, kBothPlusFiducial,      kPIDs     };   // PID
+  TString namePID[kPIDs] =  {      "All",  "Cpv", "Disp", "Fiducial",        "Both", "BothPlusFiducial"                };
+  UInt_t  PIDBIT[kPIDs]  =  { 0x00000000, BIT(0), BIT(1),     BIT(2), BIT(0)|BIT(1), BIT(0)|BIT(1)|BIT(2)              };
+                         // { kPtPi0, kEtaPi0, kPhiPi0, kAsyPi0, kAnglePi0, kInvMassPi0, kVarsPi0                      };   // pi0
+  TString tName[kVarsPi0] = {   "Pt",   "Eta",   "Phi",   "Asy",   "Angle",   "InvMass"                                };
+
+  Int_t entries = caloClArr->GetEntries();
+  Int_t iMod = 0, jMod = 0, srcMod = 0;
+  UInt_t iPIDBit = 0, jPIDBit = 0;
+  Double_t vars[kVarsPi0];
+  AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
+  TLorentzVector iGamma, jGamma;
+  for(Int_t i=0; i<entries-1; i++) { // Loop calo cluster i
+    iCaloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
+    iMod    = iCaloCluster->GetModule();
+    iPIDBit = iCaloCluster->GetPIDBit();
+    iGamma  = iCaloCluster->LorentzVector();
+    
+    for (Int_t j=0; j<entries; j++) { // Loop calo cluster j
+      jCaloCluster = (AliCaloClusterInfo*)caloClArr->At(j);
+      jMod    = jCaloCluster->GetModule();
+      jPIDBit = jCaloCluster->GetPIDBit();
+      jGamma  = jCaloCluster->LorentzVector();
+      if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
+
+      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
+      else if (iMod==2 && jMod==3)  srcMod = 23; // Mod23
+
+      vars[kPtPi0]      = (iGamma+jGamma).E();
+      vars[kEtaPi0]     = (iGamma+jGamma).Eta();
+      vars[kPhiPi0]     = (iGamma+jGamma).Phi() + TMath::TwoPi();
+      vars[kInvMassPi0] = (iGamma+jGamma).M();
+      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]);
+      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]);
+      }
+
+      for (Int_t iPID=0; iPID<kPIDs; iPID++) {
+        if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
+          ((TH2D*)list->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(),
+                                                              namePID[iPID].Data(), srcMod)))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
+          if (iPID>kFiducial) {
+            ((TH2D*)list->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_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]);
+            }
+          }
+        }
+      }
+      jCaloCluster = 0;
+    }// End loop j
+    iCaloCluster = 0;
+  } // End loop i
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *list, TClonesArray* const iCaloClArr, TList* const eventList, Int_t cent)
+{
+  // fill histograms for mixed pi0 according to evnet && pi0 candidates selection cuts
+
+  if (!list) return;
+                            // {       kAll,   kCpv,  kDisp,  kFiducial,         kBoth, kBothPlusFiducial,     kPIDs   };   // PID
+  TString namePID[kPIDs]    =  {      "All",  "Cpv", "Disp", "Fiducial",        "Both", "BothPlusFiducial"             };
+  UInt_t  PIDBIT[kPIDs]     =  { 0x00000000, BIT(0), BIT(1),     BIT(2), BIT(0)|BIT(1), BIT(0)|BIT(1)|BIT(2)           };
+                            // { kPtMixPi0, kEtaMixPi0, kPhiMixPi0, kInvMassMixPi0, kVarsMixPi0                        };   // Mix pi0
+  TString tName[kVarsMixPi0] = {      "Pt",      "Eta",      "Phi",      "InvMass"                                     };
+
+  Int_t iEntries = iCaloClArr->GetEntries();
+  Int_t iMod = 0, jMod = 0, srcMod = 0;
+  UInt_t iPIDBit = 0, jPIDBit = 0;
+  Double_t vars[kVarsMixPi0];
+  AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
+  TLorentzVector iGamma, jGamma;
+  TClonesArray *jCaloClArr = 0;
+  for (Int_t i=0; i<iEntries; i++) { // Loop calo cluster i
+    iCaloCluster=(AliCaloClusterInfo*)iCaloClArr->At(i);
+    iMod    = iCaloCluster->GetModule();
+    iPIDBit = iCaloCluster->GetPIDBit();
+    iGamma  = iCaloCluster->LorentzVector();
+
+    Int_t nev = eventList->GetSize();
+    for(Int_t ev=0; ev<nev; ev++){ // Loop events for mixing
+      jCaloClArr     = static_cast<TClonesArray*>(eventList->At(ev));
+      Int_t jEntries = jCaloClArr->GetEntries();
+      for(Int_t j=0; j<jEntries; j++){ // Loop calo cluster j
+        jCaloCluster=(AliCaloClusterInfo*)jCaloClArr->At(j);
+        jMod    = jCaloCluster->GetModule();   if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
+        jPIDBit = jCaloCluster->GetPIDBit();
+        jGamma  = jCaloCluster->LorentzVector();
+
+        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) || (iMod==2 && jMod==1))  srcMod = 12; // Mod12
+        else if ((iMod==2 && jMod==3) || (iMod==3 && jMod==2))  srcMod = 23; // Mod23
+
+        vars[kPtMixPi0]      = (iGamma+jGamma).E();
+        vars[kEtaMixPi0]     = (iGamma+jGamma).Eta();
+        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]);
+        for (Int_t iPID=0; iPID<kPIDs; iPID++) {
+          if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
+            ((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]);
+            if (iPID>kFiducial) {
+              ((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_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
+                                                                      namePID[iPID].Data(), cent)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
+            }
+          }
+        }
+        jCaloCluster = 0;
+      } // End loop j
+      jCaloClArr = 0;
+    } // End loop events
+    iCaloCluster = 0;
+  } // End loop i
+
+  return;
+}
+
+//_____________________________________________________________________________
+void AliPHOSpPbPi0Header::FillHistosMC(TList *list, AliMCEvent* const mcEvent, AliPHOSGeoUtils* const phosGeo, Int_t cent)
+{
+  // fill histograms for AOD MC
+  if (!list) return;
+
+  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
+    Double_t pt     = pMC->Pt();
+    Double_t r      = TMath::Sqrt(TMath::Power(pMC->Xv(), 2) + TMath::Power(pMC->Yv(), 2));
+    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);
+    AliAODMCParticle *gamma1=0x0, *gamma2=0x0;
+    if (r<1.) {
+      ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_RapidityWeight_cent%d_RCut", pName.Data(), cent)))->Fill(y, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_PhiWeight_cent%d_RCut", pName.Data(), cent)))->Fill(phi, weight);
+      if (y<1.)   ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_WideY_RCut", pName.Data(), cent)))->Fill(pt, weight); 
+      if (y<0.12) ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_NarrY_RCut", pName.Data(), cent)))->Fill(pt, weight); 
+      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
+        if (IsHitPHOS(gamma1, phosGeo) && IsHitPHOS(gamma2, phosGeo))
+          ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_Acc_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      }
+    }
+    Int_t mother = pMC->GetMother(); if (mother<0) continue;
+    Bool_t isPrimary = ((AliAODMCParticle*)mcEvent->GetTrack(mother))->IsPhysicalPrimary();
+    if (isPrimary) {
+      ((TH2D*)list->FindObject(Form("hMC%s_PtVertex_IsPrimary", pName.Data())))->Fill(pt, r);
+      ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_RapidityWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(y, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_PhiWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(phi, weight);
+      if (y<1.)   ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_WideY_IsPrimary", pName.Data(), cent)))->Fill(pt, weight); 
+      if (y<0.12) ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_NarrY_IsPrimary", pName.Data(), cent)))->Fill(pt, weight); 
+      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
+        if (IsHitPHOS(gamma1, phosGeo) && IsHitPHOS(gamma2, phosGeo))
+          ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_Acc_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      }
+    }
+  }
+
+  return;
+}
+
+//_____________________________________________________________________________
+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;
+
+  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);
+
+    TParticle *gamma1 = 0x0, *gamma2 = 0x0;
+    if (r<1.) {
+      ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_RapidityWeight_cent%d_RCut", pName.Data(), cent)))->Fill(y, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_PhiWeight_cent%d_RCut", pName.Data(), cent)))->Fill(phi, weight);
+      if (y<0.5)  ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_WideY_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      if (y<0.13) ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_NarrY_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      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
+        if (IsHitPHOS(gamma1, phosGeo) && IsHitPHOS(gamma2, phosGeo))
+          ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_Acc_RCut", pName.Data(), cent)))->Fill(pt, weight);
+      }
+    }
+    if (pMC->IsPrimary()) {
+      ((TH2D*)list->FindObject(Form("hMC%s_PtVertex_IsPrimary", pName.Data())))->Fill(pt, r);
+      ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_RapidityWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(y, weight);
+      ((TH2D*)list->FindObject(Form("hMC%s_PhiWeight_cent%d_IsPrimary", pName.Data(), cent)))->Fill(phi, weight);
+      if (y<0.5)  ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_WideY_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      if (y<0.13) ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_NarrY_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      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
+        if (IsHitPHOS(gamma1, phosGeo) && IsHitPHOS(gamma2, phosGeo))
+          ((TH2D*)list->FindObject(Form("hMC%s_PtWeight_cent%d_Acc_IsPrimary", pName.Data(), cent)))->Fill(pt, weight);
+      }
+    }
+
+
+  }
+
+  return;
+}
+
+//________________________________________________________________________
+Double_t AliPHOSpPbPi0Header::PrimaryParticleWeight(Int_t pdg, Double_t pt)
+{
+
+  Int_t type=0 ;
+  if (pdg == 111 || TMath::Abs(pdg)==211) type = 1;
+  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.;
+}
+
+//________________________________________________________________________
+Bool_t AliPHOSpPbPi0Header::IsHitPHOS(AliAODMCParticle* const pMC, AliPHOSGeoUtils* const phosGeo)
+{ 
+
+  Int_t mod=0; 
+  Double_t x=0., z=0.;
+  Double_t vtx[3]; pMC->XvYvZv(vtx);
+  Double_t theta = pMC->Theta();
+  Double_t phi   = pMC->Phi();
+    
+  return (phosGeo->ImpactOnEmc(vtx, theta, phi, mod, z, x) && mod!=2);
+}   
+
+//________________________________________________________________________
+Bool_t AliPHOSpPbPi0Header::IsHitPHOS(TParticle* const pMC, AliPHOSGeoUtils* const phosGeo)
+{
+
+  Int_t mod=0;
+  Double_t x=0., z=0.;
+  Double_t vtx[3] = { pMC->Vx(), pMC->Vy(), pMC->Vz() };
+  Double_t theta = pMC->Theta();
+  Double_t phi   = pMC->Phi();
+
+  return (phosGeo->ImpactOnEmc(vtx, theta, phi, mod, z, x) && mod!=2);
+}
diff --git a/PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.h b/PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.h
new file mode 100644 (file)
index 0000000..519207f
--- /dev/null
@@ -0,0 +1,99 @@
+#ifndef ALIPHOSPPBPI0HEADER_H
+#define ALIPHOSPPBPI0HEADER_H
+
+/* Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+//*************************************************************************
+// Class AliPHOSpPbPi0Header
+// class used to extract ,store info and fill histograms at event level
+// Author: H-S. Zhu, hongsheng.zhu@cern.ch
+//                   hszhu@iopp.ccnu.edu.cn
+//*************************************************************************
+
+#include <TNamed.h>
+#include <TString.h>
+
+class TList;
+class TParticle;
+class TClonesArray;
+
+class AliInputEventHander;
+class AliMCEvent;
+class AliStack;
+class AliAODMCParticle;
+class AliVCaloCells;
+class AliPHOSGeoUtils;
+
+class AliCaloClusterInfo;
+
+class AliPHOSpPbPi0Header : public TNamed {
+ public :
+
+  AliPHOSpPbPi0Header();
+  AliPHOSpPbPi0Header(const AliPHOSpPbPi0Header &src);
+  AliPHOSpPbPi0Header& operator=(const AliPHOSpPbPi0Header &src);
+  ~AliPHOSpPbPi0Header();
+
+  void     GetXYZ(Double_t *vtx)   const { for (Int_t i=3; i--;)   vtx[i]=fVtx[i]; }
+  Double_t Vx()                    const { return fVtx[0];                         }
+  Double_t Vy()                    const { return fVtx[1];                         }
+  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   IsPileupSPD()           const { return fIsPileupSPD;                    }
+  Float_t  Centrality()            const { return fCentrality;                     }
+  Double_t MagneticField()         const { return fMagneticField;                  }
+
+  Bool_t   IspAVertexOK(AliAODEvent* const event);
+  Bool_t   IspAVertexOK(AliESDEvent* const event);
+
+  void SetEventInfo(AliInputEventHandler* const handler);
+
+  void CreateHistograms(TList *list);
+  void FillHistosEvnH(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 SetNCent(Int_t ncent=10)    { fgNCent = ncent; }
+
+ private :
+
+  void CreateHistosEvnH(TList *list);
+  void CreateHistosCaloCellsQA(TList *list);
+  void CreateHistosCaloCluster(TList *list);
+  void CreateHistosPi0(TList *list);
+  void CreateHistosMixPi0(TList *list);
+  void CreateHistosMC(TList *list);
+
+  Double_t PrimaryParticleWeight(Int_t pdg, Double_t pt);
+  Bool_t   IsHitPHOS(AliAODMCParticle* const pMC, AliPHOSGeoUtils* const phosGeo);
+  Bool_t   IsHitPHOS(TParticle* const pMC, AliPHOSGeoUtils* const phosGeo);
+
+  static Bool_t fgIsMC;            // flag to use MC
+  static Int_t  fgNCent;           // # of centrality bins
+
+  enum { kAll,      kCpv,       kDisp,       kFiducial,      kBoth,      kBothPlusFiducial, 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
+
+  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   fIsPileupSPD;                  // is Pileup from SPD
+  Float_t  fCentrality;                   // event certrality
+  Double_t fMagneticField;                // magnetic field
+
+  ClassDef(AliPHOSpPbPi0Header, 1)
+};
+
+#endif
diff --git a/PWGGA/PHOSTasks/UserTasks/macros/AddTaskPHOSpPbPi0.C b/PWGGA/PHOSTasks/UserTasks/macros/AddTaskPHOSpPbPi0.C
new file mode 100644 (file)
index 0000000..4830b17
--- /dev/null
@@ -0,0 +1,56 @@
+AliAnalysisTaskSEPHOSpPbPi0* AddTaskPHOSpPbPi0(Bool_t isMCtruth=kFALSE, UInt_t triggerMask = AliVEvent::kMB, Bool_t isBadMap=kFALSE)
+{
+  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
+  if (!mgr) {
+    ::Error("AddTaskPHOSpPbPi0", "No analysis manager to connect to.");
+    return NULL;
+  }
+  TString type = mgr->GetInputEventHandler()->GetDataType();
+  if (!type.Contains("ESD") && !type.Contains("AOD")) {
+    ::Error("AddTaskPHOSpPbPi0", "PHOSpPbPi0 task needs the manager to have an ESD or AOD input handler.");
+    return NULL;
+  }
+  if (isMCtruth && type.Contains("ESD")) {
+    AliMCEventHandler *mcH = mgr->GetMCtruthEventHandler();
+    if (!mcH) {
+      ::Error("AddTaskPHOSpPbPi0", "PHOSpPbPi0 task needs the manager to have an MC evnet handler.");
+      return NULL;
+    }
+  }
+
+  const Int_t nBins = 6;
+  Double_t cBin[nBins+1] = {0., 20., 40., 60., 80., 90., 100.}; TArrayD tCent(nBins+1, cBin);
+  Int_t    buffer[nBins] = { 40,  80,  80,  100,  100,  100  }; TArrayI tBuffer(nBins, buffer);
+
+  AliAnalysisTaskSEPHOSpPbPi0 *task = new AliAnalysisTaskSEPHOSpPbPi0("TaskPHOSpPbPi0");
+  task->SetUseMC(isMCtruth);
+  task->SetXBins(tCent, tBuffer);
+  task->SetLogWeight(0.07);
+  task->SetMinNCells(2);
+  task->SetMinClusterEnergy(0.3);
+  task->SetMinM02(0.2);
+  task->SetMinDistToBad(2.5);
+  task->SetDebugLevel(-1);
+
+  // 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);
+    }
+  }
+
+  task->SelectCollisionCandidates(triggerMask); 
+  mgr->AddTask(task);
+
+  mgr->ConnectInput(task, 0, mgr->GetCommonInputContainer());
+
+  AliAnalysisDataContainer *coutput = mgr->CreateContainer("histosPHOS", TList::Class(), AliAnalysisManager::kOutputContainer, "histosPHOS.root");
+  mgr->ConnectOutput(task, 1, coutput);
+
+  return task;
+}
index 31c3432..fd8aca5 100644 (file)
@@ -34,4 +34,9 @@
 // Omega3pi
 #pragma link C++ class AliAnalysisTaskOmegaPi0PiPi+;
 
+// UserTasks
+#pragma link C++ class AliCaloClusterInfo+;
+#pragma link C++ class AliPHOSpPbPi0Header+;
+#pragma link C++ class AliAnalysisTaskSEPHOSpPbPi0+;
+
 #endif