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