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