PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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