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