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