PWGGA/PHOSTasks/UserTasks/AliPHOSpPbPi0Header.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2006, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 //
  18 // class used to extract ,store info and fill histograms at event level
  19 //
  20 // Author: H-S. Zhu, hongsheng.zhu@cern.ch
  21 //                   hszhu@iopp.ccnu.edu.cn
  22 ///////////////////////////////////////////////////////////////////////////
  23
  24 #include <iostream>
  25
  26 #include <TH1.h>
  27 #include <TH2.h>
  28 #include <TMath.h>
  29 #include <TList.h>
  30 #include <TClonesArray.h>
  31
  32 #include "AliInputEventHandler.h"
  33 #include "AliMCEvent.h"
  34 #include "AliStack.h"
  35 #include "AliAODMCParticle.h"
  36 #include "AliVEvent.h"
  37 #include "AliVVertex.h"
  38 #include "AliAnalysisUtils.h"
  39 #include "AliVCaloCells.h"
  40 #include "AliPHOSGeoUtils.h"
  41 #include "AliCentrality.h"
  42
  43 #include "AliCaloClusterInfo.h"
  44 #include "AliPHOSpPbPi0Header.h"
  45
  46 class TNamed;
  47
  48 ClassImp(AliPHOSpPbPi0Header)
  49
  50 Bool_t   AliPHOSpPbPi0Header::fgIsMC           = kFALSE;
  51 Bool_t   AliPHOSpPbPi0Header::fgUseFiducialCut = kFALSE;
  52 Int_t    AliPHOSpPbPi0Header::fgNCent          = 10;
  53 Double_t AliPHOSpPbPi0Header::fgCuts[3]        = { 10., 0., 100. };
  54
  55 //_____________________________________________________________________________
  56 AliPHOSpPbPi0Header::AliPHOSpPbPi0Header() :
  57   TNamed(), fFiredTriggerClass(), fSelMask(0), fIsVertexOK(kFALSE), fIsPileup(kFALSE), fCentrality(0.)
  58 {
  59   //
  60   // default constructor
  61   //
  62   for (Int_t i=3; i--;) fVtx[i] = 0.;
  63 }
  64
  65 //_____________________________________________________________________________
  66 AliPHOSpPbPi0Header::AliPHOSpPbPi0Header(const AliPHOSpPbPi0Header &src) :
  67   TNamed(), fFiredTriggerClass(src.fFiredTriggerClass), fSelMask(src.fSelMask), fIsVertexOK(src.fIsVertexOK),
  68   fIsPileup(src.fIsPileup), fCentrality(src.fCentrality)
  69 {
  70   //
  71   // copy constructor
  72   //
  73   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
  74 }
  75
  76 //_____________________________________________________________________________
  77 AliPHOSpPbPi0Header& AliPHOSpPbPi0Header::operator=(const AliPHOSpPbPi0Header &src)
  78 {
  79   //
  80   // assignment constructor
  81   //
  82
  83   if(&src==this) return *this;
  84
  85   fFiredTriggerClass            = src.fFiredTriggerClass;
  86   fSelMask                      = src.fSelMask;
  87   fIsVertexOK                   = src.fIsVertexOK;
  88   fIsPileup                     = src.fIsPileup;
  89   fCentrality                   = src.fCentrality;
  90   for (Int_t i=3; i--;) fVtx[i] = src.fVtx[i];
  91
  92   return *this;
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliPHOSpPbPi0Header::~AliPHOSpPbPi0Header()
  97 {
  98   //
  99   // default destructor
 100   //
 101 }
 102
 103 //_____________________________________________________________________________
 104 void AliPHOSpPbPi0Header::SetEventInfo(AliInputEventHandler* const handler)
 105 {
 106   // fill info at event level
 107
 108   AliVEvent *event = handler->GetEvent();
 109
 110   fSelMask      = handler->IsEventSelected();
 111   fIsVertexOK   = CheckEventVertex(event);
 112
 113   AliCentrality *cent = event->GetCentrality();
 114   if (cent) fCentrality = cent->GetCentralityPercentile("V0A");
 115   else      fCentrality = -1.;
 116
 117 //this->SetTitle(vertex->GetTitle());
 118   return;
 119 }
 120
 121 //_____________________________________________________________________________
 122 Bool_t AliPHOSpPbPi0Header::IsSelected()
 123 {
 124   // select event according to the event selection cuts
 125   if (!fIsVertexOK)                                    return kFALSE; // pA vertex cut
 126   if (TMath::Abs(fVtx[2])>fgCuts[0])                   return kFALSE; // Vtxz cut
 127   if (fCentrality<fgCuts[1] || fCentrality>fgCuts[2])  return kFALSE; // centrality selection
 128
 129   return kTRUE;
 130 }
 131
 132 //_____________________________________________________________________________
 133 Bool_t AliPHOSpPbPi0Header::CheckEventVertex(AliVEvent* const event)
 134 {
 135   // check event vertex
 136
 137   // set event basic info
 138   fFiredTriggerClass       = event->GetFiredTriggerClasses();
 139   const AliVVertex* trkVtx = event->GetPrimaryVertex();
 140   if (!trkVtx)                               return kFALSE;
 141   trkVtx->GetXYZ(fVtx);
 142
 143   AliAnalysisUtils *utils = new AliAnalysisUtils();
 144   if (!utils->IsVertexSelected2013pA(event)) return kFALSE;
 145   fIsPileup = utils->IsPileUpEvent(event);
 146
 147   return kTRUE;
 148 }
 149
 150 //_____________________________________________________________________________
 151 void AliPHOSpPbPi0Header::CreateHistograms(TList *listQA, TList *listRD, TList *listMC)
 152 {
 153   // create output histos of pi0 analysis according to the MC flag
 154
 155   this->CreateHistosEvent(listQA);
 156   this->CreateHistosCaloCellsQA(listQA);
 157   this->CreateHistosCaloCluster(listQA);
 158   this->CreateHistosPi0(listRD);
 159   this->CreateHistosMixPi0(listRD);
 160   if (fgIsMC) this->CreateHistosMC(listMC);
 161
 162   return;
 163 }
 164
 165 //_____________________________________________________________________________
 166 void AliPHOSpPbPi0Header::CreateHistosEvent(TList *listQA)
 167 {
 168   // create histograms at event level
 169
 170   const Int_t nhs    = 3;
 171   TString tName[nhs] = { "Centrality", "Vz", "Pileup" };
 172   Int_t   nbins[nhs] = {         110 , 500 ,     500  };
 173   Double_t xlow[nhs] = {         -10., -25.,     -25. };
 174   Double_t  xup[nhs] = {         100.,  25.,      25. };
 175
 176   TString hName;
 177   TH1D   *histo = 0;
 178   for (Int_t i=0; i<nhs; i++) {
 179     hName = Form("hEvent_%s", tName[i].Data());
 180     histo = new TH1D(hName.Data(), hName.Data(), nbins[i], xlow[i], xup[i]);
 181     histo->Sumw2(); listQA->Add(histo); histo = 0;
 182
 183     hName = Form("hEventSel_%s", tName[i].Data());
 184     histo = new TH1D(hName.Data(), hName.Data(), nbins[i], xlow[i], xup[i]);
 185     histo->Sumw2(); listQA->Add(histo); histo = 0;
 186   }
 187
 188   return;
 189 }
 190
 191 //_____________________________________________________________________________
 192 void AliPHOSpPbPi0Header::CreateHistosCaloCellsQA(TList *listQA)
 193 {
 194   // create QA histograms for CaloCells
 195
 196   TH1D *histo1       = 0;
 197   TH2D *histo2       = 0;
 198   const Int_t nhs    = 3;
 199   const Int_t nMod   = 3;
 200   TString hName = "hCaloCellsQA_E";
 201   histo1 = new TH1D(hName, hName, 300, 0., 30.);
 202   histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 203   for (Int_t i=0; i<nMod; i++) {
 204     hName  = Form("hCaloCellsQA_E_Mod%d", i+1);
 205     histo1 = new TH1D(hName.Data(), hName.Data(), 300, 0., 30.);
 206     histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 207   }
 208
 209   TString tName[nhs] = { "CentNCells", "Nxz", "Exz" };
 210   Int_t  nbinsx[nhs] = {          10 ,   64 ,  64   };
 211   Double_t xlow[nhs] = {           0.,   0.5,   0.5 };
 212   Double_t  xup[nhs] = {         100.,  64.5,  64.5 };
 213   Int_t  nbinsy[nhs] = {        1000 ,   56 ,  56   };
 214   Double_t ylow[nhs] = {           0.,   0.5,   0.5 };
 215   Double_t  yup[nhs] = {        1000.,  56.5,  56.5 };
 216
 217   for (Int_t i=0; i<nhs; i++) {
 218     if (i == 0) {
 219       hName  = Form("hCaloCellsQA_%s", tName[i].Data());
 220       histo2 = new TH2D(hName.Data(), hName.Data(), nbinsx[i], xlow[i], xup[i], nbinsy[i], ylow[i], yup[i]);
 221       histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
 222     } else {
 223       for (Int_t j=0; j<nMod; j++) {
 224         hName  = Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), j+1);
 225         histo2 = new TH2D(hName.Data(), hName.Data(), nbinsx[i], xlow[i], xup[i], nbinsy[i], ylow[i], yup[i]);
 226         histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
 227       }
 228     }
 229   }
 230
 231   return;
 232 }
 233
 234 //_____________________________________________________________________________
 235 void AliPHOSpPbPi0Header::CreateHistosCaloCluster(TList *listQA)
 236 {
 237   // create histograms for CaloCluster
 238
 239   const Int_t nMods = 3;
 240   const Int_t nTRUs = 8;
 241                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs };   // PID
 242   TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"           };
 243                          // { kPtClu, kEtaClu, kPhiClu, kM02Clu, kM20Clu, kTOFClu, kNCellsClu, kNClustersClu, kVarsClu };   // clusters
 244   TString tName[kVarsClu] = {   "Pt",   "Eta",   "Phi",   "M02",   "M20",   "TOF",   "NCells",   "NClusters"           };
 245   Int_t    bins[kVarsClu] = {   500 ,  300   ,   120  ,    200 ,    200 ,    400 ,       200 ,          20             };
 246   Double_t xMin[kVarsClu] = {     0.,   -0.15,     4.4,      0.,      0.,   -2e-6,         0.,           0.            };
 247   Double_t xMax[kVarsClu] = {    50.,    0.15,     5.6,      2.,      2.,    2e-6,      1000.,          20.            };
 248
 249   TH1I *hncells = 0;
 250   TH1D *histo1  = 0;
 251   TH2D *histo2  = 0;
 252   TString hName;
 253
 254   hName   = Form("hCaloCluster_%s", tName[kNCellsClu].Data());
 255   hncells = new TH1I(hName.Data(), hName.Data(), bins[kNCellsClu], (Int_t)xMin[kNCellsClu], (Int_t)xMax[kNCellsClu]);
 256   hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 257
 258   hName   = Form("hCaloCluster_%s", tName[kNClustersClu].Data());
 259   hncells = new TH1I(hName.Data(), hName.Data(), bins[kNClustersClu], (Int_t)xMin[kNClustersClu], (Int_t)xMax[kNClustersClu]);
 260   hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 261
 262   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 263     hName  = Form("hCaloCluster_%s_%s", tName[kPtClu].Data(), namePID[iPID].Data());
 264     histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
 265     histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 266
 267     for (Int_t icent=0; icent<fgNCent; icent++) {
 268       hName  = Form("hCaloCluster_%s_%s_cent%d", tName[kPtClu].Data(), namePID[iPID].Data(), icent);
 269       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
 270       histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 271     }
 272
 273     hName  = Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(), namePID[iPID].Data());
 274     histo2 = new TH2D(hName.Data(), hName.Data(), bins[kEtaClu], xMin[kEtaClu], xMax[kEtaClu], bins[kPhiClu], xMin[kPhiClu], xMax[kPhiClu]);
 275     histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
 276     for (Int_t i=0; i<kVarsClu-3; i++) {
 277       hName  = Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
 278       histo2 = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
 279       histo2->Sumw2(); listQA->Add(histo2); histo2 = 0;
 280     }
 281   }
 282
 283   for (Int_t iMod=0; iMod<nMods; iMod++) {
 284     hName   = Form("hCaloCluster_%s_Mod%d", tName[kNCellsClu].Data(), iMod+1);
 285     hncells = new TH1I(hName.Data(), hName.Data(), bins[kNCellsClu], (Int_t)xMin[kNCellsClu], (Int_t)xMax[kNCellsClu]);
 286     hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 287
 288     hName   = Form("hCaloCluster_%s_Mod%d", tName[kNClustersClu].Data(), iMod+1);
 289     hncells = new TH1I(hName.Data(), hName.Data(), bins[kNClustersClu], (Int_t)xMin[kNClustersClu], (Int_t)xMax[kNClustersClu]);
 290     hncells->Sumw2(); listQA->Add(hncells); hncells = 0;
 291
 292     for (Int_t iTRU=0; iTRU<nTRUs; iTRU++) {
 293       hName  = Form("hCaloCluster_%s_Mod%d_TRU%d", tName[kPtClu].Data(), iMod+1, iTRU+1);
 294       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
 295       histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 296     }
 297
 298     for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 299       hName  = Form("hCaloCluster_%s_Mod%d_%s", tName[kPtClu].Data(), iMod+1, namePID[iPID].Data());
 300       histo1 = new TH1D(hName.Data(), hName.Data(), bins[kPtClu], xMin[kPtClu], xMax[kPtClu]);
 301       histo1->Sumw2(); listQA->Add(histo1); histo1 = 0;
 302     }
 303   }
 304
 305   return;
 306 }
 307
 308 //_____________________________________________________________________________
 309 void AliPHOSpPbPi0Header::CreateHistosPi0(TList *listRD)
 310 {
 311   // create histograms for Pi0
 312
 313   const Int_t nComb = 2;
 314   TString srcMod[nComb]  =  { "Mod11", "Mod33" };
 315                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs      };  // PID
 316   TString namePID[kPIDs] =  { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"                };
 317                          // { kPtPi0, kEtaPi0, kPhiPi0, kAsyPi0,   kAnglePi0, kInvMassPi0, kVarsPi0 };  // pi0
 318   TString tName[kVarsPi0] = {   "Pt",   "Eta",   "Phi",   "Asy",     "Angle",   "InvMass"           };
 319   Int_t    bins[kVarsPi0] = {   500 ,  300   ,   120  ,    100 ,        180 ,       1000            };
 320   Double_t xMin[kVarsPi0] = {     0.,   -0.15,     4.4,      0.,          0.,          0.           };
 321   Double_t xMax[kVarsPi0] = {    50.,    0.15,     5.6,      1., TMath::Pi(),          1.           };
 322
 323   const Int_t    nPtbins     = 32;
 324   Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
 325                                   2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
 326                                   6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
 327                                  30.0, 36.0, 44.0 };
 328
 329   TH2D   *histo = 0;
 330   TString hName;
 331
 332   hName = Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data());
 333   histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaPi0], xMin[kEtaPi0], xMax[kEtaPi0], bins[kPhiPi0], xMin[kPhiPi0], xMax[kPhiPi0]);
 334   histo->Sumw2(); listRD->Add(histo); histo = 0;
 335   for (Int_t i=0; i<kVarsPi0-2; i++) {
 336     hName = Form("hPi0_%s%s", tName[0].Data(), tName[i+1].Data());
 337     histo = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
 338     histo->Sumw2(); listRD->Add(histo); histo = 0;
 339   }
 340
 341   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 342     hName = Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data());
 343     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
 344     histo->Sumw2(); listRD->Add(histo); histo = 0;
 345
 346     hName = Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(), namePID[iPID].Data());
 347     histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaPi0], xMin[kEtaPi0], xMax[kEtaPi0], bins[kPhiPi0], xMin[kPhiPi0], xMax[kPhiPi0]);
 348     histo->Sumw2(); listRD->Add(histo); histo = 0;
 349     for (Int_t icent=0; icent<fgNCent; icent++) {
 350       hName = Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), icent);
 351       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
 352       histo->Sumw2(); listRD->Add(histo); histo = 0;
 353     }
 354
 355     for (Int_t i=0; i<kVarsPi0-2; i++) {
 356       hName = Form("hPi0_%s%s_%s", tName[0].Data(), tName[i+1].Data(), namePID[iPID].Data());
 357       histo = new TH2D(hName.Data(), hName.Data(), bins[0], xMin[0], xMax[0], bins[i+1], xMin[i+1], xMax[i+1]);
 358       histo->Sumw2(); listRD->Add(histo); histo = 0;
 359     }
 360     for (Int_t iComb=0; iComb<nComb; iComb++) {
 361       hName = Form("hPi0_%s%s_%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
 362       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassPi0], xMin[kInvMassPi0], xMax[kInvMassPi0]);
 363       histo->Sumw2(); listRD->Add(histo); histo = 0;
 364     }
 365   }
 366
 367   return;
 368 }
 369
 370 //_____________________________________________________________________________
 371 void AliPHOSpPbPi0Header::CreateHistosMixPi0(TList *listRD)
 372 {
 373   // create Histograms for Mixed Pi0
 374
 375   const Int_t nComb = 2;
 376   TString  srcMod[nComb]     = { "Mod11", "Mod33" };
 377                              //{ kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs    };  // PID
 378   TString namePID[kPIDs]     = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"              };
 379                              //{ kPtMixPi0, kEtaMixPi0, kPhiMixPi0, kInvMassMixPi0, kVarsMixPi0      };  // Mix Pi0
 380   TString tName[kVarsMixPi0] = {      "Pt",      "Eta",      "Phi",      "InvMass"                   };
 381   Int_t    bins[kVarsMixPi0] = {      500 ,     300   ,      120  ,          1000                    };
 382   Double_t xMin[kVarsMixPi0] = {        0.,      -0.15,        4.4,             0.                   };
 383   Double_t xMax[kVarsMixPi0] = {       50.,       0.15,        5.6,             1.                   };
 384
 385   const Int_t    nPtbins     = 32;
 386   Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
 387                                   2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
 388                                   6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
 389                                  30.0, 36.0, 44.0 };
 390
 391   TH2D   *histo = 0;
 392   TString hName;
 393
 394   hName = Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data());
 395   histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaMixPi0], xMin[kEtaMixPi0], xMax[kEtaMixPi0],
 396                                                bins[kPhiMixPi0], xMin[kPhiMixPi0], xMax[kPhiMixPi0]);
 397   histo->Sumw2(); listRD->Add(histo); histo = 0;
 398
 399   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 400     hName = Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data());
 401     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
 402     histo->Sumw2(); listRD->Add(histo); histo = 0;
 403     hName = Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(), namePID[iPID].Data());
 404     histo = new TH2D(hName.Data(), hName.Data(), bins[kEtaMixPi0], xMin[kEtaMixPi0], xMax[kEtaMixPi0],
 405                                                  bins[kPhiMixPi0], xMin[kPhiMixPi0], xMax[kPhiMixPi0]);
 406     histo->Sumw2(); listRD->Add(histo); histo = 0;
 407     for (Int_t icent=0; icent<fgNCent; icent++) {
 408       hName = Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), icent);
 409       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
 410       histo->Sumw2(); listRD->Add(histo); histo = 0;
 411     }
 412     for (Int_t iComb=0; iComb<nComb; iComb++) {
 413       hName = Form("hMixPi0_%s%s_%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(), namePID[iPID].Data(), srcMod[iComb].Data());
 414       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMixPi0], xMin[kInvMassMixPi0], xMax[kInvMassMixPi0]);
 415       histo->Sumw2(); listRD->Add(histo); histo = 0;
 416     }
 417   }
 418
 419   return;
 420 }
 421
 422 //_____________________________________________________________________________
 423 void AliPHOSpPbPi0Header::CreateHistosMC(TList *listMC)
 424 {
 425   // create Histograms for MC
 426
 427   const Int_t acc    = 2;
 428   const char *aName[acc]   = {"Pi0InAcc", "GammaInAcc"};
 429   const Int_t types  = 5;
 430   const char *pName[types] = {"InclusivePi0", "PrimaryPi0", "2ndPi0FromK0s", "2ndPi0FromMaterials", "2ndPi0FromOtherDecays"};
 431                          //  { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs    };  // PID
 432   TString namePID[kPIDs]   = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"              };
 433                          //  { kPtMC, kRapidityMC, kRadiusMC,         kPhiMC, kInvMassMC, kVarsMC  };  // MC
 434   TString tName[kVarsMC]   = {  "Pt",  "Rapidity",  "Radius",          "Phi",  "InvMass"           };
 435   Int_t    bins[kVarsMC]   = {  500 ,        600 ,      1000,           360 ,       500            };
 436   Double_t xMin[kVarsMC]   = {    0.,         -3.,        0.,             0.,         0.           };
 437   Double_t xMax[kVarsMC]   = {   50.,          3.,      500., TMath::TwoPi(),         0.5          };
 438
 439   const Int_t    nPtbins     = 32;
 440   Double_t ptbins[nPtbins+1] = {  0.8,  1.0,  1.2,  1.4,  1.6,  1.8,  2.0,  2.2,  2.4,  2.6,
 441                                   2.8,  3.0,  3.2,  3.4,  3.6,  3.8,  4.0,  4.5,  5.0,  5.5,
 442                                   6.0,  7.0,  8.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0,
 443                                  30.0, 36.0, 44.0 };
 444
 445   TH2D   *histo = 0;
 446   TString hName;
 447
 448   for (Int_t iType=0; iType<types; iType++) {
 449     hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kRapidityMC].Data());
 450     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRapidityMC], xMin[kRapidityMC], xMax[kRapidityMC]);
 451     histo->Sumw2(); listMC->Add(histo); histo = 0;
 452     hName  = Form("hMC%s_%s%s_Check", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data());
 453     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 454     histo->Sumw2(); listMC->Add(histo); histo = 0;
 455     hName  = Form("hMC%s_%s%s", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data());
 456     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 457     histo->Sumw2(); listMC->Add(histo); histo = 0;
 458     hName  = Form("hMC%s_%s%s_MergedCluster", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data());
 459     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 460     histo->Sumw2(); listMC->Add(histo); histo = 0;
 461     hName  = Form("hMC%s_%s%s_Reco", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data());
 462     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 463     histo->Sumw2(); listMC->Add(histo); histo = 0;
 464     if (strncmp(pName[iType], "2ndPi0FromK0s", 13) == 0) {
 465       hName  = Form("hMC%s_%sMatrix", pName[iType], tName[kPtMC].Data());
 466       histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
 467       histo->Sumw2(); listMC->Add(histo); histo = 0;
 468     }
 469     for (Int_t iAcc=0; iAcc<acc; iAcc++) {
 470       hName  = Form("hMC%s_%s%s_%s", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), aName[iAcc]);
 471       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 472       histo->Sumw2(); listMC->Add(histo); histo = 0;
 473       if (strncmp(pName[iType], "2ndPi0FromK0s", 13) == 0) {
 474         hName  = Form("hMC%s_%sMatrix_%s", pName[iType], tName[kPtMC].Data(), aName[iAcc]);
 475         histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
 476         histo->Sumw2(); listMC->Add(histo); histo = 0;
 477       }
 478     }
 479     for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 480       hName  = Form("hMC%s_%s%s_%s_Reco", pName[iType], tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data());
 481       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 482       histo->Sumw2(); listMC->Add(histo); histo = 0;
 483     }
 484
 485     for (Int_t iCent=0; iCent<fgNCent; iCent++) {
 486       hName  = Form("hMC%s_%s%s_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRapidityMC].Data(), iCent);
 487       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRapidityMC], xMin[kRapidityMC], xMax[kRapidityMC]);
 488       histo->Sumw2(); listMC->Add(histo); histo = 0;
 489       hName  = Form("hMC%s_%s%s_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), iCent);
 490       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 491       histo->Sumw2(); listMC->Add(histo); histo = 0;
 492       hName  = Form("hMC%s_%s%s_MergedCluster_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), iCent);
 493       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 494       histo->Sumw2(); listMC->Add(histo); histo = 0;
 495       hName  = Form("hMC%s_%s%s_Reco_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), iCent);
 496       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 497       histo->Sumw2(); listMC->Add(histo); histo = 0;
 498       if (strncmp(pName[iType], "2ndPi0FromK0s", 13) == 0) {
 499         hName  = Form("hMC%s_%sMatrix_cent%d", pName[iType], tName[kPtMC].Data(), iCent);
 500         histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
 501         histo->Sumw2(); listMC->Add(histo); histo = 0;
 502       }
 503       for (Int_t iAcc=0; iAcc<acc; iAcc++) {
 504         hName  = Form("hMC%s_%s%s_%s_cent%d", pName[iType], tName[kPtMC].Data(), tName[kRadiusMC].Data(), aName[iAcc], iCent);
 505         histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kRadiusMC], xMin[kRadiusMC], xMax[kRadiusMC]);
 506         histo->Sumw2(); listMC->Add(histo); histo = 0;
 507         if (strncmp(pName[iType], "2ndPi0FromK0s", 13) == 0) {
 508           hName  = Form("hMC%s_%sMatrix_%s_cent%d", pName[iType], tName[kPtMC].Data(), aName[iAcc], iCent);
 509           histo = new TH2D(hName.Data(), hName.Data(), bins[kPtMC], xMin[kPtMC], xMax[kPtMC], bins[kPtMC], xMin[kPtMC], xMax[kPtMC]);
 510           histo->Sumw2(); listMC->Add(histo); histo = 0;
 511         }
 512       }
 513       for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 514         hName  = Form("hMC%s_%s%s_%s_Reco_cent%d", pName[iType], tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data(), iCent);
 515         histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 516         histo->Sumw2(); listMC->Add(histo); histo = 0;
 517       }
 518     }
 519   }
 520   for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 521     hName  = Form("hMCInclusiveCvtedPi0_%s%s_%s_Reco", tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data());
 522     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 523     histo->Sumw2(); listMC->Add(histo); histo = 0;
 524     hName  = Form("hMCInclusivePurePi0_%s%s_%s_Reco", tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data());
 525     histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 526     histo->Sumw2(); listMC->Add(histo); histo = 0;
 527     for (Int_t iCent=0; iCent<fgNCent; iCent++) {
 528       hName  = Form("hMCInclusiveCvtedPi0_%s%s_%s_Reco_cent%d", tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data(), iCent);
 529       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 530       histo->Sumw2(); listMC->Add(histo); histo = 0;
 531       hName  = Form("hMCInclusivePurePi0_%s%s_%s_Reco_cent%d", tName[kPtMC].Data(), tName[kInvMassMC].Data(), namePID[iPID].Data(), iCent);
 532       histo = new TH2D(hName.Data(), hName.Data(), nPtbins, ptbins, bins[kInvMassMC], xMin[kInvMassMC], xMax[kInvMassMC]);
 533       histo->Sumw2(); listMC->Add(histo); histo = 0;
 534     }
 535   }
 536
 537   return;
 538 }
 539
 540 //_____________________________________________________________________________
 541 void AliPHOSpPbPi0Header::FillHistosEvent(TList *listQA)
 542 {
 543   // fill histograms at event level according to event selection cuts
 544
 545   if (!listQA) return;
 546
 547   const Int_t nhs      = 2;
 548   TString tName[nhs+1] = { "Centrality",    "Vz", "Pileup" };
 549   Double_t dist[nhs]   = {  fCentrality, fVtx[2]           };
 550   for (Int_t i=nhs; i--;) {
 551     ((TH1D*)listQA->FindObject(Form("hEvent_%s", tName[i].Data())))->Fill(dist[i]);
 552     if (this->IsSelected()) ((TH1D*)listQA->FindObject(Form("hEventSel_%s", tName[i].Data())))->Fill(dist[i]);
 553   }
 554   if (fIsPileup) {
 555     ((TH1D*)listQA->FindObject(Form("hEvent_%s", tName[nhs].Data())))->Fill(dist[nhs-1]);
 556     if (this->IsSelected())  ((TH1D*)listQA->FindObject(Form("hEventSel_%s", tName[nhs].Data())))->Fill(dist[nhs-1]);
 557   }
 558
 559   return;
 560 }
 561
 562 //_____________________________________________________________________________
 563 void AliPHOSpPbPi0Header::FillHistosCaloCellsQA(TList *listQA, AliVCaloCells* const cells, AliPHOSGeoUtils* const phosGeo)
 564 {
 565   // fill QA histograms for calocells according to evnet selection cuts
 566
 567   if (!listQA) return;
 568
 569   const Int_t nhs    = 3;
 570   TString tName[nhs] = { "CentNCells", "Nxz", "Exz" };
 571
 572   ((TH2D*)listQA->FindObject(Form("hCaloCellsQA_%s", tName[0].Data())))->Fill(fCentrality, cells->GetNumberOfCells());
 573   for (Int_t iCell=0; iCell<cells->GetNumberOfCells(); iCell++) {
 574     Int_t relID[4] = {0,0,0,0};
 575     phosGeo->AbsToRelNumbering(cells->GetCellNumber(iCell), relID); // Int_t mod   = relID[0]; Int_t cellX = relID[2]; Int_t cellZ = relID[3];
 576     Double_t wight[2] = { 1., cells->GetAmplitude(iCell)};
 577     ((TH1D*)listQA->FindObject(Form("hCaloCellsQA_%s", "E")))->Fill(wight[1]);
 578     ((TH1D*)listQA->FindObject(Form("hCaloCellsQA_%s_Mod%d", "E", relID[0])))->Fill(wight[1]);
 579     for (Int_t i=1; i<nhs; i++)
 580       ((TH2D*)listQA->FindObject(Form("hCaloCellsQA_%s_Mod%d", tName[i].Data(), relID[0])))->Fill(relID[2], relID[3], wight[i-1]);
 581   }
 582
 583   return;
 584 }
 585
 586 //_____________________________________________________________________________
 587 void AliPHOSpPbPi0Header::FillHistosCaloCluster(TList *listQA, TClonesArray* const caloClArr, Int_t cent)
 588 {
 589   // fill histograms for calocluster according to evnet && calocluster candidates selection cuts
 590
 591   if (!listQA)  return;
 592                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs       };  // PID
 593   TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"                 };
 594   UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)     };
 595                          // { kPtClu, kEtaClu, kPhiClu, kM02Clu, kM20Clu, kTOFClu, kNCellsClu, kNClustersClu, kVarsClu };   // clusters
 596   TString tName[kVarsClu] = {   "Pt",   "Eta",   "Phi",   "M02",   "M20",   "TOF",   "NCells",   "NClusters"           };
 597
 598   Int_t entries = caloClArr->GetEntriesFast();
 599   Int_t iMod = 0, iTRU = 0, nclusters[3] = { 0, 0, 0 };
 600   UInt_t pidBit = 0;
 601   AliCaloClusterInfo *caloCluster = 0;
 602   TLorentzVector gamma;
 603   for(Int_t i=0; i<entries; i++) { // Loop calo cluster
 604     caloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
 605     iMod   = caloCluster->GetModule();   nclusters[iMod-1]++;
 606     iTRU   = caloCluster->GetTRUNumber();
 607     pidBit = caloCluster->GetPIDBit();
 608     gamma  = caloCluster->LorentzVector();
 609     Double_t vars[kVarsClu-1] = { gamma.E(),
 610                                   gamma.Eta(),
 611                                   (gamma.Phi()+TMath::TwoPi()),
 612                                   caloCluster->GetM02(),
 613                                   caloCluster->GetM20(),
 614                                   caloCluster->GetTOF(),
 615                                   (Double_t)caloCluster->GetNCells() };
 616
 617     ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s", tName[kNCellsClu].Data())))->Fill((Int_t)vars[kNCellsClu]);
 618     ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNCellsClu].Data(), iMod)))->Fill((Int_t)vars[kNCellsClu]);
 619     ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d_TRU%d", tName[kPtClu].Data(), iMod, iTRU)))->Fill(vars[kPtClu]);
 620
 621     for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 622       if ((pidBit&PIDBIT[iPID])==PIDBIT[iPID]) {
 623         ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_%s", tName[kPtClu].Data(), namePID[iPID].Data())))->Fill(vars[kPtClu]);
 624         ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d_%s", tName[kPtClu].Data(), iMod, namePID[iPID].Data())))->Fill(vars[kPtClu]);
 625         ((TH1D*)listQA->FindObject(Form("hCaloCluster_%s_%s_cent%d", tName[kPtClu].Data(), namePID[iPID].Data(), cent)))->Fill(vars[kPtClu]);
 626         ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s_%s", tName[kEtaClu].Data(), tName[kPhiClu].Data(),
 627                                                               namePID[iPID].Data())))->Fill(vars[kEtaClu], vars[kPhiClu]);
 628         for (Int_t j=0; j<kVarsClu-3; j++)
 629           ((TH2D*)listQA->FindObject(Form("hCaloCluster_%s%s_%s", tName[0].Data(), tName[j+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[j+1]);
 630       }
 631     }
 632   } // End loop calo cluster
 633
 634   ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s", tName[kNClustersClu].Data())))->Fill(entries);
 635   for (Int_t jMod=1; jMod<4; jMod++) {
 636     if (jMod==2) continue;
 637     ((TH1I*)listQA->FindObject(Form("hCaloCluster_%s_Mod%d", tName[kNClustersClu].Data(), jMod)))->Fill(nclusters[jMod-1]);
 638   }
 639
 640   return;
 641 }
 642
 643 //_____________________________________________________________________________
 644 void AliPHOSpPbPi0Header::FillHistosPi0(TList *listRD, TClonesArray* const caloClArr, Int_t cent)
 645 {
 646   // fill histograms for pi0 according to evnet && pi0 candidates selection cuts
 647
 648   if (!listRD) return;
 649                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
 650   TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
 651   UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
 652                          // { kPtPi0, kEtaPi0, kPhiPi0, kAsyPi0, kAnglePi0, kInvMassPi0, kVarsPi0  };  // Pi0
 653   TString tName[kVarsPi0] = {   "Pt",   "Eta",   "Phi",   "Asy",   "Angle",   "InvMass"            };
 654
 655   Int_t entries = caloClArr->GetEntriesFast();
 656   Int_t iMod = 0, jMod = 0, srcMod = 0;
 657   UInt_t iPIDBit = 0, jPIDBit = 0;
 658   Double_t vars[kVarsPi0];
 659   AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
 660   TLorentzVector iGamma, jGamma;
 661   for(Int_t i=0; i<entries-1; i++) { // Loop calo cluster i
 662     iCaloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
 663     iMod    = iCaloCluster->GetModule();
 664     iPIDBit = iCaloCluster->GetPIDBit();
 665     iGamma  = iCaloCluster->LorentzVector();
 666
 667     for (Int_t j=i+1; j<entries; j++) { // Loop calo cluster j
 668       jCaloCluster = (AliCaloClusterInfo*)caloClArr->At(j);
 669       jMod    = jCaloCluster->GetModule();
 670       jPIDBit = jCaloCluster->GetPIDBit();
 671       jGamma  = jCaloCluster->LorentzVector();
 672       if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
 673
 674       if (iMod==1 && jMod==1)       srcMod = 11; // Mod11
 675       else if (iMod==3 && jMod==3)  srcMod = 33; // Mod33
 676       else if (iMod==2 && jMod==2)  srcMod = 22; // Mod22
 677       else if (iMod==1 && jMod==2)  srcMod = 12; // Mod12
 678       else if (iMod==2 && jMod==3)  srcMod = 23; // Mod23
 679
 680       vars[kPtPi0]      = (iGamma+jGamma).E();
 681       vars[kEtaPi0]     = (iGamma+jGamma).Eta();
 682       vars[kPhiPi0]     = (iGamma+jGamma).Phi() + TMath::TwoPi();
 683       vars[kInvMassPi0] = (iGamma+jGamma).M();
 684       vars[kAsyPi0]     = TMath::Abs((iGamma-jGamma).E())/(iGamma+jGamma).E();
 685       vars[kAnglePi0]   = jGamma.Angle(iGamma.Vect());
 686
 687       ((TH2D*)listRD->FindObject(Form("hPi0_%s%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data())))->Fill(vars[kEtaPi0], vars[kPhiPi0]);
 688       for (Int_t k=0; k<kVarsPi0-2; k++) {
 689         ((TH2D*)listRD->FindObject(Form("hPi0_%s%s", tName[0].Data(), tName[k+1].Data())))->Fill(vars[0], vars[k+1]);
 690       }
 691
 692       for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 693         if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
 694           ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[kEtaPi0].Data(), tName[kPhiPi0].Data(),
 695                                                         namePID[iPID].Data())))->Fill(vars[kEtaPi0], vars[kPhiPi0]);
 696           ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
 697                                                         namePID[iPID].Data())))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
 698           ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s_Mod%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
 699                                                               namePID[iPID].Data(), srcMod)))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
 700           ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s_cent%d", tName[kPtPi0].Data(), tName[kInvMassPi0].Data(),
 701                                                                namePID[iPID].Data(), cent)))->Fill(vars[kPtPi0], vars[kInvMassPi0]);
 702           for (Int_t k=0; k<kVarsPi0-2; k++) {
 703             ((TH2D*)listRD->FindObject(Form("hPi0_%s%s_%s", tName[0].Data(), tName[k+1].Data(), namePID[iPID].Data())))->Fill(vars[0], vars[k+1]);
 704           }
 705         }
 706       }
 707       jCaloCluster = 0;
 708     }// End loop j
 709     iCaloCluster = 0;
 710   } // End loop i
 711
 712   return;
 713 }
 714
 715 //_____________________________________________________________________________
 716 void AliPHOSpPbPi0Header::FillHistosMixPi0(TList *listRD, TClonesArray* const iCaloClArr, TList* const eventList, Int_t cent)
 717 {
 718   // fill histograms for mixed pi0 according to evnet && pi0 candidates selection cuts
 719
 720   if (!listRD) return;
 721                         // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
 722   TString namePID[kPIDs] = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
 723   UInt_t  PIDBIT[kPIDs]  = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
 724                             // { kPtMixPi0, kEtaMixPi0, kPhiMixPi0, kInvMassMixPi0, kVarsMixPi0   };   // Mix pi0
 725   TString tName[kVarsMixPi0] = {      "Pt",      "Eta",      "Phi",      "InvMass"                };
 726
 727   Int_t iEntries = iCaloClArr->GetEntriesFast();
 728   Int_t iMod = 0, jMod = 0, srcMod = 0;
 729   UInt_t iPIDBit = 0, jPIDBit = 0;
 730   Double_t vars[kVarsMixPi0];
 731   AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
 732   TLorentzVector iGamma, jGamma;
 733   TClonesArray *jCaloClArr = 0;
 734   for (Int_t i=0; i<iEntries; i++) { // Loop calo cluster i
 735     iCaloCluster=(AliCaloClusterInfo*)iCaloClArr->At(i);
 736     iMod    = iCaloCluster->GetModule();
 737     iPIDBit = iCaloCluster->GetPIDBit();
 738     iGamma  = iCaloCluster->LorentzVector();
 739
 740     Int_t nev = eventList->GetSize();
 741     for(Int_t ev=0; ev<nev; ev++){ // Loop events for mixing
 742       jCaloClArr     = static_cast<TClonesArray*>(eventList->At(ev));
 743       Int_t jEntries = jCaloClArr->GetEntriesFast();
 744       for(Int_t j=0; j<jEntries; j++){ // Loop calo cluster j
 745         jCaloCluster=(AliCaloClusterInfo*)jCaloClArr->At(j);
 746         jMod    = jCaloCluster->GetModule();   if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
 747         jPIDBit = jCaloCluster->GetPIDBit();
 748         jGamma  = jCaloCluster->LorentzVector();
 749
 750         if (iMod==1 && jMod==1 )                                srcMod = 11; // Mod11
 751         else if (iMod==3 && jMod==3)                            srcMod = 33; // Mod33
 752         else if (iMod==2 && jMod==2)                            srcMod = 22; // Mod22
 753         else if ((iMod==1 && jMod==2) || (iMod==2 && jMod==1))  srcMod = 12; // Mod12
 754         else if ((iMod==2 && jMod==3) || (iMod==3 && jMod==2))  srcMod = 23; // Mod23
 755
 756         vars[kPtMixPi0]      = (iGamma+jGamma).E();
 757         vars[kEtaMixPi0]     = (iGamma+jGamma).Eta();
 758         vars[kPhiMixPi0]     = (iGamma+jGamma).Phi() + TMath::TwoPi();
 759         vars[kInvMassMixPi0] = (iGamma+jGamma).M();
 760
 761         ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
 762         for (Int_t iPID=0; iPID<kPIDs; iPID++) {
 763           if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
 764             ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s", tName[kEtaMixPi0].Data(), tName[kPhiMixPi0].Data(),
 765                                                                namePID[iPID].Data())))->Fill(vars[kEtaMixPi0], vars[kPhiMixPi0]);
 766             ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
 767                                                                namePID[iPID].Data())))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
 768             ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s_Mod%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
 769                                                                      namePID[iPID].Data(), srcMod)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
 770             ((TH2D*)listRD->FindObject(Form("hMixPi0_%s%s_%s_cent%d", tName[kPtMixPi0].Data(), tName[kInvMassMixPi0].Data(),
 771                                                                       namePID[iPID].Data(), cent)))->Fill(vars[kPtMixPi0], vars[kInvMassMixPi0]);
 772           }
 773         }
 774         jCaloCluster = 0;
 775       } // End loop j
 776       jCaloClArr = 0;
 777     } // End loop events
 778     iCaloCluster = 0;
 779   } // End loop i
 780
 781   return;
 782 }
 783
 784 //_____________________________________________________________________________
 785 void AliPHOSpPbPi0Header::FillHistosMC(TList *listMC, AliStack* const stack, TClonesArray* const caloClArr, AliPHOSGeoUtils* const phosGeo, Int_t cent)
 786 {
 787   // fill histograms for ESD MC
 788   if (!listMC)  return;
 789   if (!stack)   return;
 790
 791                          // { kAll,  kCpv,  kDisp,  kBoth,  kCpv2,  kDisp2,  kBoth2,     kPIDs     };  // PID
 792   TString namePID[kPIDs]  = { "All", "Cpv", "Disp", "Both", "Cpv2", "Disp2", "Both2"               };
 793   UInt_t  PIDBIT[kPIDs]   = { 0x0,  BIT(0), BIT(1), BIT(0)|BIT(1), BIT(2), BIT(3), BIT(2)|BIT(3)   };
 794
 795   // MC truth info
 796   TString pName;
 797   const Int_t nParticles = stack->GetNtrack();
 798   Int_t iMod = 0, jMod = 0;
 799   for(Int_t iParticle=0; iParticle<nParticles; iParticle++) {
 800     TParticle* pMC  = stack->Particle(iParticle);
 801     if (pMC->GetPdgCode() != 111)  continue;   // Pi0
 802     if (pMC->GetNDaughters() != 2) continue;   // Do not account Dalitz decays
 803     TParticle *iGamma = stack->Particle(pMC->GetFirstDaughter());
 804     TParticle *jGamma = stack->Particle(pMC->GetLastDaughter());
 805     if (!(iGamma->GetPdgCode()==22 && jGamma->GetPdgCode()==22)) continue; // Get into the branch: pi0->gamma+gamma
 806 //  if (!(iGamma->GetPdgCode()==22 && jGamma->GetPdgCode()==22)) { Printf("PDG: P1=%d, P2=%d\n", iGamma->GetPdgCode(), jGamma->GetPdgCode()); continue; }
 807
 808     pName = ClassifyMCPi0(iParticle, stack);   // Classify pi0
 809
 810     //Pi0 Information
 811     Double_t pt       = pMC->Pt();
 812     Double_t radius   = pMC->R();
 813     Double_t rapidity = pMC->Y();
 814     Double_t eta      = pMC->Eta();
 815     Double_t phi      = pMC->Phi(); while (phi<0.) phi += TMath::TwoPi(); while (phi>TMath::TwoPi()) phi -= TMath::TwoPi();
 816
 817     ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRapidity"))->Fill(pt, rapidity);
 818     ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRapidity_cent%d", cent)))->Fill(pt, rapidity);
 819
 820     ((TH2D*)listMC->FindObject(Form("hMC%s_PtRapidity", pName.Data())))->Fill(pt, rapidity);
 821     ((TH2D*)listMC->FindObject(Form("hMC%s_PtRapidity_cent%d", pName.Data(), cent)))->Fill(pt, rapidity);
 822
 823     ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius_Check"))->Fill(pt, radius);
 824     ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_Check", pName.Data())))->Fill(pt, radius);
 825     if (TMath::Abs(rapidity)<0.5) {
 826       ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius"))->Fill(pt, radius);
 827       ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_cent%d", cent)))->Fill(pt, radius);
 828
 829       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius", pName.Data())))->Fill(pt, radius);
 830       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_cent%d", pName.Data(), cent)))->Fill(pt, radius);
 831
 832       if (pName.Contains("K0s") && stack->IsPhysicalPrimary(pMC->GetFirstMother()) ) {
 833         TParticle* pMother = stack->Particle(pMC->GetFirstMother());
 834         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix", pName.Data())))->Fill(pt, pMother->Pt());
 835         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix_cent%d", pName.Data(), cent)))->Fill(pt, pMother->Pt());
 836       }
 837     }
 838
 839     if ((((phi>260.*TMath::DegToRad()) && (phi<280.*TMath::DegToRad())) || ((phi>300.*TMath::DegToRad())
 840            && (phi<320.*TMath::DegToRad()))) && (TMath::Abs(eta)<0.135)) {
 841       ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius_Pi0InAcc"))->Fill(pt, radius);
 842       ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_Pi0InAcc_cent%d", cent)))->Fill(pt, radius);
 843
 844       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_Pi0InAcc", pName.Data())))->Fill(pt, radius);
 845       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_Pi0InAcc_cent%d", pName.Data(), cent)))->Fill(pt, radius);
 846
 847       if (pName.Contains("K0s") && stack->IsPhysicalPrimary(pMC->GetFirstMother())) {
 848         TParticle* pMother = stack->Particle(pMC->GetFirstMother());
 849         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix_Pi0InAcc", pName.Data())))->Fill(pt, pMother->Pt());
 850         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix_Pi0InAcc_cent%d", pName.Data(), cent)))->Fill(pt, pMother->Pt());
 851       }
 852     }
 853
 854     //Pi0's daughter particles in PHOS acceptance
 855     iMod = HitPHOSModule(iGamma, phosGeo);
 856     jMod = HitPHOSModule(jGamma, phosGeo);
 857     if (iMod!=-1 && iMod!=2 && jMod!=-1 && jMod!=2 && TMath::Abs(iMod-jMod)<2) { // !remove module 2
 858       ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius_GammaInAcc"))->Fill(pt, radius);
 859       ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_GammaInAcc_cent%d", cent)))->Fill(pt, radius);
 860
 861       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_GammaInAcc", pName.Data())))->Fill(pt, radius);
 862       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_GammaInAcc_cent%d", pName.Data(), cent)))->Fill(pt, radius);
 863
 864       if (pName.Contains("K0s") && stack->IsPhysicalPrimary(pMC->GetFirstMother())) {
 865         TParticle* pMother = stack->Particle(pMC->GetFirstMother());
 866         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix_GammaInAcc", pName.Data())))->Fill(pt, pMother->Pt());
 867         ((TH2D*)listMC->FindObject(Form("hMC%s_PtMatrix_GammaInAcc_cent%d", pName.Data(), cent)))->Fill(pt, pMother->Pt());
 868       }
 869     }
 870   }
 871
 872   // Reconstruction info
 873   Bool_t isConverted_i = kFALSE, isConverted_j = kFALSE;
 874   Int_t  entries  = caloClArr->GetEntriesFast();
 875   Int_t  iPi0Indx =-1,  jPi0Indx   =-1;
 876   UInt_t iPIDBit  = 0,  jPIDBit    = 0;
 877   Double_t pi0Pt  = 0., pi0InvMass = 0., radius = 0.;
 878   AliCaloClusterInfo *iCaloCluster = 0, *jCaloCluster = 0;
 879   TLorentzVector iLorentz, jLorentz;
 880   for(Int_t i=0; i<entries; i++) { // Loop calo cluster i
 881     iCaloCluster = (AliCaloClusterInfo*)caloClArr->At(i);
 882     if (iCaloCluster->IsMergedClusterFromPi0(stack, iPi0Indx)) {
 883       pName  = ClassifyMCPi0(iPi0Indx, stack);
 884       pi0Pt  = (iCaloCluster->LorentzVector()).E();
 885       radius = stack->Particle(iPi0Indx)->R();
 886       ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius_MergedCluster"))->Fill(pi0Pt, radius);
 887       ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_MergedCluster_cent%d", cent)))->Fill(pi0Pt, radius);
 888
 889       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_MergedCluster", pName.Data())))->Fill(pi0Pt, radius);
 890       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_MergedCluster_cent%d", pName.Data(), cent)))->Fill(pi0Pt, radius);
 891     }
 892     else if (!iCaloCluster->IsClusterFromCvtedPi0(stack, isConverted_i, iPi0Indx) || i==entries-1) { iCaloCluster = 0; continue; }
 893
 894     for (Int_t j=i+1; j<entries; j++) { // Loop calo cluster j
 895       jCaloCluster = (AliCaloClusterInfo*)caloClArr->At(j);
 896
 897       if (!jCaloCluster->IsClusterFromCvtedPi0(stack, isConverted_j, jPi0Indx)) { jCaloCluster = 0; continue; }
 898       if (jPi0Indx != iPi0Indx) { jCaloCluster = 0; continue; }  // coming from the same pi0
 899       iMod     = iCaloCluster->GetModule();
 900       iPIDBit  = iCaloCluster->GetPIDBit();
 901       iLorentz = iCaloCluster->LorentzVector();
 902
 903       pName    = ClassifyMCPi0(iPi0Indx, stack);
 904       jMod     = jCaloCluster->GetModule();
 905       jPIDBit  = jCaloCluster->GetPIDBit();
 906       jLorentz = jCaloCluster->LorentzVector();
 907       if (TMath::Abs(iMod-jMod)>1) { jCaloCluster = 0; continue; }
 908
 909       pi0Pt      = (iLorentz+jLorentz).E();
 910       pi0InvMass = (iLorentz+jLorentz).M();
 911       radius     = stack->Particle(iPi0Indx)->R();
 912       ((TH2D*)listMC->FindObject("hMCInclusivePi0_PtRadius_Reco"))->Fill(pi0Pt, radius);
 913       ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtRadius_Reco_cent%d", cent)))->Fill(pi0Pt, radius);
 914
 915       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_Reco", pName.Data())))->Fill(pi0Pt, radius);
 916       ((TH2D*)listMC->FindObject(Form("hMC%s_PtRadius_Reco_cent%d", pName.Data(), cent)))->Fill(pi0Pt, radius);
 917
 918       for (Int_t iPID=0; iPID<kPIDs; iPID++)
 919         if ((iPIDBit & jPIDBit & PIDBIT[iPID])==PIDBIT[iPID]) {
 920           ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtInvMass_%s_Reco", namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
 921           ((TH2D*)listMC->FindObject(Form("hMCInclusivePi0_PtInvMass_%s_Reco_cent%d", namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
 922
 923           ((TH2D*)listMC->FindObject(Form("hMC%s_PtInvMass_%s_Reco", pName.Data(), namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
 924           ((TH2D*)listMC->FindObject(Form("hMC%s_PtInvMass_%s_Reco_cent%d", pName.Data(), namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
 925
 926           if (isConverted_i || isConverted_j) {
 927             ((TH2D*)listMC->FindObject(Form("hMCInclusiveCvtedPi0_PtInvMass_%s_Reco", namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
 928             ((TH2D*)listMC->FindObject(Form("hMCInclusiveCvtedPi0_PtInvMass_%s_Reco_cent%d", namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
 929           } else {
 930             ((TH2D*)listMC->FindObject(Form("hMCInclusivePurePi0_PtInvMass_%s_Reco", namePID[iPID].Data())))->Fill(pi0Pt, pi0InvMass);
 931             ((TH2D*)listMC->FindObject(Form("hMCInclusivePurePi0_PtInvMass_%s_Reco_cent%d", namePID[iPID].Data(), cent)))->Fill(pi0Pt, pi0InvMass);
 932           }
 933         }
 934     } // End loop calo cluster j
 935   } // End loop calo cluster i
 936
 937   return;
 938 }
 939
 940 //________________________________________________________________________
 941 TString AliPHOSpPbPi0Header::ClassifyMCPi0(Int_t index, AliStack* const stack)
 942 {
 943
 944     TString  pName;
 945     TParticle* pMC  = stack->Particle(index);
 946     Int_t      mpdg = 0;
 947
 948     if (pMC->GetFirstMother()>-1) {
 949       TParticle* pMother = stack->Particle(pMC->GetFirstMother());
 950       mpdg = TMath::Abs(pMother->GetPdgCode());
 951     }
 952
 953     // Classify Pi0 sources
 954     if (index<stack->GetNprimary())                              pName = "PrimaryPi0";
 955     else if (mpdg==310)                                          pName = "2ndPi0FromK0s";
 956     else if (mpdg==321 || mpdg==211 || mpdg==2212 || mpdg==2112) pName = "2ndPi0FromMaterials";     // K+-, pi+-, p(pbar), n(nbar)
 957     else                                                         pName = "2ndPi0FromOtherDecays";   // all other secondary pi0s
 958
 959     return pName;
 960 }
 961
 962 //________________________________________________________________________
 963 Int_t AliPHOSpPbPi0Header::HitPHOSModule(TParticle* const pMC, AliPHOSGeoUtils* const phosGeo)
 964 {
 965
 966   Int_t mod=0, relID[4]={0,0,0,0};
 967   Double_t x=0., z=0.;
 968   Double_t vtx[3] = { pMC->Vx(), pMC->Vy(), pMC->Vz() };
 969   Double_t theta = pMC->Theta();
 970   Double_t phi   = pMC->Phi();
 971
 972   if (!phosGeo->ImpactOnEmc(vtx, theta, phi, mod, z, x)) return -1;
 973   phosGeo->RelPosToRelId(mod, x, z, relID);
 974   if (fgUseFiducialCut) {
 975     const Int_t edgeX = 2;
 976     const Int_t edgeZ = 2;
 977     if (relID[2] >edgeX && relID[2]<(65-edgeX) && relID[3]>edgeZ && relID[3] <(57-edgeZ)) return relID[0];
 978     else return -1;
 979   }
 980
 981   return relID[0];
 982 }