PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, 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 // Class for the analysis of high pT pi0 event by event
  18 // Pi0/Eta identified by one of the following:
  19 //  -Invariant mass of 2 cluster in calorimeter
  20 //  -Shower shape analysis in calorimeter
  21 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
  22 //
  23 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TList.h>
  29 #include <TClonesArray.h>
  30 #include <TObjString.h>
  31
  32 // --- Analysis system ---
  33 #include "AliAnaPi0EbE.h"
  34 #include "AliCaloTrackReader.h"
  35 #include "AliIsolationCut.h"
  36 #include "AliNeutralMesonSelection.h"
  37 #include "AliCaloPID.h"
  38 #include "AliMCAnalysisUtils.h"
  39 #include "AliStack.h"
  40 #include "AliFiducialCut.h"
  41 #include "TParticle.h"
  42 #include "AliVCluster.h"
  43 #include "AliESDEvent.h"
  44 #include "AliAODEvent.h"
  45 #include "AliAODMCParticle.h"
  46
  47 ClassImp(AliAnaPi0EbE)
  48
  49 //____________________________
  50 AliAnaPi0EbE::AliAnaPi0EbE() :
  51     AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
  52     fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),
  53     fNLMCutMin(-1),                fNLMCutMax(10),
  54     fTimeCutMin(-10000),           fTimeCutMax(10000),
  55     fRejectTrackMatch(kTRUE),
  56     fFillPileUpHistograms(0),
  57     fFillWeightHistograms(kFALSE), fFillTMHisto(0),
  58     fFillSelectClHisto(0),         fFillOnlySimpleSSHisto(1),    fFillEMCALBCHistograms(0),
  59     fInputAODGammaConvName(""),
  60     // Histograms
  61     fhPt(0),                       fhE(0),
  62     fhEEta(0),                     fhEPhi(0),
  63     fhPtEta(0),                    fhPtPhi(0),                   fhEtaPhi(0),
  64     fhEtaPhiEMCALBC0(0),           fhEtaPhiEMCALBC1(0),          fhEtaPhiEMCALBCN(0),
  65     fhPtCentrality(),              fhPtEventPlane(0),
  66     fhPtReject(0),                 fhEReject(0),
  67     fhEEtaReject(0),               fhEPhiReject(0),              fhEtaPhiReject(0),
  68     fhMass(0),                     fhMassPt(0),                  fhMassSplitPt(0),
  69     fhSelectedMass(0),             fhSelectedMassPt(0),          fhSelectedMassSplitPt(0),
  70     fhAsymmetry(0),                fhSelectedAsymmetry(0),
  71     fhSplitE(0),                   fhSplitPt(0),
  72     fhSplitPtEta(0),               fhSplitPtPhi(0),
  73     fhNLocMaxSplitPt(0),
  74     fhPtDecay(0),                  fhEDecay(0),
  75     // Shower shape histos
  76     fhEDispersion(0),              fhELambda0(0),                fhELambda1(0),
  77     fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),
  78     fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(0),
  79     fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
  80     fhDispEtaE(0),                 fhDispPhiE(0),
  81     fhSumEtaE(0),                  fhSumPhiE(0),                 fhSumEtaPhiE(0),
  82     fhDispEtaPhiDiffE(0),          fhSphericityE(0),
  83
  84     // MC histos
  85     fhMCE(),                       fhMCPt(),
  86     fhMCPhi(),                     fhMCEta(),
  87     fhMCEReject(),                 fhMCPtReject(),
  88     fhMCPtCentrality(),
  89     fhMCPi0PtGenRecoFraction(0),   fhMCEtaPtGenRecoFraction(0),
  90     fhMCPi0DecayPt(0),             fhMCPi0DecayPtFraction(0),
  91     fhMCEtaDecayPt(0),             fhMCEtaDecayPtFraction(0),
  92     fhMCOtherDecayPt(0),
  93     fhMassPairMCPi0(0),            fhMassPairMCEta(0),
  94     fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
  95     // Weight studies
  96     fhECellClusterRatio(0),        fhECellClusterLogRatio(0),
  97     fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
  98     fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
  99     fhTrackMatchedMCParticleE(0),
 100     fhTrackMatchedMCParticleDEta(0), fhTrackMatchedMCParticleDPhi(0),
 101     fhdEdx(0),                     fhEOverP(0),                 fhEOverPNoTRD(0),
 102     // Number of local maxima in cluster
 103     fhNLocMaxE(0),                 fhNLocMaxPt(0),
 104     // PileUp
 105     fhTimeENoCut(0),                    fhTimeESPD(0),           fhTimeESPDMulti(0),
 106     fhTimeNPileUpVertSPD(0),            fhTimeNPileUpVertTrack(0),
 107     fhTimeNPileUpVertContributors(0),
 108     fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
 109 {
 110   //default ctor
 111
 112   for(Int_t i = 0; i < 6; i++)
 113   {
 114     fhMCE              [i] = 0;
 115     fhMCPt             [i] = 0;
 116     fhMCNLocMaxPt      [i] = 0;
 117     fhMCPhi            [i] = 0;
 118     fhMCEta            [i] = 0;
 119     fhMCPtCentrality   [i] = 0;
 120
 121     fhMCSplitE         [i] = 0;
 122     fhMCSplitPt        [i] = 0;
 123     fhMCSplitPtPhi     [i] = 0;
 124     fhMCSplitPtEta     [i] = 0;
 125     fhMCNLocMaxSplitPt [i] = 0;
 126
 127     fhEMCLambda0       [i] = 0;
 128     fhEMCLambda0NoTRD  [i] = 0;
 129     fhEMCLambda0FracMaxCellCut[i]= 0;
 130     fhEMCFracMaxCell   [i] = 0;
 131     fhEMCLambda1       [i] = 0;
 132     fhEMCDispersion    [i] = 0;
 133
 134     fhMCEDispEta       [i] = 0;
 135     fhMCEDispPhi       [i] = 0;
 136     fhMCESumEtaPhi     [i] = 0;
 137     fhMCEDispEtaPhiDiff[i] = 0;
 138     fhMCESphericity    [i] = 0;
 139     fhMCEAsymmetry     [i] = 0;
 140
 141     fhMCMassPt             [i]=0;
 142     fhMCMassSplitPt        [i]=0;
 143     fhMCSelectedMassPt     [i]=0;
 144     fhMCSelectedMassSplitPt[i]=0;
 145
 146     for(Int_t j = 0; j < 7; j++)
 147     {
 148       fhMCLambda0DispEta    [j][i] = 0;
 149       fhMCLambda0DispPhi    [j][i] = 0;
 150       fhMCDispEtaDispPhi    [j][i] = 0;
 151       fhMCAsymmetryLambda0  [j][i] = 0;
 152       fhMCAsymmetryDispEta  [j][i] = 0;
 153       fhMCAsymmetryDispPhi  [j][i] = 0;
 154     }
 155   }
 156
 157   for(Int_t j = 0; j < 7; j++)
 158   {
 159     fhLambda0DispEta    [j] = 0;
 160     fhLambda0DispPhi    [j] = 0;
 161     fhDispEtaDispPhi    [j] = 0;
 162     fhAsymmetryLambda0  [j] = 0;
 163     fhAsymmetryDispEta  [j] = 0;
 164     fhAsymmetryDispPhi  [j] = 0;
 165
 166     fhPtPi0PileUp       [j] = 0;
 167   }
 168
 169   for(Int_t i = 0; i < 3; i++)
 170   {
 171     fhELambda0LocMax       [i] = 0;
 172     fhELambda1LocMax       [i] = 0;
 173     fhEDispersionLocMax    [i] = 0;
 174     fhEDispEtaLocMax       [i] = 0;
 175     fhEDispPhiLocMax       [i] = 0;
 176     fhESumEtaPhiLocMax     [i] = 0;
 177     fhEDispEtaPhiDiffLocMax[i] = 0;
 178     fhESphericityLocMax    [i] = 0;
 179     fhEAsymmetryLocMax     [i] = 0;
 180   }
 181
 182   //Weight studies
 183   for(Int_t i =0; i < 14; i++){
 184     fhLambda0ForW0[i] = 0;
 185     //fhLambda1ForW0[i] = 0;
 186     if(i<8)fhMassPairLocMax[i] = 0;
 187   }
 188
 189   for(Int_t i = 0; i < 11; i++)
 190   {
 191     fhEtaPhiTriggerEMCALBC       [i] = 0 ;
 192     fhTimeTriggerEMCALBC         [i] = 0 ;
 193     fhTimeTriggerEMCALBCPileUpSPD[i] = 0 ;
 194
 195     fhEtaPhiTriggerEMCALBCUM     [i] = 0 ;
 196     fhTimeTriggerEMCALBCUM       [i] = 0 ;
 197
 198   }
 199
 200   //Initialize parameters
 201   InitParameters();
 202
 203 }
 204
 205 //_______________________________________________________________________________
 206 void AliAnaPi0EbE::FillPileUpHistograms(const Float_t energy, const Float_t time)
 207 {
 208   // Fill some histograms to understand pile-up
 209   if(!fFillPileUpHistograms) return;
 210
 211   //printf("E %f, time %f\n",energy,time);
 212   AliVEvent * event = GetReader()->GetInputEvent();
 213
 214   fhTimeENoCut->Fill(energy,time);
 215   if(GetReader()->IsPileUpFromSPD())     fhTimeESPD     ->Fill(energy,time);
 216   if(event->IsPileupFromSPDInMultBins()) fhTimeESPDMulti->Fill(energy,time);
 217
 218   if(energy < 8) return; // Fill time figures for high energy clusters not too close to trigger threshold
 219
 220   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 221   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 222
 223   // N pile up vertices
 224   Int_t nVerticesSPD    = -1;
 225   Int_t nVerticesTracks = -1;
 226
 227   if      (esdEv)
 228   {
 229     nVerticesSPD    = esdEv->GetNumberOfPileupVerticesSPD();
 230     nVerticesTracks = esdEv->GetNumberOfPileupVerticesTracks();
 231
 232   }//ESD
 233   else if (aodEv)
 234   {
 235     nVerticesSPD    = aodEv->GetNumberOfPileupVerticesSPD();
 236     nVerticesTracks = aodEv->GetNumberOfPileupVerticesTracks();
 237   }//AOD
 238
 239   fhTimeNPileUpVertSPD  ->Fill(time,nVerticesSPD);
 240   fhTimeNPileUpVertTrack->Fill(time,nVerticesTracks);
 241
 242   //printf("Is SPD %d, Is SPD Multi %d, n spd %d, n track %d\n",
 243   //       GetReader()->IsPileUpFromSPD(),event->IsPileupFromSPDInMultBins(),nVerticesSPD,nVerticesTracks);
 244
 245   Int_t ncont = -1;
 246   Float_t z1 = -1, z2 = -1;
 247   Float_t diamZ = -1;
 248   for(Int_t iVert=0; iVert<nVerticesSPD;iVert++)
 249   {
 250     if      (esdEv)
 251     {
 252       const AliESDVertex* pv=esdEv->GetPileupVertexSPD(iVert);
 253       ncont=pv->GetNContributors();
 254       z1 = esdEv->GetPrimaryVertexSPD()->GetZ();
 255       z2 = pv->GetZ();
 256       diamZ = esdEv->GetDiamondZ();
 257     }//ESD
 258     else if (aodEv)
 259     {
 260       AliAODVertex *pv=aodEv->GetVertex(iVert);
 261       if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
 262       ncont=pv->GetNContributors();
 263       z1=aodEv->GetPrimaryVertexSPD()->GetZ();
 264       z2=pv->GetZ();
 265       diamZ = aodEv->GetDiamondZ();
 266     }// AOD
 267
 268     Double_t distZ  = TMath::Abs(z2-z1);
 269     diamZ  = TMath::Abs(z2-diamZ);
 270
 271     fhTimeNPileUpVertContributors  ->Fill(time,ncont);
 272     fhTimePileUpMainVertexZDistance->Fill(time,distZ);
 273     fhTimePileUpMainVertexZDiamond ->Fill(time,diamZ);
 274
 275   }// loop
 276 }
 277
 278
 279 //___________________________________________________________________________________________
 280 void AliAnaPi0EbE::FillRejectedClusterHistograms(const TLorentzVector mom, const Int_t mctag)
 281 {
 282  // Fill histograms that do not pass the identification (SS case only)
 283
 284   Float_t ener  = mom.E();
 285   Float_t pt    = mom.Pt();
 286   Float_t phi   = mom.Phi();
 287   if(phi < 0) phi+=TMath::TwoPi();
 288   Float_t eta = mom.Eta();
 289
 290   fhPtReject     ->Fill(pt);
 291   fhEReject      ->Fill(ener);
 292
 293   fhEEtaReject   ->Fill(ener,eta);
 294   fhEPhiReject   ->Fill(ener,phi);
 295   fhEtaPhiReject ->Fill(eta,phi);
 296
 297   if(IsDataMC())
 298   {
 299     Int_t mcIndex = GetMCIndex(mctag);
 300     fhMCEReject  [mcIndex] ->Fill(ener);
 301     fhMCPtReject [mcIndex] ->Fill(pt);
 302   }
 303 }
 304
 305 //_____________________________________________________________________________________
 306 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
 307                                                  const Int_t nMaxima,
 308                                                  const Int_t tag,
 309                                                  const Float_t asy)
 310 {
 311   // Fill shower shape, timing and other histograms for selected clusters from decay
 312
 313   Float_t e    = cluster->E();
 314   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
 315   Float_t l0   = cluster->GetM02();
 316   Float_t l1   = cluster->GetM20();
 317   Int_t   nSM  = GetModuleNumber(cluster);
 318
 319   Int_t ebin = -1;
 320   if      (e < 2 ) ebin = 0;
 321   else if (e < 4 ) ebin = 1;
 322   else if (e < 6 ) ebin = 2;
 323   else if (e < 10) ebin = 3;
 324   else if (e < 15) ebin = 4;
 325   else if (e < 20) ebin = 5;
 326   else             ebin = 6;
 327
 328   Int_t indexMax = -1;
 329   if     (nMaxima==1) indexMax = 0 ;
 330   else if(nMaxima==2) indexMax = 1 ;
 331   else                indexMax = 2 ;
 332
 333
 334   AliVCaloCells * cell = 0x0;
 335   if(fCalorimeter == "PHOS")
 336     cell = GetPHOSCells();
 337   else
 338     cell = GetEMCALCells();
 339
 340   Float_t maxCellFraction = 0;
 341   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
 342   fhEFracMaxCell->Fill(e,maxCellFraction);
 343
 344   FillWeightHistograms(cluster);
 345
 346   fhEDispersion->Fill(e, disp);
 347   fhELambda0   ->Fill(e, l0  );
 348   fhELambda1   ->Fill(e, l1  );
 349
 350   Float_t ll0  = 0., ll1  = 0.;
 351   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
 352   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 353   if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
 354   {
 355     GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 356                                                                                  ll0, ll1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
 357
 358     fhDispEtaE        -> Fill(e,dEta);
 359     fhDispPhiE        -> Fill(e,dPhi);
 360     fhSumEtaE         -> Fill(e,sEta);
 361     fhSumPhiE         -> Fill(e,sPhi);
 362     fhSumEtaPhiE      -> Fill(e,sEtaPhi);
 363     fhDispEtaPhiDiffE -> Fill(e,dPhi-dEta);
 364     if(dEta+dPhi>0)fhSphericityE -> Fill(e,(dPhi-dEta)/(dEta+dPhi));
 365
 366     fhDispEtaDispPhi[ebin]->Fill(dEta,dPhi);
 367     fhLambda0DispEta[ebin]->Fill(l0  ,dEta);
 368     fhLambda0DispPhi[ebin]->Fill(l0  ,dPhi);
 369
 370     if (fAnaType==kSSCalo)
 371     {
 372       // Asymmetry histograms
 373       fhAsymmetryLambda0[ebin]->Fill(l0 ,asy);
 374       fhAsymmetryDispEta[ebin]->Fill(dEta,asy);
 375       fhAsymmetryDispPhi[ebin]->Fill(dPhi,asy);
 376     }
 377   }
 378
 379   fhNLocMaxE ->Fill(e ,nMaxima);
 380
 381   fhELambda0LocMax   [indexMax]->Fill(e,l0);
 382   fhELambda1LocMax   [indexMax]->Fill(e,l1);
 383   fhEDispersionLocMax[indexMax]->Fill(e,disp);
 384
 385   if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
 386   {
 387     fhEDispEtaLocMax       [indexMax]-> Fill(e,dEta);
 388     fhEDispPhiLocMax       [indexMax]-> Fill(e,dPhi);
 389     fhESumEtaPhiLocMax     [indexMax]-> Fill(e,sEtaPhi);
 390     fhEDispEtaPhiDiffLocMax[indexMax]-> Fill(e,dPhi-dEta);
 391     if(dEta+dPhi>0)       fhESphericityLocMax[indexMax]->Fill(e,(dPhi-dEta)/(dEta+dPhi));
 392     if(fAnaType==kSSCalo) fhEAsymmetryLocMax [indexMax]->Fill(e  ,asy);
 393
 394   }
 395
 396   if(fCalorimeter=="EMCAL" && nSM < 6)
 397   {
 398     fhELambda0NoTRD->Fill(e, l0  );
 399     fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
 400   }
 401
 402   if(maxCellFraction < 0.5)
 403     fhELambda0FracMaxCellCut->Fill(e, l0  );
 404
 405   fhETime  ->Fill(e, cluster->GetTOF()*1.e9);
 406   fhENCells->Fill(e, cluster->GetNCells());
 407
 408   // Fill Track matching control histograms
 409   if(fFillTMHisto)
 410   {
 411     Float_t dZ  = cluster->GetTrackDz();
 412     Float_t dR  = cluster->GetTrackDx();
 413
 414     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 415     {
 416       dR = 2000., dZ = 2000.;
 417       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 418     }
 419     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
 420
 421     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
 422     {
 423       fhTrackMatchedDEta->Fill(e,dZ);
 424       fhTrackMatchedDPhi->Fill(e,dR);
 425       if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
 426     }
 427
 428     // Check dEdx and E/p of matched clusters
 429
 430     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 431     {
 432       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
 433
 434       if(track)
 435       {
 436         Float_t dEdx = track->GetTPCsignal();
 437         fhdEdx->Fill(e, dEdx);
 438
 439         Float_t eOverp = e/track->P();
 440         fhEOverP->Fill(e,  eOverp);
 441
 442         if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e,  eOverp);
 443
 444       }
 445       //else
 446       //  printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
 447
 448
 449
 450       if(IsDataMC())
 451       {
 452         Float_t mctag = -1;
 453         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
 454         {
 455           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 456                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  2.5 ;
 457           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  0.5 ;
 458           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  1.5 ;
 459           else                                                                                 mctag =  3.5 ;
 460
 461         }
 462         else
 463         {
 464           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 465                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) mctag =  6.5 ;
 466           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) mctag =  4.5 ;
 467           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) mctag =  5.5 ;
 468           else                                                                                 mctag =  7.5 ;
 469         }
 470
 471         fhTrackMatchedMCParticleE   ->Fill(e , mctag);
 472         fhTrackMatchedMCParticleDEta->Fill(dZ, mctag);
 473         fhTrackMatchedMCParticleDPhi->Fill(dR, mctag);
 474
 475       }  // MC
 476     }
 477   }// Track matching histograms
 478
 479   if(IsDataMC())
 480   {
 481     Int_t mcIndex = GetMCIndex(tag);
 482
 483     fhEMCLambda0[mcIndex]    ->Fill(e, l0);
 484     fhEMCLambda1[mcIndex]    ->Fill(e, l1);
 485     fhEMCDispersion[mcIndex] ->Fill(e, disp);
 486     fhEMCFracMaxCell[mcIndex]->Fill(e,maxCellFraction);
 487
 488     if(fCalorimeter=="EMCAL" && nSM < 6)
 489       fhEMCLambda0NoTRD[mcIndex]->Fill(e, l0  );
 490
 491     if(maxCellFraction < 0.5)
 492       fhEMCLambda0FracMaxCellCut[mcIndex]->Fill(e, l0  );
 493
 494     if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
 495     {
 496       fhMCEDispEta        [mcIndex]-> Fill(e,dEta);
 497       fhMCEDispPhi        [mcIndex]-> Fill(e,dPhi);
 498       fhMCESumEtaPhi      [mcIndex]-> Fill(e,sEtaPhi);
 499       fhMCEDispEtaPhiDiff [mcIndex]-> Fill(e,dPhi-dEta);
 500       if(dEta+dPhi>0)fhMCESphericity[mcIndex]-> Fill(e,(dPhi-dEta)/(dEta+dPhi));
 501
 502       if (fAnaType==kSSCalo)
 503       {
 504         fhMCAsymmetryLambda0[ebin][mcIndex]->Fill(l0 ,asy);
 505         fhMCAsymmetryDispEta[ebin][mcIndex]->Fill(dEta,asy);
 506         fhMCAsymmetryDispPhi[ebin][mcIndex]->Fill(dPhi,asy);
 507       }
 508
 509       fhMCDispEtaDispPhi[ebin][mcIndex]->Fill(dEta,dPhi);
 510       fhMCLambda0DispEta[ebin][mcIndex]->Fill(l0  ,dEta);
 511       fhMCLambda0DispPhi[ebin][mcIndex]->Fill(l0  ,dPhi);
 512
 513     }
 514
 515   }//MC
 516
 517 }
 518
 519 //________________________________________________________
 520 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
 521 {
 522   // Calculate weights and fill histograms
 523
 524   if(!fFillWeightHistograms || GetMixedEvent()) return;
 525
 526   AliVCaloCells* cells = 0;
 527   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 528   else                        cells = GetPHOSCells();
 529
 530   // First recalculate energy in case non linearity was applied
 531   Float_t  energy = 0;
 532   Float_t  ampMax = 0;
 533   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 534   {
 535
 536     Int_t id       = clus->GetCellsAbsId()[ipos];
 537
 538     //Recalibrate cell energy if needed
 539     Float_t amp = cells->GetCellAmplitude(id);
 540     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 541
 542     energy    += amp;
 543
 544     if(amp> ampMax)
 545       ampMax = amp;
 546
 547   } // energy loop
 548
 549   if(energy <=0 )
 550   {
 551     printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
 552     return;
 553   }
 554
 555   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
 556   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
 557
 558   //Get the ratio and log ratio to all cells in cluster
 559   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 560   {
 561     Int_t id       = clus->GetCellsAbsId()[ipos];
 562
 563     //Recalibrate cell energy if needed
 564     Float_t amp = cells->GetCellAmplitude(id);
 565     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 566
 567     fhECellClusterRatio   ->Fill(energy,amp/energy);
 568     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
 569   }
 570
 571   //Recalculate shower shape for different W0
 572   if(fCalorimeter=="EMCAL"){
 573
 574     Float_t l0org = clus->GetM02();
 575     Float_t l1org = clus->GetM20();
 576     Float_t dorg  = clus->GetDispersion();
 577
 578     for(Int_t iw = 0; iw < 14; iw++)
 579     {
 580       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
 581       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
 582
 583       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
 584       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
 585
 586     } // w0 loop
 587
 588     // Set the original values back
 589     clus->SetM02(l0org);
 590     clus->SetM20(l1org);
 591     clus->SetDispersion(dorg);
 592
 593   }// EMCAL
 594 }
 595
 596 //__________________________________________
 597 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
 598 {
 599         //Save parameters used for analysis
 600   TString parList ; //this will be list of parameters used for this analysis.
 601   const Int_t buffersize = 255;
 602   char onePar[buffersize] ;
 603
 604   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
 605   parList+=onePar ;
 606   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
 607   parList+=onePar ;
 608
 609   if(fAnaType == kSSCalo)
 610   {
 611     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 612     parList+=onePar ;
 613     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 614     parList+=onePar ;
 615     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 616     parList+=onePar ;
 617     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 618     parList+=onePar ;
 619   }
 620
 621   //Get parameters set in base class.
 622   parList += GetBaseParametersList() ;
 623
 624   //Get parameters set in PID class.
 625   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
 626
 627   return new TObjString(parList) ;
 628 }
 629
 630 //_____________________________________________
 631 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 632 {
 633   // Create histograms to be saved in output file and
 634   // store them in outputContainer
 635   TList * outputContainer = new TList() ;
 636   outputContainer->SetName("Pi0EbEHistos") ;
 637
 638   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 639   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 640   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 641   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
 642   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
 643   Int_t tbins    = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin   = GetHistogramRanges()->GetHistoTimeMin();
 644   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax   = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin   = GetHistogramRanges()->GetHistoNClusterCellMin();
 645
 646   Int_t   nmassbins   = GetHistogramRanges()->GetHistoMassBins();
 647   Float_t massmin     = GetHistogramRanges()->GetHistoMassMin();
 648   Float_t massmax     = GetHistogramRanges()->GetHistoMassMax();
 649
 650   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 651   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 652   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 653   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 654   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 655   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 656
 657   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
 658   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
 659   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
 660   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
 661   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
 662   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
 663
 664   Int_t   ntimebins= GetHistogramRanges()->GetHistoTimeBins();
 665   Float_t timemax  = GetHistogramRanges()->GetHistoTimeMax();
 666   Float_t timemin  = GetHistogramRanges()->GetHistoTimeMin();
 667
 668   TString nlm[]   ={"1 Local Maxima","2 Local Maxima", "NLM > 2"};
 669   TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 670   TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 671   Int_t   bin[]   = {0,2,4,6,10,15,20,100}; // energy bins
 672
 673   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
 674   fhPt->SetYTitle("N");
 675   fhPt->SetXTitle("p_{T} (GeV/c)");
 676   outputContainer->Add(fhPt) ;
 677
 678   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
 679   fhE->SetYTitle("N");
 680   fhE->SetXTitle("E (GeV)");
 681   outputContainer->Add(fhE) ;
 682
 683   fhEPhi  = new TH2F
 684   ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 685   fhEPhi->SetYTitle("#phi (rad)");
 686   fhEPhi->SetXTitle("E (GeV)");
 687   outputContainer->Add(fhEPhi) ;
 688
 689   fhEEta  = new TH2F
 690   ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 691   fhEEta->SetYTitle("#eta");
 692   fhEEta->SetXTitle("E (GeV)");
 693   outputContainer->Add(fhEEta) ;
 694
 695   fhPtPhi  = new TH2F
 696   ("hPtPhi","Selected #pi^{0} (#eta) pairs: p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 697   fhPtPhi->SetYTitle("#phi (rad)");
 698   fhPtPhi->SetXTitle("p_{T} (GeV/c)");
 699   outputContainer->Add(fhPtPhi) ;
 700
 701   fhPtEta  = new TH2F
 702   ("hPtEta","Selected #pi^{0} (#eta) pairs: p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 703   fhPtEta->SetYTitle("#eta");
 704   fhPtEta->SetXTitle("p_{T} (GeV/c)");
 705   outputContainer->Add(fhPtEta) ;
 706
 707   fhEtaPhi  = new TH2F
 708   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
 709   fhEtaPhi->SetYTitle("#phi (rad)");
 710   fhEtaPhi->SetXTitle("#eta");
 711   outputContainer->Add(fhEtaPhi) ;
 712
 713   if(fCalorimeter=="EMCAL" && fFillEMCALBCHistograms)
 714   {
 715     fhEtaPhiEMCALBC0  = new TH2F
 716     ("hEtaPhiEMCALBC0","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| < 25 ns, EMCAL-BC=0",netabins,etamin,etamax,nphibins,phimin,phimax);
 717     fhEtaPhiEMCALBC0->SetYTitle("#phi (rad)");
 718     fhEtaPhiEMCALBC0->SetXTitle("#eta");
 719     outputContainer->Add(fhEtaPhiEMCALBC0) ;
 720
 721     fhEtaPhiEMCALBC1  = new TH2F
 722     ("hEtaPhiEMCALBC1","cluster,E > 2 GeV, #eta vs #phi, for clusters with 25 < |time| < 75 ns, EMCAL-BC=1",netabins,etamin,etamax,nphibins,phimin,phimax);
 723     fhEtaPhiEMCALBC1->SetYTitle("#phi (rad)");
 724     fhEtaPhiEMCALBC1->SetXTitle("#eta");
 725     outputContainer->Add(fhEtaPhiEMCALBC1) ;
 726
 727     fhEtaPhiEMCALBCN  = new TH2F
 728     ("hEtaPhiEMCALBCN","cluster,E > 2 GeV, #eta vs #phi, for clusters with |time| > 75 ns, EMCAL-BC>1",netabins,etamin,etamax,nphibins,phimin,phimax);
 729     fhEtaPhiEMCALBCN->SetYTitle("#phi (rad)");
 730     fhEtaPhiEMCALBCN->SetXTitle("#eta");
 731     outputContainer->Add(fhEtaPhiEMCALBCN) ;
 732
 733     for(Int_t i = 0; i < 11; i++)
 734     {
 735       fhEtaPhiTriggerEMCALBC[i] = new TH2F
 736       (Form("hEtaPhiTriggerEMCALBC%d",i-5),
 737        Form("meson E > 2 GeV, #eta vs #phi, Trigger EMCAL-BC=%d",i-5),
 738        netabins,etamin,etamax,nphibins,phimin,phimax);
 739       fhEtaPhiTriggerEMCALBC[i]->SetYTitle("#phi (rad)");
 740       fhEtaPhiTriggerEMCALBC[i]->SetXTitle("#eta");
 741       outputContainer->Add(fhEtaPhiTriggerEMCALBC[i]) ;
 742
 743       fhTimeTriggerEMCALBC[i] = new TH2F
 744       (Form("hTimeTriggerEMCALBC%d",i-5),
 745        Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
 746        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 747       fhTimeTriggerEMCALBC[i]->SetXTitle("E (GeV)");
 748       fhTimeTriggerEMCALBC[i]->SetYTitle("time (ns)");
 749       outputContainer->Add(fhTimeTriggerEMCALBC[i]);
 750
 751       fhTimeTriggerEMCALBCPileUpSPD[i] = new TH2F
 752       (Form("hTimeTriggerEMCALBC%dPileUpSPD",i-5),
 753        Form("meson time vs E, Trigger EMCAL-BC=%d",i-5),
 754        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 755       fhTimeTriggerEMCALBCPileUpSPD[i]->SetXTitle("E (GeV)");
 756       fhTimeTriggerEMCALBCPileUpSPD[i]->SetYTitle("time (ns)");
 757       outputContainer->Add(fhTimeTriggerEMCALBCPileUpSPD[i]);
 758
 759       fhEtaPhiTriggerEMCALBCUM[i] = new TH2F
 760       (Form("hEtaPhiTriggerEMCALBC%d_UnMatch",i-5),
 761        Form("meson E > 2 GeV, #eta vs #phi, unmatched trigger EMCAL-BC=%d",i-5),
 762        netabins,etamin,etamax,nphibins,phimin,phimax);
 763       fhEtaPhiTriggerEMCALBCUM[i]->SetYTitle("#phi (rad)");
 764       fhEtaPhiTriggerEMCALBCUM[i]->SetXTitle("#eta");
 765       outputContainer->Add(fhEtaPhiTriggerEMCALBCUM[i]) ;
 766
 767       fhTimeTriggerEMCALBCUM[i] = new TH2F
 768       (Form("hTimeTriggerEMCALBC%d_UnMatch",i-5),
 769        Form("meson time vs E, unmatched trigger EMCAL-BC=%d",i-5),
 770        nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
 771       fhTimeTriggerEMCALBCUM[i]->SetXTitle("E (GeV)");
 772       fhTimeTriggerEMCALBCUM[i]->SetYTitle("time (ns)");
 773       outputContainer->Add(fhTimeTriggerEMCALBCUM[i]);
 774
 775     }
 776   }
 777
 778   fhPtCentrality  = new TH2F("hPtCentrality","centrality vs p_{T}",nptbins,ptmin,ptmax, 100,0,100);
 779   fhPtCentrality->SetYTitle("centrality");
 780   fhPtCentrality->SetXTitle("p_{T}(GeV/c)");
 781   outputContainer->Add(fhPtCentrality) ;
 782
 783   fhPtEventPlane  = new TH2F("hPtEventPlane","event plane angle vs p_{T}",nptbins,ptmin,ptmax, 100,0,TMath::Pi());
 784   fhPtEventPlane->SetYTitle("Event plane angle (rad)");
 785   fhPtEventPlane->SetXTitle("p_{T} (GeV/c)");
 786   outputContainer->Add(fhPtEventPlane) ;
 787
 788   if(fAnaType == kSSCalo)
 789   {
 790     fhPtReject  = new TH1F("hPtReject","Number of rejected as #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
 791     fhPtReject->SetYTitle("N");
 792     fhPtReject->SetXTitle("p_{T} (GeV/c)");
 793     outputContainer->Add(fhPtReject) ;
 794
 795     fhEReject  = new TH1F("hEReject","Number of rejected as  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
 796     fhEReject->SetYTitle("N");
 797     fhEReject->SetXTitle("E (GeV)");
 798     outputContainer->Add(fhEReject) ;
 799
 800     fhEPhiReject  = new TH2F
 801     ("hEPhiReject","Rejected #pi^{0} (#eta) cluster: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 802     fhEPhiReject->SetYTitle("#phi (rad)");
 803     fhEPhiReject->SetXTitle("E (GeV)");
 804     outputContainer->Add(fhEPhiReject) ;
 805
 806     fhEEtaReject  = new TH2F
 807     ("hEEtaReject","Rejected #pi^{0} (#eta) cluster: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 808     fhEEtaReject->SetYTitle("#eta");
 809     fhEEtaReject->SetXTitle("E (GeV)");
 810     outputContainer->Add(fhEEtaReject) ;
 811
 812     fhEtaPhiReject  = new TH2F
 813     ("hEtaPhiReject","Rejected #pi^{0} (#eta) cluster: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
 814     fhEtaPhiReject->SetYTitle("#phi (rad)");
 815     fhEtaPhiReject->SetXTitle("#eta");
 816     outputContainer->Add(fhEtaPhiReject) ;
 817   }
 818
 819   fhMass  = new TH2F
 820   ("hMass","all pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 821   fhMass->SetYTitle("mass (GeV/c^{2})");
 822   fhMass->SetXTitle("E (GeV)");
 823   outputContainer->Add(fhMass) ;
 824
 825   fhSelectedMass  = new TH2F
 826   ("hSelectedMass","Selected #pi^{0} (#eta) pairs mass: E vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 827   fhSelectedMass->SetYTitle("mass (GeV/c^{2})");
 828   fhSelectedMass->SetXTitle("E (GeV)");
 829   outputContainer->Add(fhSelectedMass) ;
 830
 831   fhMassPt  = new TH2F
 832   ("hMassPt","all pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 833   fhMassPt->SetYTitle("mass (GeV/c^{2})");
 834   fhMassPt->SetXTitle("p_{T} (GeV/c)");
 835   outputContainer->Add(fhMassPt) ;
 836
 837   fhSelectedMassPt  = new TH2F
 838   ("hSelectedMassPt","Selected #pi^{0} (#eta) pairs mass: p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
 839   fhSelectedMassPt->SetYTitle("mass (GeV/c^{2})");
 840   fhSelectedMassPt->SetXTitle("p_{T} (GeV/c)");
 841   outputContainer->Add(fhSelectedMassPt) ;
 842
 843   if(fAnaType != kSSCalo)
 844   {
 845     fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 846     fhPtDecay->SetYTitle("N");
 847     fhPtDecay->SetXTitle("p_{T} (GeV/c)");
 848     outputContainer->Add(fhPtDecay) ;
 849
 850     fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 851     fhEDecay->SetYTitle("N");
 852     fhEDecay->SetXTitle("E (GeV)");
 853     outputContainer->Add(fhEDecay) ;
 854   }
 855
 856   ////////
 857
 858   if( fFillSelectClHisto )
 859   {
 860
 861     fhEDispersion  = new TH2F
 862     ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 863     fhEDispersion->SetYTitle("D^{2}");
 864     fhEDispersion->SetXTitle("E (GeV)");
 865     outputContainer->Add(fhEDispersion) ;
 866
 867     fhELambda0  = new TH2F
 868     ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 869     fhELambda0->SetYTitle("#lambda_{0}^{2}");
 870     fhELambda0->SetXTitle("E (GeV)");
 871     outputContainer->Add(fhELambda0) ;
 872
 873     fhELambda1  = new TH2F
 874     ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 875     fhELambda1->SetYTitle("#lambda_{1}^{2}");
 876     fhELambda1->SetXTitle("E (GeV)");
 877     outputContainer->Add(fhELambda1) ;
 878
 879     fhELambda0FracMaxCellCut  = new TH2F
 880     ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 881     fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
 882     fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
 883     outputContainer->Add(fhELambda0FracMaxCellCut) ;
 884
 885     fhEFracMaxCell  = new TH2F
 886     ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
 887     fhEFracMaxCell->SetYTitle("Fraction");
 888     fhEFracMaxCell->SetXTitle("E (GeV)");
 889     outputContainer->Add(fhEFracMaxCell) ;
 890
 891     if(fCalorimeter=="EMCAL")
 892     {
 893       fhELambda0NoTRD  = new TH2F
 894       ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 895       fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
 896       fhELambda0NoTRD->SetXTitle("E (GeV)");
 897       outputContainer->Add(fhELambda0NoTRD) ;
 898
 899       fhEFracMaxCellNoTRD  = new TH2F
 900       ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
 901       fhEFracMaxCellNoTRD->SetYTitle("Fraction");
 902       fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
 903       outputContainer->Add(fhEFracMaxCellNoTRD) ;
 904
 905       if(!fFillOnlySimpleSSHisto)
 906       {
 907         fhDispEtaE  = new TH2F ("hDispEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 908         fhDispEtaE->SetXTitle("E (GeV)");
 909         fhDispEtaE->SetYTitle("#sigma^{2}_{#eta #eta}");
 910         outputContainer->Add(fhDispEtaE);
 911
 912         fhDispPhiE  = new TH2F ("hDispPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 913         fhDispPhiE->SetXTitle("E (GeV)");
 914         fhDispPhiE->SetYTitle("#sigma^{2}_{#phi #phi}");
 915         outputContainer->Add(fhDispPhiE);
 916
 917         fhSumEtaE  = new TH2F ("hSumEtaE","#sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",  nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 918         fhSumEtaE->SetXTitle("E (GeV)");
 919         fhSumEtaE->SetYTitle("#delta^{2}_{#eta #eta}");
 920         outputContainer->Add(fhSumEtaE);
 921
 922         fhSumPhiE  = new TH2F ("hSumPhiE","#sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",
 923                                nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 924         fhSumPhiE->SetXTitle("E (GeV)");
 925         fhSumPhiE->SetYTitle("#delta^{2}_{#phi #phi}");
 926         outputContainer->Add(fhSumPhiE);
 927
 928         fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",
 929                                   nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
 930         fhSumEtaPhiE->SetXTitle("E (GeV)");
 931         fhSumEtaPhiE->SetYTitle("#delta^{2}_{#eta #phi}");
 932         outputContainer->Add(fhSumEtaPhiE);
 933
 934         fhDispEtaPhiDiffE  = new TH2F ("hDispEtaPhiDiffE","#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",
 935                                        nptbins,ptmin,ptmax,200, -10,10);
 936         fhDispEtaPhiDiffE->SetXTitle("E (GeV)");
 937         fhDispEtaPhiDiffE->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
 938         outputContainer->Add(fhDispEtaPhiDiffE);
 939
 940         fhSphericityE  = new TH2F ("hSphericityE","(#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",
 941                                    nptbins,ptmin,ptmax, 200, -1,1);
 942         fhSphericityE->SetXTitle("E (GeV)");
 943         fhSphericityE->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
 944         outputContainer->Add(fhSphericityE);
 945
 946         for(Int_t i = 0; i < 7; i++)
 947         {
 948           fhDispEtaDispPhi[i] = new TH2F (Form("hDispEtaDispPhi_EBin%d",i),Form("#sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
 949                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
 950           fhDispEtaDispPhi[i]->SetXTitle("#sigma^{2}_{#eta #eta}");
 951           fhDispEtaDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
 952           outputContainer->Add(fhDispEtaDispPhi[i]);
 953
 954           fhLambda0DispEta[i] = new TH2F (Form("hLambda0DispEta_EBin%d",i),Form("#lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",bin[i],bin[i+1]),
 955                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
 956           fhLambda0DispEta[i]->SetXTitle("#lambda^{2}_{0}");
 957           fhLambda0DispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
 958           outputContainer->Add(fhLambda0DispEta[i]);
 959
 960           fhLambda0DispPhi[i] = new TH2F (Form("hLambda0DispPhi_EBin%d",i),Form("#lambda^{2}_{0}} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",bin[i],bin[i+1]),
 961                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
 962           fhLambda0DispPhi[i]->SetXTitle("#lambda^{2}_{0}");
 963           fhLambda0DispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
 964           outputContainer->Add(fhLambda0DispPhi[i]);
 965
 966         }
 967       }
 968     }
 969
 970     fhNLocMaxE = new TH2F("hNLocMaxE","Number of local maxima in cluster",
 971                           nptbins,ptmin,ptmax,10,0,10);
 972     fhNLocMaxE ->SetYTitle("N maxima");
 973     fhNLocMaxE ->SetXTitle("E (GeV)");
 974     outputContainer->Add(fhNLocMaxE) ;
 975
 976     if(fAnaType == kSSCalo)
 977     {
 978       fhNLocMaxPt = new TH2F("hNLocMaxPt","Number of local maxima in cluster",
 979                             nptbins,ptmin,ptmax,10,0,10);
 980       fhNLocMaxPt ->SetYTitle("N maxima");
 981       fhNLocMaxPt ->SetXTitle("p_{T} (GeV/c)");
 982       outputContainer->Add(fhNLocMaxPt) ;
 983     }
 984
 985     for (Int_t i = 0; i < 3; i++)
 986     {
 987       fhELambda0LocMax[i]  = new TH2F(Form("hELambda0LocMax%d",i+1),
 988                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, %s",nlm[i].Data()),
 989                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 990       fhELambda0LocMax[i]->SetYTitle("#lambda_{0}^{2}");
 991       fhELambda0LocMax[i]->SetXTitle("E (GeV)");
 992       outputContainer->Add(fhELambda0LocMax[i]) ;
 993
 994       fhELambda1LocMax[i]  = new TH2F(Form("hELambda1LocMax%d",i+1),
 995                                       Form("Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, %s",nlm[i].Data()),
 996                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 997       fhELambda1LocMax[i]->SetYTitle("#lambda_{1}^{2}");
 998       fhELambda1LocMax[i]->SetXTitle("E (GeV)");
 999       outputContainer->Add(fhELambda1LocMax[i]) ;
1000
1001       fhEDispersionLocMax[i]  = new TH2F(Form("hEDispersionLocMax%d",i+1),
1002                                          Form("Selected #pi^{0} (#eta) pairs: E vs dispersion^{2}, %s",nlm[i].Data()),
1003                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1004       fhEDispersionLocMax[i]->SetYTitle("dispersion^{2}");
1005       fhEDispersionLocMax[i]->SetXTitle("E (GeV)");
1006       outputContainer->Add(fhEDispersionLocMax[i]) ;
1007
1008       if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
1009       {
1010         fhEDispEtaLocMax[i]  = new TH2F(Form("hEDispEtaLocMax%d",i+1),
1011                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #eta}, %s",nlm[i].Data()),
1012                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1013         fhEDispEtaLocMax[i]->SetYTitle("#sigma_{#eta #eta}");
1014         fhEDispEtaLocMax[i]->SetXTitle("E (GeV)");
1015         outputContainer->Add(fhEDispEtaLocMax[i]) ;
1016
1017         fhEDispPhiLocMax[i]  = new TH2F(Form("hEDispPhiLocMax%d",i+1),
1018                                         Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi}, %s",nlm[i].Data()),
1019                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1020         fhEDispPhiLocMax[i]->SetYTitle("#sigma_{#phi #phi}");
1021         fhEDispPhiLocMax[i]->SetXTitle("E (GeV)");
1022         outputContainer->Add(fhEDispPhiLocMax[i]) ;
1023
1024         fhESumEtaPhiLocMax[i]  = new TH2F(Form("hESumEtaPhiLocMax%d",i+1),
1025                                           Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#eta #phi}, %s",nlm[i].Data()),
1026                                           nptbins,ptmin,ptmax,2*ssbins,-ssmax,ssmax);
1027         fhESumEtaPhiLocMax[i]->SetYTitle("#sigma_{#eta #phi}");
1028         fhESumEtaPhiLocMax[i]->SetXTitle("E (GeV)");
1029         outputContainer->Add(fhESumEtaPhiLocMax[i]) ;
1030
1031         fhEDispEtaPhiDiffLocMax[i]  = new TH2F(Form("hEDispEtaPhiDiffLocMax%d",i+1),
1032                                                Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta}, %s",nlm[i].Data()),
1033                                                nptbins,ptmin,ptmax,200, -10,10);
1034         fhEDispEtaPhiDiffLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta}");
1035         fhEDispEtaPhiDiffLocMax[i]->SetXTitle("E (GeV)");
1036         outputContainer->Add(fhEDispEtaPhiDiffLocMax[i]) ;
1037
1038         fhESphericityLocMax[i]  = new TH2F(Form("hESphericityLocMax%d",i+1),
1039                                            Form("Selected #pi^{0} (#eta) pairs: E vs #sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta}), %s",nlm[i].Data()),
1040                                            nptbins,ptmin,ptmax,200, -1,1);
1041         fhESphericityLocMax[i]->SetYTitle("#sigma_{#phi #phi} - #sigma_{#eta #eta} / (#sigma_{#phi #phi} + #sigma_{#eta #eta})");
1042         fhESphericityLocMax[i]->SetXTitle("E (GeV)");
1043         outputContainer->Add(fhESphericityLocMax[i]) ;
1044       }
1045
1046     }
1047
1048     fhENCells  = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
1049     fhENCells->SetXTitle("E (GeV)");
1050     fhENCells->SetYTitle("# of cells in cluster");
1051     outputContainer->Add(fhENCells);
1052
1053     fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
1054     fhETime->SetXTitle("E (GeV)");
1055     fhETime->SetYTitle("t (ns)");
1056     outputContainer->Add(fhETime);
1057
1058   }
1059
1060
1061     fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
1062     fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
1063     fhEPairDiffTime->SetYTitle("#Delta t (ns)");
1064     outputContainer->Add(fhEPairDiffTime);
1065
1066   if(fAnaType == kIMCalo)
1067   {
1068     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
1069     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
1070       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
1071       "2 Local Maxima paired with more than 2 Local Maxima",
1072       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
1073
1074     for (Int_t i = 0; i < 8 ; i++)
1075     {
1076
1077       if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
1078
1079       fhMassPairLocMax[i]  = new TH2F
1080       (Form("MassPairLocMax%s",combiName[i].Data()),
1081        Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
1082        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1083       fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
1084       fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
1085       outputContainer->Add(fhMassPairLocMax[i]) ;
1086     }
1087   }
1088
1089   if(fFillTMHisto)
1090   {
1091     fhTrackMatchedDEta  = new TH2F
1092     ("hTrackMatchedDEta",
1093      "d#eta of cluster-track vs cluster energy",
1094      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
1095     fhTrackMatchedDEta->SetYTitle("d#eta");
1096     fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
1097
1098     fhTrackMatchedDPhi  = new TH2F
1099     ("hTrackMatchedDPhi",
1100      "d#phi of cluster-track vs cluster energy",
1101      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
1102     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
1103     fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
1104
1105     fhTrackMatchedDEtaDPhi  = new TH2F
1106     ("hTrackMatchedDEtaDPhi",
1107      "d#eta vs d#phi of cluster-track vs cluster energy",
1108      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
1109     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
1110     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
1111
1112     outputContainer->Add(fhTrackMatchedDEta) ;
1113     outputContainer->Add(fhTrackMatchedDPhi) ;
1114     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
1115
1116     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
1117     fhdEdx->SetXTitle("E (GeV)");
1118     fhdEdx->SetYTitle("<dE/dx>");
1119     outputContainer->Add(fhdEdx);
1120
1121     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1122     fhEOverP->SetXTitle("E (GeV)");
1123     fhEOverP->SetYTitle("E/p");
1124     outputContainer->Add(fhEOverP);
1125
1126     if(fCalorimeter=="EMCAL")
1127     {
1128       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
1129       fhEOverPNoTRD->SetXTitle("E (GeV)");
1130       fhEOverPNoTRD->SetYTitle("E/p");
1131       outputContainer->Add(fhEOverPNoTRD);
1132     }
1133
1134     if(IsDataMC() && fFillTMHisto)
1135     {
1136       fhTrackMatchedMCParticleE  = new TH2F
1137       ("hTrackMatchedMCParticleE",
1138        "Origin of particle vs energy",
1139        nptbins,ptmin,ptmax,8,0,8);
1140       fhTrackMatchedMCParticleE->SetXTitle("E (GeV)");
1141       //fhTrackMatchedMCParticleE->SetYTitle("Particle type");
1142
1143       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(1 ,"Photon");
1144       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(2 ,"Electron");
1145       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1146       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(4 ,"Rest");
1147       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1148       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1149       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1150       fhTrackMatchedMCParticleE->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1151
1152       outputContainer->Add(fhTrackMatchedMCParticleE);
1153
1154       fhTrackMatchedMCParticleDEta  = new TH2F
1155       ("hTrackMatchedMCParticleDEta",
1156        "Origin of particle vs #eta residual",
1157        nresetabins,resetamin,resetamax,8,0,8);
1158       fhTrackMatchedMCParticleDEta->SetXTitle("#Delta #eta");
1159       //fhTrackMatchedMCParticleDEta->SetYTitle("Particle type");
1160
1161       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(1 ,"Photon");
1162       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(2 ,"Electron");
1163       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1164       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(4 ,"Rest");
1165       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1166       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1167       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1168       fhTrackMatchedMCParticleDEta->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1169
1170       outputContainer->Add(fhTrackMatchedMCParticleDEta);
1171
1172       fhTrackMatchedMCParticleDPhi  = new TH2F
1173       ("hTrackMatchedMCParticleDPhi",
1174        "Origin of particle vs #phi residual",
1175        nresphibins,resphimin,resphimax,8,0,8);
1176       fhTrackMatchedMCParticleDPhi->SetXTitle("#Delta #phi");
1177       //fhTrackMatchedMCParticleDPhi->SetYTitle("Particle type");
1178
1179       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(1 ,"Photon");
1180       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(2 ,"Electron");
1181       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
1182       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(4 ,"Rest");
1183       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
1184       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
1185       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
1186       fhTrackMatchedMCParticleDPhi->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
1187
1188       outputContainer->Add(fhTrackMatchedMCParticleDPhi);
1189
1190
1191     }
1192   }
1193
1194   if(fFillWeightHistograms)
1195   {
1196     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1197                                      nptbins,ptmin,ptmax, 100,0,1.);
1198     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1199     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
1200     outputContainer->Add(fhECellClusterRatio);
1201
1202     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1203                                         nptbins,ptmin,ptmax, 100,-10,0);
1204     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1205     fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1206     outputContainer->Add(fhECellClusterLogRatio);
1207
1208     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
1209                                         nptbins,ptmin,ptmax, 100,0,1.);
1210     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
1211     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
1212     outputContainer->Add(fhEMaxCellClusterRatio);
1213
1214     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
1215                                            nptbins,ptmin,ptmax, 100,-10,0);
1216     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
1217     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
1218     outputContainer->Add(fhEMaxCellClusterLogRatio);
1219
1220     for(Int_t iw = 0; iw < 14; iw++)
1221     {
1222       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected decay photons from neutral meson",1+0.5*iw),
1223                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1224       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
1225       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
1226       outputContainer->Add(fhLambda0ForW0[iw]);
1227
1228 //      fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected decay photons from neutral meson",0.5+0.5*iw),
1229 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1230 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
1231 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
1232 //      outputContainer->Add(fhLambda1ForW0[iw]);
1233
1234     }
1235   }
1236
1237   if(IsDataMC())
1238   {
1239     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC && fAnaType==kSSCalo)
1240     {
1241       fhMCPi0PtGenRecoFraction = new TH2F("hMCPi0PtGenRecoFraction","Number of clusters from #pi^{0} (2 #gamma) identified as #pi^{0} (#eta), pT versus E primary #pi^{0} / E reco",
1242                                    nptbins,ptmin,ptmax,200,0,2);
1243       fhMCPi0PtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1244       fhMCPi0PtGenRecoFraction->SetYTitle("E^{ #pi^{0} mother} / E^{rec}");
1245       outputContainer->Add(fhMCPi0PtGenRecoFraction) ;
1246
1247       fhMCEtaPtGenRecoFraction = new TH2F("hMCEtaPtGenRecoFraction","Number of clusters from #eta (2 #gamma) identified as #pi^{0} (#eta),pT versus E primary #eta / E reco",
1248                                    nptbins,ptmin,ptmax,200,0,2);
1249       fhMCEtaPtGenRecoFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1250       fhMCEtaPtGenRecoFraction->SetYTitle("E^{ #eta mother} / E^{rec}");
1251       outputContainer->Add(fhMCEtaPtGenRecoFraction) ;
1252
1253       fhMCPi0DecayPt = new TH1F("hMCPi0DecayPt","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1254       fhMCPi0DecayPt->SetYTitle("N");
1255       fhMCPi0DecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1256       outputContainer->Add(fhMCPi0DecayPt) ;
1257
1258       fhMCPi0DecayPtFraction = new TH2F("hMCPi0DecayPtFraction","Number of #gamma from #pi^{0} decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #pi^{0}",
1259                                         nptbins,ptmin,ptmax,100,0,1);
1260       fhMCPi0DecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1261       fhMCPi0DecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1262       outputContainer->Add(fhMCPi0DecayPtFraction) ;
1263
1264       fhMCEtaDecayPt = new TH1F("hMCEtaDecayPt","Number of #gamma from #eta decay  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1265       fhMCEtaDecayPt->SetYTitle("N");
1266       fhMCEtaDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1267       outputContainer->Add(fhMCEtaDecayPt) ;
1268
1269       fhMCEtaDecayPtFraction = new TH2F("hMCEtaDecayPtFraction","Number of #gamma from #eta decay  identified as #pi^{0} (#eta), pT versus E primary  #gamma / E primary #eta",
1270                                         nptbins,ptmin,ptmax,100,0,1);
1271       fhMCEtaDecayPtFraction->SetXTitle("p^{rec}_{T} (GeV/c)");
1272       fhMCEtaDecayPtFraction->SetYTitle("E^{gen} / E^{gen-mother}");
1273       outputContainer->Add(fhMCEtaDecayPtFraction) ;
1274
1275       fhMCOtherDecayPt = new TH1F("hMCOtherDecayPt","Number of #gamma decay (not #eta or #pi^{0})  identified as #pi^{0} (#eta)",nptbins,ptmin,ptmax);
1276       fhMCOtherDecayPt->SetYTitle("N");
1277       fhMCOtherDecayPt->SetXTitle("p^{rec}_{T} (GeV/c)");
1278       outputContainer->Add(fhMCOtherDecayPt) ;
1279
1280     }
1281
1282     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1283        GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1284     {
1285
1286       fhAnglePairMCPi0  = new TH2F
1287       ("AnglePairMCPi0",
1288        "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
1289       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
1290       fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
1291       outputContainer->Add(fhAnglePairMCPi0) ;
1292
1293       if (fAnaType!= kSSCalo)
1294       {
1295         fhAnglePairMCEta  = new TH2F
1296         ("AnglePairMCEta",
1297          "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1298         fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1299         fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1300         outputContainer->Add(fhAnglePairMCEta) ;
1301
1302         fhMassPairMCPi0  = new TH2F
1303         ("MassPairMCPi0",
1304          "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1305         fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1306         fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1307         outputContainer->Add(fhMassPairMCPi0) ;
1308
1309         fhMassPairMCEta  = new TH2F
1310         ("MassPairMCEta",
1311          "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1312         fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1313         fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1314         outputContainer->Add(fhMassPairMCEta) ;
1315       }
1316
1317       for(Int_t i = 0; i < 6; i++)
1318       {
1319
1320         fhMCE[i]  = new TH1F
1321         (Form("hE_MC%s",pname[i].Data()),
1322          Form("Identified as #pi^{0} (#eta), cluster from %s",
1323               ptype[i].Data()),
1324          nptbins,ptmin,ptmax);
1325         fhMCE[i]->SetYTitle("N");
1326         fhMCE[i]->SetXTitle("E (GeV)");
1327         outputContainer->Add(fhMCE[i]) ;
1328
1329         fhMCPt[i]  = new TH1F
1330         (Form("hPt_MC%s",pname[i].Data()),
1331          Form("Identified as #pi^{0} (#eta), cluster from %s",
1332               ptype[i].Data()),
1333          nptbins,ptmin,ptmax);
1334         fhMCPt[i]->SetYTitle("N");
1335         fhMCPt[i]->SetXTitle("p_{T} (GeV/c)");
1336         outputContainer->Add(fhMCPt[i]) ;
1337
1338         fhMCPtCentrality[i]  = new TH2F
1339         (Form("hPtCentrality_MC%s",pname[i].Data()),
1340          Form("Identified as #pi^{0} (#eta), cluster from %s",
1341               ptype[i].Data()),
1342          nptbins,ptmin,ptmax, 100,0,100);
1343         fhMCPtCentrality[i]->SetYTitle("centrality");
1344         fhMCPtCentrality[i]->SetXTitle("p_{T} (GeV/c)");
1345         outputContainer->Add(fhMCPtCentrality[i]) ;
1346
1347         if(fAnaType == kSSCalo)
1348         {
1349
1350           fhMCNLocMaxPt[i] = new TH2F
1351           (Form("hNLocMaxPt_MC%s",pname[i].Data()),
1352            Form("cluster from %s, pT of cluster, for NLM",ptype[i].Data()),
1353            nptbins,ptmin,ptmax,10,0,10);
1354           fhMCNLocMaxPt[i] ->SetYTitle("N maxima");
1355           fhMCNLocMaxPt[i] ->SetXTitle("p_{T} (GeV/c)");
1356           outputContainer->Add(fhMCNLocMaxPt[i]) ;
1357
1358           fhMCEReject[i]  = new TH1F
1359           (Form("hEReject_MC%s",pname[i].Data()),
1360            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1361                 ptype[i].Data()),
1362            nptbins,ptmin,ptmax);
1363           fhMCEReject[i]->SetYTitle("N");
1364           fhMCEReject[i]->SetXTitle("E (GeV)");
1365           outputContainer->Add(fhMCEReject[i]) ;
1366
1367           fhMCPtReject[i]  = new TH1F
1368           (Form("hPtReject_MC%s",pname[i].Data()),
1369            Form("Rejected as #pi^{0} (#eta), cluster from %s",
1370                 ptype[i].Data()),
1371            nptbins,ptmin,ptmax);
1372           fhMCPtReject[i]->SetYTitle("N");
1373           fhMCPtReject[i]->SetXTitle("p_{T} (GeV/c)");
1374           outputContainer->Add(fhMCPtReject[i]) ;
1375         }
1376
1377         fhMCPhi[i]  = new TH2F
1378         (Form("hPhi_MC%s",pname[i].Data()),
1379          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1380          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1381         fhMCPhi[i]->SetYTitle("#phi");
1382         fhMCPhi[i]->SetXTitle("p_{T} (GeV/c)");
1383         outputContainer->Add(fhMCPhi[i]) ;
1384
1385         fhMCEta[i]  = new TH2F
1386         (Form("hEta_MC%s",pname[i].Data()),
1387          Form("Identified as #pi^{0} (#eta), cluster from %s",
1388               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1389         fhMCEta[i]->SetYTitle("#eta");
1390         fhMCEta[i]->SetXTitle("p_{T} (GeV/c)");
1391         outputContainer->Add(fhMCEta[i]) ;
1392
1393         fhMCMassPt[i]  = new TH2F
1394         (Form("hMassPt_MC%s",pname[i].Data()),
1395          Form("all pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1396          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1397         fhMCMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1398         fhMCMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1399         outputContainer->Add(fhMCMassPt[i]) ;
1400
1401         fhMCSelectedMassPt[i]  = new TH2F
1402         (Form("hSelectedMassPt_MC%s",pname[i].Data()),
1403          Form("Selected #pi^{0} (#eta) pairs mass: p_{T} vs massfrom %s",ptype[i].Data()),
1404          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1405         fhMCSelectedMassPt[i]->SetYTitle("mass (GeV/c^{2})");
1406         fhMCSelectedMassPt[i]->SetXTitle("p_{T} (GeV/c)");
1407         outputContainer->Add(fhMCSelectedMassPt[i]) ;
1408
1409
1410         if( fFillSelectClHisto )
1411         {
1412           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1413                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1414                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1415           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1416           fhEMCLambda0[i]->SetXTitle("E (GeV)");
1417           outputContainer->Add(fhEMCLambda0[i]) ;
1418
1419           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1420                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1421                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1422           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1423           fhEMCLambda1[i]->SetXTitle("E (GeV)");
1424           outputContainer->Add(fhEMCLambda1[i]) ;
1425
1426           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1427                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1428                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1429           fhEMCDispersion[i]->SetYTitle("D^{2}");
1430           fhEMCDispersion[i]->SetXTitle("E (GeV)");
1431           outputContainer->Add(fhEMCDispersion[i]) ;
1432
1433           if(fCalorimeter=="EMCAL")
1434           {
1435             fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1436                                              Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1437                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1438             fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1439             fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1440             outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
1441
1442             if(!fFillOnlySimpleSSHisto)
1443             {
1444               fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1445                                            Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1446                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1447               fhMCEDispEta[i]->SetXTitle("E (GeV)");
1448               fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1449               outputContainer->Add(fhMCEDispEta[i]);
1450
1451               fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1452                                            Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1453                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1454               fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1455               fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1456               outputContainer->Add(fhMCEDispPhi[i]);
1457
1458               fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1459                                              Form("cluster from %s : #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
1460                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1461               fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1462               fhMCESumEtaPhi[i]->SetYTitle("#delta^{2}_{#eta #phi}");
1463               outputContainer->Add(fhMCESumEtaPhi[i]);
1464
1465               fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1466                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1467                                                   nptbins,ptmin,ptmax,200,-10,10);
1468               fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1469               fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1470               outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1471
1472               fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1473                                               Form("cluster from %s : (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptype[i].Data()),
1474                                               nptbins,ptmin,ptmax, 200,-1,1);
1475               fhMCESphericity[i]->SetXTitle("E (GeV)");
1476               fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1477               outputContainer->Add(fhMCESphericity[i]);
1478
1479               for(Int_t ie = 0; ie < 7; ie++)
1480               {
1481                 fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1482                                                       Form("cluster from %s : #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1483                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1484                 fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1485                 fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1486                 outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1487
1488                 fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1489                                                       Form("cluster from %s : #lambda^{2}_{0} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1490                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1491                 fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1492                 fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1493                 outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1494
1495                 fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1496                                                       Form("cluster from %s :#lambda^{2}_{0} vs #sigma^{2}_{#phi #phi} for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1497                                                       ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1498                 fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1499                 fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1500                 outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1501
1502               }
1503             }
1504           }
1505
1506           fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1507                                                     Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1508                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1509           fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1510           fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1511           outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1512
1513           fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1514                                           Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1515                                           nptbins,ptmin,ptmax,100,0,1);
1516           fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1517           fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1518           outputContainer->Add(fhEMCFracMaxCell[i]) ;
1519
1520         }//
1521       } // shower shape histo
1522
1523     } //Not MC reader
1524   }//Histos with MC
1525
1526   if(fAnaType==kSSCalo)
1527   {
1528     fhAsymmetry  = new TH2F ("hAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1529                               nptbins,ptmin,ptmax, 200, -1,1);
1530     fhAsymmetry->SetXTitle("E (GeV)");
1531     fhAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1532     outputContainer->Add(fhAsymmetry);
1533
1534     fhSelectedAsymmetry  = new TH2F ("hSelectedAsymmetry","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1535                              nptbins,ptmin,ptmax, 200, -1,1);
1536     fhSelectedAsymmetry->SetXTitle("E (GeV)");
1537     fhSelectedAsymmetry->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1538     outputContainer->Add(fhSelectedAsymmetry);
1539
1540     fhSplitE  = new TH1F
1541     ("hSplitE","Selected #pi^{0} (#eta) pairs energy sum of split sub-clusters",nptbins,ptmin,ptmax);
1542     fhSplitE->SetYTitle("counts");
1543     fhSplitE->SetXTitle("E (GeV)");
1544     outputContainer->Add(fhSplitE) ;
1545
1546     fhSplitPt  = new TH1F
1547     ("hSplitPt","Selected #pi^{0} (#eta) pairs pT sum of split sub-clusters",nptbins,ptmin,ptmax);
1548     fhSplitPt->SetYTitle("counts");
1549     fhSplitPt->SetXTitle("p_{T} (GeV/c)");
1550     outputContainer->Add(fhSplitPt) ;
1551
1552
1553     fhSplitPtPhi  = new TH2F
1554     ("hSplitPtPhi","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
1555     fhSplitPtPhi->SetYTitle("#phi (rad)");
1556     fhSplitPtPhi->SetXTitle("p_{T} (GeV/c)");
1557     outputContainer->Add(fhSplitPtPhi) ;
1558
1559     fhSplitPtEta  = new TH2F
1560     ("hSplitPtEta","Selected #pi^{0} (#eta) pairs: sum split sub-cluster p_{T} vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
1561     fhSplitPtEta->SetYTitle("#eta");
1562     fhSplitPtEta->SetXTitle("p_{T} (GeV/c)");
1563     outputContainer->Add(fhSplitPtEta) ;
1564
1565
1566     fhNLocMaxSplitPt = new TH2F("hNLocMaxSplitPt","Number of local maxima in cluster",
1567                          nptbins,ptmin,ptmax,10,0,10);
1568     fhNLocMaxSplitPt ->SetYTitle("N maxima");
1569     fhNLocMaxSplitPt ->SetXTitle("p_{T} (GeV/c)");
1570     outputContainer->Add(fhNLocMaxSplitPt) ;
1571
1572
1573     fhMassSplitPt  = new TH2F
1574     ("hMassSplitPt","all pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1575     fhMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1576     fhMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1577     outputContainer->Add(fhMassSplitPt) ;
1578
1579     fhSelectedMassSplitPt  = new TH2F
1580     ("hSelectedMassSplitPt","Selected #pi^{0} (#eta) pairs mass: sum split sub-cluster p_{T} vs mass",nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1581     fhSelectedMassSplitPt->SetYTitle("mass (GeV/c^{2})");
1582     fhSelectedMassSplitPt->SetXTitle("p_{T} (GeV/c)");
1583     outputContainer->Add(fhSelectedMassSplitPt) ;
1584
1585
1586
1587     if(IsDataMC())
1588     {
1589       for(Int_t i = 0; i< 6; i++)
1590       {
1591         fhMCEAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1592                                        Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1593                                        nptbins,ptmin,ptmax, 200,-1,1);
1594         fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1595         fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1596         outputContainer->Add(fhMCEAsymmetry[i]);
1597
1598         fhMCSplitE[i]  = new TH1F
1599         (Form("hSplitE_MC%s",pname[i].Data()),
1600          Form("cluster from %s, energy sum of split sub-clusters",ptype[i].Data()),
1601          nptbins,ptmin,ptmax);
1602         fhMCSplitE[i]->SetYTitle("counts");
1603         fhMCSplitE[i]->SetXTitle("E (GeV)");
1604         outputContainer->Add(fhMCSplitE[i]) ;
1605
1606         fhMCSplitPt[i]  = new TH1F
1607         (Form("hSplitPt_MC%s",pname[i].Data()),
1608          Form("cluster from %s, pT sum of split sub-clusters",ptype[i].Data()),
1609          nptbins,ptmin,ptmax);
1610         fhMCSplitPt[i]->SetYTitle("counts");
1611         fhMCSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1612         outputContainer->Add(fhMCSplitPt[i]) ;
1613
1614
1615         fhMCSplitPtPhi[i]  = new TH2F
1616         (Form("hSplitPtPhi_MC%s",pname[i].Data()),
1617          Form("Identified as #pi^{0} (#eta), cluster from %s",ptype[i].Data()),
1618          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1619         fhMCSplitPtPhi[i]->SetYTitle("#phi");
1620         fhMCSplitPtPhi[i]->SetXTitle("p_{T} (GeV/c)");
1621         outputContainer->Add(fhMCSplitPtPhi[i]) ;
1622
1623         fhMCSplitPtEta[i]  = new TH2F
1624         (Form("hSplitPtEta_MC%s",pname[i].Data()),
1625          Form("Identified as #pi^{0} (#eta), cluster from %s",
1626               ptype[i].Data()),nptbins,ptmin,ptmax,netabins,etamin,etamax);
1627         fhMCSplitPtEta[i]->SetYTitle("#eta");
1628         fhMCSplitPtEta[i]->SetXTitle("p_{T} (GeV/c)");
1629         outputContainer->Add(fhMCSplitPtEta[i]) ;
1630
1631
1632         fhMCNLocMaxSplitPt[i] = new TH2F
1633         (Form("hNLocMaxSplitPt_MC%s",pname[i].Data()),
1634          Form("cluster from %s, pT sum of split sub-clusters, for NLM",ptype[i].Data()),
1635          nptbins,ptmin,ptmax,10,0,10);
1636         fhMCNLocMaxSplitPt[i] ->SetYTitle("N maxima");
1637         fhMCNLocMaxSplitPt[i] ->SetXTitle("p_{T} (GeV/c)");
1638         outputContainer->Add(fhMCNLocMaxSplitPt[i]) ;
1639
1640         fhMCMassSplitPt[i]  = new TH2F
1641         (Form("hMassSplitPt_MC%s",pname[i].Data()),
1642          Form("all pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
1643          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1644         fhMCMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
1645         fhMCMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1646         outputContainer->Add(fhMCMassSplitPt[i]) ;
1647
1648         fhMCSelectedMassSplitPt[i]  = new TH2F
1649         (Form("hSelectedMassSplitPt_MC%s",pname[i].Data()),
1650          Form("Selected #pi^{0} (#eta) pairs mass: split p_{T} vs mass from %s",ptype[i].Data()),
1651          nptbins,ptmin,ptmax, nmassbins,massmin,massmax);
1652         fhMCSelectedMassSplitPt[i]->SetYTitle("mass (GeV/c^{2})");
1653         fhMCSelectedMassSplitPt[i]->SetXTitle("p_{T} (GeV/c)");
1654         outputContainer->Add(fhMCSelectedMassSplitPt[i]) ;
1655
1656       }
1657     }
1658   }
1659
1660   if(fAnaType==kSSCalo && fFillSelectClHisto && !fFillOnlySimpleSSHisto )
1661   {
1662
1663
1664     for(Int_t i = 0; i< 3; i++)
1665     {
1666       fhEAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1667                                         Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1668                                         nptbins,ptmin,ptmax,200, -1,1);
1669       fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1670       fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1671       outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1672     }
1673
1674     for(Int_t ie = 0; ie< 7; ie++)
1675     {
1676
1677       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1678                                          Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1679                                          ssbins,ssmin,ssmax , 200,-1,1);
1680       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1681       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1682       outputContainer->Add(fhAsymmetryLambda0[ie]);
1683
1684       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1685                                          Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1686                                          ssbins,ssmin,ssmax , 200,-1,1);
1687       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1688       fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1689       outputContainer->Add(fhAsymmetryDispEta[ie]);
1690
1691       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1692                                          Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1693                                          ssbins,ssmin,ssmax , 200,-1,1);
1694       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1695       fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1696       outputContainer->Add(fhAsymmetryDispPhi[ie]);
1697     }
1698
1699
1700     if(IsDataMC())
1701     {
1702       for(Int_t i = 0; i< 6; i++)
1703       {
1704         for(Int_t ie = 0; ie < 7; ie++)
1705         {
1706           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1707                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1708                                                   ssbins,ssmin,ssmax , 200,-1,1);
1709           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1710           fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1711           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1712
1713           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1714                                                   Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1715                                                   ssbins,ssmin,ssmax , 200,-1,1);
1716           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1717           fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1718           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1719
1720           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1721                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1722                                                   ssbins,ssmin,ssmax , 200,-1,1);
1723           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1724           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1725           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1726         }
1727       }
1728     }
1729   }
1730
1731   if(fFillPileUpHistograms)
1732   {
1733
1734     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
1735
1736     for(Int_t i = 0 ; i < 7 ; i++)
1737     {
1738       fhPtPi0PileUp[i]  = new TH1F(Form("hPtPi0PileUp%s",pileUpName[i].Data()),
1739                                       Form("Selected #pi^{0} (#eta) p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()), nptbins,ptmin,ptmax);
1740       fhPtPi0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
1741       outputContainer->Add(fhPtPi0PileUp[i]);
1742     }
1743
1744     fhTimeENoCut  = new TH2F ("hTimeE_NoCut","time of cluster vs E of clusters, no cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1745     fhTimeENoCut->SetXTitle("E (GeV)");
1746     fhTimeENoCut->SetYTitle("time (ns)");
1747     outputContainer->Add(fhTimeENoCut);
1748
1749     fhTimeESPD  = new TH2F ("hTimeE_SPD","time of cluster vs E of clusters, SPD cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1750     fhTimeESPD->SetXTitle("E (GeV)");
1751     fhTimeESPD->SetYTitle("time (ns)");
1752     outputContainer->Add(fhTimeESPD);
1753
1754     fhTimeESPDMulti  = new TH2F ("hTimeE_SPDMulti","time of cluster vs E of clusters, SPD multi cut", nptbins,ptmin,ptmax, ntimebins,timemin,timemax);
1755     fhTimeESPDMulti->SetXTitle("E (GeV)");
1756     fhTimeESPDMulti->SetYTitle("time (ns)");
1757     outputContainer->Add(fhTimeESPDMulti);
1758
1759     fhTimeNPileUpVertSPD  = new TH2F ("hTime_NPileUpVertSPD","time of cluster vs N pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1760     fhTimeNPileUpVertSPD->SetYTitle("# vertex ");
1761     fhTimeNPileUpVertSPD->SetXTitle("time (ns)");
1762     outputContainer->Add(fhTimeNPileUpVertSPD);
1763
1764     fhTimeNPileUpVertTrack  = new TH2F ("hTime_NPileUpVertTracks","time of cluster vs N pile-up Tracks vertex", ntimebins,timemin,timemax, 50,0,50 );
1765     fhTimeNPileUpVertTrack->SetYTitle("# vertex ");
1766     fhTimeNPileUpVertTrack->SetXTitle("time (ns)");
1767     outputContainer->Add(fhTimeNPileUpVertTrack);
1768
1769     fhTimeNPileUpVertContributors  = new TH2F ("hTime_NPileUpVertContributors","time of cluster vs N constributors to pile-up SPD vertex", ntimebins,timemin,timemax,50,0,50);
1770     fhTimeNPileUpVertContributors->SetYTitle("# vertex ");
1771     fhTimeNPileUpVertContributors->SetXTitle("time (ns)");
1772     outputContainer->Add(fhTimeNPileUpVertContributors);
1773
1774     fhTimePileUpMainVertexZDistance  = new TH2F ("hTime_PileUpMainVertexZDistance","time of cluster vs distance in Z pile-up SPD vertex - main SPD vertex",ntimebins,timemin,timemax,100,0,50);
1775     fhTimePileUpMainVertexZDistance->SetYTitle("distance Z (cm) ");
1776     fhTimePileUpMainVertexZDistance->SetXTitle("time (ns)");
1777     outputContainer->Add(fhTimePileUpMainVertexZDistance);
1778
1779     fhTimePileUpMainVertexZDiamond  = new TH2F ("hTime_PileUpMainVertexZDiamond","time of cluster vs distance in Z pile-up SPD vertex - z diamond",ntimebins,timemin,timemax,100,0,50);
1780     fhTimePileUpMainVertexZDiamond->SetYTitle("diamond distance Z (cm) ");
1781     fhTimePileUpMainVertexZDiamond->SetXTitle("time (ns)");
1782     outputContainer->Add(fhTimePileUpMainVertexZDiamond);
1783
1784   }
1785
1786   //Keep neutral meson selection histograms if requiered
1787   //Setting done in AliNeutralMesonSelection
1788
1789   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1790   {
1791     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
1792
1793     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1794       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
1795
1796     delete nmsHistos;
1797   }
1798
1799   return outputContainer ;
1800
1801 }
1802
1803 //_____________________________________________
1804 Int_t AliAnaPi0EbE::GetMCIndex(const Int_t tag)
1805 {
1806
1807   // Assign mc index depending on MC bit set
1808
1809   if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
1810   {
1811     return kmcPi0 ;
1812   }//pi0
1813   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
1814   {
1815     return kmcEta ;
1816   }//eta
1817   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
1818              GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
1819   {
1820     return kmcConversion ;
1821   }//conversion photon
1822   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
1823   {
1824     return kmcPhoton ;
1825   }//photon   no conversion
1826   else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
1827   {
1828     return kmcElectron ;
1829   }//electron
1830   else
1831   {
1832     return kmcHadron ;
1833   }//other particles
1834
1835 }
1836
1837 //__________________________________________________________________
1838 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
1839                                        AliAODPWG4Particle * photon2,
1840                                        Int_t & label, Int_t & tag)
1841 {
1842   // Check the labels of pare in case mother was same pi0 or eta
1843   // Set the new AOD accordingly
1844
1845   Int_t  label1 = photon1->GetLabel();
1846   Int_t  label2 = photon2->GetLabel();
1847
1848   if(label1 < 0 || label2 < 0 ) return ;
1849
1850   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader());
1851   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader());
1852   Int_t tag1 = photon1->GetTag();
1853   Int_t tag2 = photon2->GetTag();
1854
1855   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
1856   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
1857        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
1858      (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
1859       GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
1860      )
1861   {
1862
1863     //Check if pi0/eta mother is the same
1864     if(GetReader()->ReadStack())
1865     {
1866       if(label1>=0)
1867       {
1868         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
1869         label1 = mother1->GetFirstMother();
1870         //mother1 = GetMCStack()->Particle(label1);//pi0
1871       }
1872       if(label2>=0)
1873       {
1874         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
1875         label2 = mother2->GetFirstMother();
1876         //mother2 = GetMCStack()->Particle(label2);//pi0
1877       }
1878     } // STACK
1879     else if(GetReader()->ReadAODMCParticles())
1880     {//&& (input > -1)){
1881       if(label1>=0)
1882       {
1883         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label1);//photon in kine tree
1884         label1 = mother1->GetMother();
1885         //mother1 = GetMCStack()->Particle(label1);//pi0
1886       }
1887       if(label2>=0)
1888       {
1889         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles())->At(label2);//photon in kine tree
1890         label2 = mother2->GetMother();
1891         //mother2 = GetMCStack()->Particle(label2);//pi0
1892       }
1893     }// AOD
1894
1895     //printf("mother1 %d, mother2 %d\n",label1,label2);
1896     if( label1 == label2 && label1>=0 )
1897     {
1898
1899       label = label1;
1900
1901       TLorentzVector mom1 = *(photon1->Momentum());
1902       TLorentzVector mom2 = *(photon2->Momentum());
1903
1904       Double_t angle = mom2.Angle(mom1.Vect());
1905       Double_t mass  = (mom1+mom2).M();
1906       Double_t epair = (mom1+mom2).E();
1907
1908       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
1909       {
1910         fhMassPairMCPi0 ->Fill(epair,mass);
1911         fhAnglePairMCPi0->Fill(epair,angle);
1912         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
1913       }
1914       else
1915       {
1916         fhMassPairMCEta ->Fill(epair,mass);
1917         fhAnglePairMCEta->Fill(epair,angle);
1918         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
1919       }
1920
1921     } // same label
1922   } // both from eta or pi0 decay
1923
1924 }
1925
1926 //____________________________________________________________________________
1927 void AliAnaPi0EbE::Init()
1928 {
1929   //Init
1930   //Do some checks
1931   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1932     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1933     abort();
1934   }
1935   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1936     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1937     abort();
1938   }
1939
1940 }
1941
1942 //____________________________________________________________________________
1943 void AliAnaPi0EbE::InitParameters()
1944 {
1945   //Initialize the parameters of the analysis.
1946   AddToHistogramsName("AnaPi0EbE_");
1947
1948   fInputAODGammaConvName = "PhotonsCTS" ;
1949   fAnaType = kIMCalo ;
1950   fCalorimeter = "EMCAL" ;
1951   fMinDist  = 2.;
1952   fMinDist2 = 4.;
1953   fMinDist3 = 5.;
1954
1955   fNLMECutMin[0] = 10.;
1956   fNLMECutMin[1] = 6. ;
1957   fNLMECutMin[2] = 6. ;
1958
1959 }
1960
1961 //__________________________________________________________________
1962 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
1963 {
1964   //Do analysis and fill aods
1965
1966   switch(fAnaType)
1967   {
1968     case kIMCalo:
1969       MakeInvMassInCalorimeter();
1970       break;
1971
1972     case kSSCalo:
1973       MakeShowerShapeIdentification();
1974       break;
1975
1976     case kIMCaloTracks:
1977       MakeInvMassInCalorimeterAndCTS();
1978       break;
1979
1980   }
1981 }
1982
1983 //____________________________________________
1984 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
1985 {
1986   //Do analysis and fill aods
1987   //Search for the photon decay in calorimeters
1988   //Read photon list from AOD, produced in class AliAnaPhoton
1989   //Check if 2 photons have the mass of the pi0.
1990
1991   TLorentzVector mom1;
1992   TLorentzVector mom2;
1993   TLorentzVector mom ;
1994
1995   Int_t tag   = 0;
1996   Int_t label = 0;
1997
1998   if(!GetInputAODBranch()){
1999     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
2000     abort();
2001   }
2002
2003   //Get shower shape information of clusters
2004   TObjArray *clusters = 0;
2005   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2006   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2007
2008   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
2009     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2010
2011     //Vertex cut in case of mixed events
2012     Int_t evtIndex1 = 0 ;
2013     if(GetMixedEvent())
2014       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
2015     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
2016     mom1 = *(photon1->Momentum());
2017
2018     //Get original cluster, to recover some information
2019     Int_t iclus = -1;
2020     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2021
2022     if(!cluster1){
2023       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
2024       return;
2025     }
2026
2027     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
2028     {
2029       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
2030
2031       Int_t evtIndex2 = 0 ;
2032       if(GetMixedEvent())
2033         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2034
2035       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
2036         continue ;
2037
2038       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
2039
2040       mom2 = *(photon2->Momentum());
2041
2042       //Get original cluster, to recover some information
2043       Int_t iclus2;
2044       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
2045
2046       if(!cluster2)
2047       {
2048         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
2049         return;
2050       }
2051
2052       Float_t e1    = photon1->E();
2053       Float_t e2    = photon2->E();
2054
2055       //Select clusters with good time window difference
2056       Float_t tof1  = cluster1->GetTOF()*1e9;;
2057       Float_t tof2  = cluster2->GetTOF()*1e9;;
2058       Double_t t12diff = tof1-tof2;
2059       fhEPairDiffTime->Fill(e1+e2,    t12diff);
2060       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2061
2062       //Play with the MC stack if available
2063       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
2064
2065       // Check the invariant mass for different selection on the local maxima
2066       // Name of AOD method TO BE FIXED
2067       Int_t nMaxima1 = photon1->GetFiducialArea();
2068       Int_t nMaxima2 = photon2->GetFiducialArea();
2069
2070       Double_t mass  = (mom1+mom2).M();
2071       Double_t epair = (mom1+mom2).E();
2072
2073       if(nMaxima1==nMaxima2)
2074       {
2075         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
2076         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
2077         else                 fhMassPairLocMax[2]->Fill(epair,mass);
2078       }
2079       else if(nMaxima1==1 || nMaxima2==1)
2080       {
2081         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
2082         else                             fhMassPairLocMax[4]->Fill(epair,mass);
2083       }
2084       else
2085         fhMassPairLocMax[5]->Fill(epair,mass);
2086
2087       // combinations with SS axis cut and NLM cut
2088       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2089       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
2090       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2091       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
2092
2093       //Skip events with too few or too many  NLM
2094       if((nMaxima1 < fNLMCutMin || nMaxima1 > fNLMCutMax) || (nMaxima2 < fNLMCutMin || nMaxima2 > fNLMCutMax)) continue ;
2095
2096       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - NLM of out of range: cluster1 %d, cluster2 %d \n",nMaxima1, nMaxima2);
2097
2098       //Mass of all pairs
2099       fhMass->Fill(epair,(mom1+mom2).M());
2100
2101       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2102       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2103       {
2104         if(GetDebug()>1)
2105           printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Selected gamma pair: pt %f, phi %f, eta%f \n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2106
2107         //Fill some histograms about shower shape
2108         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2109         {
2110           FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
2111           FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
2112         }
2113
2114         // Tag both photons as decay
2115         photon1->SetTagged(kTRUE);
2116         photon2->SetTagged(kTRUE);
2117
2118         fhPtDecay->Fill(photon1->Pt());
2119         fhEDecay ->Fill(photon1->E() );
2120
2121         fhPtDecay->Fill(photon2->Pt());
2122         fhEDecay ->Fill(photon2->E() );
2123
2124         //Create AOD for analysis
2125         mom = mom1+mom2;
2126
2127         //Mass of selected pairs
2128         fhSelectedMass->Fill(epair,mom.M());
2129
2130         // Fill histograms to undertand pile-up before other cuts applied
2131         // Remember to relax time cuts in the reader
2132         FillPileUpHistograms(mom.E(),((cluster1->GetTOF()+cluster2->GetTOF())*1e9) /2);
2133
2134         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2135
2136         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2137         pi0.SetDetector(photon1->GetDetector());
2138
2139         // MC
2140         pi0.SetLabel(label);
2141         pi0.SetTag(tag);
2142
2143         //Set the indeces of the original caloclusters
2144         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
2145         //pi0.SetInputFileIndex(input);
2146
2147         AddAODParticle(pi0);
2148
2149       }//pi0
2150
2151     }//2n photon loop
2152
2153   }//1st photon loop
2154
2155   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
2156
2157 }
2158
2159 //__________________________________________________
2160 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
2161 {
2162   //Do analysis and fill aods
2163   //Search for the photon decay in calorimeters
2164   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
2165   //Check if 2 photons have the mass of the pi0.
2166
2167   TLorentzVector mom1;
2168   TLorentzVector mom2;
2169   TLorentzVector mom ;
2170   Int_t tag   = 0;
2171   Int_t label = 0;
2172   Int_t evtIndex = 0;
2173
2174   // Check calorimeter input
2175   if(!GetInputAODBranch()){
2176     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
2177     abort();
2178   }
2179
2180   // Get the array with conversion photons
2181   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
2182   if(!inputAODGammaConv) {
2183
2184     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
2185
2186     if(!inputAODGammaConv) {
2187       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
2188
2189       return;
2190     }
2191   }
2192
2193   //Get shower shape information of clusters
2194   TObjArray *clusters = 0;
2195   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
2196   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
2197
2198   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
2199   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
2200   if(nCTS<=0 || nCalo <=0)
2201   {
2202     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
2203     return;
2204   }
2205
2206   if(GetDebug() > 1)
2207     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
2208
2209   // Do the loop, first calo, second CTS
2210   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
2211     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
2212     mom1 = *(photon1->Momentum());
2213
2214     //Get original cluster, to recover some information
2215     Int_t iclus = -1;
2216     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
2217
2218     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
2219       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
2220       if(GetMixedEvent())
2221         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
2222       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2223
2224       mom2 = *(photon2->Momentum());
2225
2226       Double_t mass  = (mom1+mom2).M();
2227       Double_t epair = (mom1+mom2).E();
2228
2229       Int_t nMaxima = photon1->GetFiducialArea();
2230       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
2231       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
2232       else                fhMassPairLocMax[2]->Fill(epair,mass);
2233
2234       if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) continue ;
2235       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - NLM %d of out of range \n",nMaxima);
2236
2237       //Play with the MC stack if available
2238       if(IsDataMC())
2239       {
2240         Int_t   label2 = photon2->GetLabel();
2241         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader()));
2242
2243         HasPairSameMCMother(photon1, photon2, label, tag) ;
2244       }
2245
2246       //Mass of selected pairs
2247       fhMass->Fill(epair,(mom1+mom2).M());
2248
2249       //Select good pair (good phi, pt cuts, aperture and invariant mass)
2250       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
2251       {
2252         if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Selected gamma pair: pt %f, phi %f, eta%f\n",(mom1+mom2).Pt(), (mom1+mom2).Phi()*180./3.1416, (mom1+mom2).Eta());
2253
2254         //Fill some histograms about shower shape
2255         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2256         {
2257           FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
2258         }
2259
2260         // Tag both photons as decay
2261         photon1->SetTagged(kTRUE);
2262         photon2->SetTagged(kTRUE);
2263
2264         fhPtDecay->Fill(photon1->Pt());
2265         fhEDecay ->Fill(photon1->E() );
2266
2267         //Create AOD for analysis
2268
2269         mom = mom1+mom2;
2270
2271         //Mass of selected pairs
2272         fhSelectedMass->Fill(epair,mom.M());
2273
2274         // Fill histograms to undertand pile-up before other cuts applied
2275         // Remember to relax time cuts in the reader
2276         if(cluster)FillPileUpHistograms(mom.E(),cluster->GetTOF()*1e9);
2277
2278         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
2279
2280         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
2281         pi0.SetDetector(photon1->GetDetector());
2282
2283         // MC
2284         pi0.SetLabel(label);
2285         pi0.SetTag(tag);
2286
2287         //Set the indeces of the original tracks or caloclusters
2288         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
2289         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
2290         //pi0.SetInputFileIndex(input);
2291
2292         AddAODParticle(pi0);
2293
2294       }//pi0
2295     }//2n photon loop
2296
2297   }//1st photon loop
2298
2299   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
2300
2301 }
2302
2303
2304 //_________________________________________________
2305 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
2306 {
2307   //Search for pi0 in fCalorimeter with shower shape analysis
2308
2309   TObjArray * pl        = 0x0;
2310   AliVCaloCells * cells = 0x0;
2311   //Select the Calorimeter of the photon
2312   if      (fCalorimeter == "PHOS" )
2313   {
2314     pl    = GetPHOSClusters();
2315     cells = GetPHOSCells();
2316   }
2317   else if (fCalorimeter == "EMCAL")
2318   {
2319     pl    = GetEMCALClusters();
2320     cells = GetEMCALCells();
2321   }
2322
2323   if(!pl)
2324   {
2325     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
2326     return;
2327   }
2328
2329   TLorentzVector mom ;
2330   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
2331   {
2332     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
2333
2334     Int_t evtIndex = 0 ;
2335     if (GetMixedEvent())
2336     {
2337       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2338     }
2339
2340     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
2341
2342     //Get Momentum vector,
2343     Double_t vertex[]={0,0,0};
2344     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
2345     {
2346       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
2347     }//Assume that come from vertex in straight line
2348     else
2349     {
2350       calo->GetMomentum(mom,vertex) ;
2351     }
2352
2353     //If too small or big pt, skip it
2354     if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) continue ;
2355
2356     //Check acceptance selection
2357     if(IsFiducialCutOn())
2358     {
2359       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
2360       if(! in ) continue ;
2361     }
2362
2363     if(GetDebug() > 1)
2364       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Min pt cut and fiducial cut passed: pt %3.2f, phi %2.2f, eta %1.2f\n",mom.Pt(),mom.Phi(),mom.Eta());
2365
2366     //Check Distance to Bad channel, set bit.
2367     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2368     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
2369     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
2370       continue ;
2371
2372     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
2373
2374     //.......................................
2375     // TOF cut, BE CAREFUL WITH THIS CUT
2376     Double_t tof = calo->GetTOF()*1e9;
2377     if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
2378
2379     //Play with the MC stack if available
2380     //Check origin of the candidates
2381     Int_t tag   = 0 ;
2382     if(IsDataMC())
2383     {
2384       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
2385       //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
2386       if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",tag);
2387     }
2388
2389     //Skip matched clusters with tracks
2390     if(fRejectTrackMatch && IsTrackMatched(calo, GetReader()->GetInputEvent()))
2391     {
2392       FillRejectedClusterHistograms(mom,tag);
2393       continue ;
2394     }
2395
2396     //Check PID
2397     //PID selection or bit setting
2398     Int_t    nMaxima = 0 ;
2399     Double_t mass    = 0 , angle = 0;
2400     TLorentzVector    l1, l2;
2401     Int_t    absId1 = -1; Int_t absId2 = -1;
2402
2403     Int_t idPartType = GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
2404                                                                                    GetVertex(evtIndex),nMaxima,
2405                                                                                    mass,angle,l1,l2,absId1,absId2) ;
2406
2407     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",idPartType);
2408
2409
2410     //Skip events with too few or too many  NLM
2411     if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax)
2412     {
2413       FillRejectedClusterHistograms(mom,tag);
2414       continue ;
2415     }
2416
2417     if(nMaxima == 1 && fNLMECutMin[0] > mom.E()) continue;
2418     if(nMaxima == 2 && fNLMECutMin[1] > mom.E()) continue;
2419     if(nMaxima >  2 && fNLMECutMin[2] > mom.E()) continue;
2420
2421     if(GetDebug() > 1)
2422       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - NLM %d accepted \n",nMaxima);
2423
2424     Float_t e1 = l1.Energy();
2425     Float_t e2 = l2.Energy();
2426     TLorentzVector l12 = l1+l2;
2427     Float_t ptSplit = l12.Pt();
2428     Float_t  eSplit = e1+e2;
2429     Int_t mcIndex = GetMCIndex(tag);
2430
2431     //mass of all clusters
2432     fhMass       ->Fill(mom.E(),mass);
2433     fhMassPt     ->Fill(mom.Pt(),mass);
2434     fhMassSplitPt->Fill(ptSplit,mass);
2435
2436     if(IsDataMC())
2437     {
2438       fhMCMassPt[mcIndex]     ->Fill(mom.Pt(),mass);
2439       fhMCMassSplitPt[mcIndex]->Fill(ptSplit,mass);
2440     }
2441
2442     // Asymmetry of all clusters
2443     Float_t asy =-10;
2444
2445     if(e1+e2 > 0) asy = (e1-e2) / (e1+e2);
2446     fhAsymmetry->Fill(mom.E(),asy);
2447
2448
2449     if(IsDataMC())
2450     {
2451       fhMCEAsymmetry[mcIndex]->Fill(mom.E(),asy);
2452     }
2453
2454     // If cluster does not pass pid, not pi0/eta, skip it.
2455     if     (GetOutputAODName().Contains("Pi0") && idPartType != AliCaloPID::kPi0)
2456     {
2457       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Pi0\n");
2458       FillRejectedClusterHistograms(mom,tag);
2459       continue ;
2460     }
2461
2462     else if(GetOutputAODName().Contains("Eta") && idPartType != AliCaloPID::kEta)
2463     {
2464       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Cluster is not Eta\n");
2465       FillRejectedClusterHistograms(mom,tag);
2466       continue ;
2467     }
2468
2469     if(GetDebug() > 1)
2470       printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0/Eta selection cuts passed: pT %3.2f, pdg %d\n",
2471                               mom.Pt(), idPartType);
2472
2473     //Mass and asymmetry of selected pairs
2474     fhSelectedAsymmetry  ->Fill(mom.E() ,asy );
2475     fhSelectedMass       ->Fill(mom.E() ,mass);
2476     fhSelectedMassPt     ->Fill(mom.Pt(),mass);
2477     fhSelectedMassSplitPt->Fill(ptSplit ,mass);
2478
2479     fhSplitE        ->Fill( eSplit);
2480     fhSplitPt       ->Fill(ptSplit);
2481     Float_t phi = mom.Phi();
2482     if(phi<0) phi+=TMath::TwoPi();
2483     fhSplitPtPhi    ->Fill(ptSplit,phi);
2484     fhSplitPtEta    ->Fill(ptSplit,mom.Eta());
2485     fhNLocMaxSplitPt->Fill(ptSplit ,nMaxima);
2486     fhNLocMaxPt     ->Fill(mom.Pt(),nMaxima);
2487
2488     //Check split-clusters with good time window difference
2489     Double_t tof1  = cells->GetCellTime(absId1);
2490     GetCaloUtils()->RecalibrateCellTime(tof1, fCalorimeter, absId1,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2491     tof1*=1.e9;
2492
2493     Double_t tof2  = cells->GetCellTime(absId2);
2494     GetCaloUtils()->RecalibrateCellTime(tof2, fCalorimeter, absId2,GetReader()->GetInputEvent()->GetBunchCrossNumber());
2495     tof2*=1.e9;
2496
2497     Double_t t12diff = tof1-tof2;
2498     fhEPairDiffTime->Fill(e1+e2,    t12diff);
2499
2500     if(IsDataMC())
2501     {
2502       fhMCSplitE        [mcIndex]->Fill( eSplit);
2503       fhMCSplitPt       [mcIndex]->Fill(ptSplit);
2504       fhMCSplitPtPhi    [mcIndex]->Fill(ptSplit,phi);
2505       fhMCSplitPtEta    [mcIndex]->Fill(ptSplit,mom.Eta());
2506       fhMCNLocMaxSplitPt[mcIndex]->Fill(ptSplit ,nMaxima);
2507       fhMCNLocMaxPt     [mcIndex]->Fill(mom.Pt(),nMaxima);
2508
2509       fhMCSelectedMassPt     [mcIndex]->Fill(mom.Pt(),mass);
2510       fhMCSelectedMassSplitPt[mcIndex]->Fill(ptSplit,mass);
2511
2512     }
2513
2514     //-----------------------
2515     //Create AOD for analysis
2516
2517     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
2518     aodpi0.SetLabel(calo->GetLabel());
2519
2520     //Set the indeces of the original caloclusters
2521     aodpi0.SetCaloLabel(calo->GetID(),-1);
2522     aodpi0.SetDetector(fCalorimeter);
2523
2524     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
2525     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
2526     else                         aodpi0.SetDistToBad(0) ;
2527
2528     // Check if cluster is pi0 via cluster splitting
2529     aodpi0.SetIdentifiedParticleType(idPartType);
2530
2531     // Add number of local maxima to AOD, method name in AOD to be FIXED
2532     aodpi0.SetFiducialArea(nMaxima);
2533
2534     aodpi0.SetTag(tag);
2535
2536     //Fill some histograms about shower shape
2537     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
2538     {
2539       FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
2540     }
2541
2542     // Fill histograms to undertand pile-up before other cuts applied
2543     // Remember to relax time cuts in the reader
2544     Double_t tofcluster   = calo->GetTOF()*1e9;
2545     Double_t tofclusterUS = TMath::Abs(tofcluster);
2546
2547     FillPileUpHistograms(calo->E(),tofcluster);
2548
2549     Int_t id = GetReader()->GetTriggerClusterId();
2550     if(fFillEMCALBCHistograms && fCalorimeter=="EMCAL" && id >=0 )
2551     {
2552       Float_t phicluster = aodpi0.Phi();
2553       if(phicluster < 0) phicluster+=TMath::TwoPi();
2554
2555       if(calo->E() > 2)
2556       {
2557         if      (tofclusterUS < 25) fhEtaPhiEMCALBC0->Fill(aodpi0.Eta(), phicluster);
2558         else if (tofclusterUS < 75) fhEtaPhiEMCALBC1->Fill(aodpi0.Eta(), phicluster);
2559         else                        fhEtaPhiEMCALBCN->Fill(aodpi0.Eta(), phicluster);
2560       }
2561
2562       Int_t bc = GetReader()->GetTriggerClusterBC();
2563       if(TMath::Abs(bc) < 6  && !GetReader()->IsBadCellTriggerEvent() && !GetReader()->IsExoticEvent() )
2564       {
2565         if(GetReader()->IsTriggerMatched())
2566         {
2567           if(calo->E() > 2) fhEtaPhiTriggerEMCALBC[bc+5]->Fill(aodpi0.Eta(), phicluster);
2568           fhTimeTriggerEMCALBC[bc+5]->Fill(calo->E(), tofcluster);
2569           if(GetReader()->IsPileUpFromSPD()) fhTimeTriggerEMCALBCPileUpSPD[bc+5]->Fill(calo->E(), tofcluster);
2570         }
2571         else
2572         {
2573           if(calo->E() > 2) fhEtaPhiTriggerEMCALBCUM[bc+5]->Fill(aodpi0.Eta(), phicluster);
2574           fhTimeTriggerEMCALBCUM[bc+5]->Fill(calo->E(), tofcluster);
2575         }
2576       }
2577       else if(TMath::Abs(bc) >= 6)
2578         printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Trigger BC not expected = %d\n",bc);
2579     }
2580
2581     //Add AOD with pi0 object to aod branch
2582     AddAODParticle(aodpi0);
2583
2584   }//loop
2585
2586   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
2587
2588 }
2589 //______________________________________________
2590 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
2591 {
2592   //Do analysis and fill histograms
2593
2594   if(!GetOutputAODBranch())
2595   {
2596     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
2597     abort();
2598   }
2599   //Loop on stored AOD pi0
2600   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2601   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2602
2603   Float_t cen = GetEventCentrality();
2604   Float_t ep  = GetEventPlaneAngle();
2605
2606   for(Int_t iaod = 0; iaod < naod ; iaod++)
2607   {
2608
2609     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2610     Int_t pdg = pi0->GetIdentifiedParticleType();
2611
2612     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
2613
2614     //Fill pi0 histograms
2615     Float_t ener  = pi0->E();
2616     Float_t pt    = pi0->Pt();
2617     Float_t phi   = pi0->Phi();
2618     if(phi < 0) phi+=TMath::TwoPi();
2619     Float_t eta = pi0->Eta();
2620
2621     fhPt     ->Fill(pt  );
2622     fhE      ->Fill(ener);
2623
2624     fhEEta   ->Fill(ener,eta);
2625     fhEPhi   ->Fill(ener,phi);
2626     fhPtEta  ->Fill(pt  ,eta);
2627     fhPtPhi  ->Fill(pt  ,phi);
2628     fhEtaPhi ->Fill(eta ,phi);
2629
2630     fhPtCentrality ->Fill(pt,cen) ;
2631     fhPtEventPlane ->Fill(pt,ep ) ;
2632
2633     if(fFillPileUpHistograms)
2634     {
2635       if(GetReader()->IsPileUpFromSPD())               fhPtPi0PileUp[0]->Fill(pt);
2636       if(GetReader()->IsPileUpFromEMCal())             fhPtPi0PileUp[1]->Fill(pt);
2637       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtPi0PileUp[2]->Fill(pt);
2638       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtPi0PileUp[3]->Fill(pt);
2639       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtPi0PileUp[4]->Fill(pt);
2640       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtPi0PileUp[5]->Fill(pt);
2641       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtPi0PileUp[6]->Fill(pt);
2642     }
2643
2644
2645     if(IsDataMC())
2646     {
2647       Int_t tag     = pi0->GetTag();
2648       Int_t mcIndex = GetMCIndex(tag);
2649
2650       fhMCE  [mcIndex] ->Fill(ener);
2651       fhMCPt [mcIndex] ->Fill(pt);
2652       fhMCPhi[mcIndex] ->Fill(pt,phi);
2653       fhMCEta[mcIndex] ->Fill(pt,eta);
2654
2655       fhMCPtCentrality[mcIndex]->Fill(pt,cen);
2656
2657       if((mcIndex==kmcPhoton || mcIndex==kmcPi0 || mcIndex==kmcEta) && fAnaType==kSSCalo)
2658       {
2659         Float_t efracMC = 0;
2660         Int_t label = pi0->GetLabel();
2661
2662         Bool_t ok = kFALSE;
2663         TLorentzVector mom   = GetMCAnalysisUtils()->GetMother(label,GetReader(),ok);
2664         if(!ok) continue;
2665
2666         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2667         {
2668           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2669           if(grandmom.E() > 0 && ok)
2670           {
2671             efracMC =  grandmom.E()/ener;
2672             fhMCPi0PtGenRecoFraction ->Fill(pt,efracMC);
2673           }
2674         }
2675         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
2676         {
2677           fhMCPi0DecayPt->Fill(pt);
2678           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,111,GetReader(),ok);
2679           if(grandmom.E() > 0 && ok)
2680           {
2681             efracMC =  mom.E()/grandmom.E();
2682             fhMCPi0DecayPtFraction ->Fill(pt,efracMC);
2683           }
2684         }
2685         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta))
2686         {
2687           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2688           if(grandmom.E() > 0 && ok)
2689           {
2690             efracMC =  grandmom.E()/ener;
2691             fhMCEtaPtGenRecoFraction ->Fill(pt,efracMC);
2692           }
2693         }
2694         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay))
2695         {
2696           fhMCEtaDecayPt->Fill(pt);
2697           TLorentzVector grandmom = GetMCAnalysisUtils()->GetMotherWithPDG(label,221,GetReader(),ok);
2698           if(grandmom.E() > 0 && ok)
2699           {
2700             efracMC =  mom.E()/grandmom.E();
2701             fhMCEtaDecayPtFraction ->Fill(pt,efracMC);
2702           }
2703         }
2704         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2705         {
2706           fhMCOtherDecayPt->Fill(pt);
2707         }
2708
2709       }
2710
2711     }//Histograms with MC
2712
2713   }// aod loop
2714
2715 }
2716
2717 //__________________________________________________________________
2718 void AliAnaPi0EbE::Print(const Option_t * opt) const
2719 {
2720   //Print some relevant parameters set for the analysis
2721   if(! opt)
2722     return;
2723
2724   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2725   AliAnaCaloTrackCorrBaseClass::Print("");
2726   printf("Analysis Type = %d \n",  fAnaType) ;
2727   if(fAnaType == kSSCalo){
2728     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2729     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2730     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2731     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2732   }
2733   printf("    \n") ;
2734
2735 }
2736
2737