PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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