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("#sigma'^{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("#sigma'^{2}_{#phi #phi}");
 729       outputContainer->Add(fhSumPhiE);
 730
 731       fhSumEtaPhiE  = new TH2F ("hSumEtaPhiE","#sigma'^{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("#sigma'^{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}_{#phi #phi}");
 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       fhAnglePairMCEta  = new TH2F
1040       ("AnglePairMCEta",
1041        "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
1042       fhAnglePairMCEta->SetYTitle("#alpha (rad)");
1043       fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
1044       outputContainer->Add(fhAnglePairMCEta) ;
1045
1046       fhMassPairMCPi0  = new TH2F
1047       ("MassPairMCPi0",
1048        "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1049       fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
1050       fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
1051       outputContainer->Add(fhMassPairMCPi0) ;
1052
1053       fhMassPairMCEta  = new TH2F
1054       ("MassPairMCEta",
1055        "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
1056       fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
1057       fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
1058       outputContainer->Add(fhMassPairMCEta) ;
1059
1060       if( fFillSelectClHisto )
1061       {
1062         for(Int_t i = 0; i < 6; i++)
1063         {
1064           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
1065                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
1066                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1067           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
1068           fhEMCLambda0[i]->SetXTitle("E (GeV)");
1069           outputContainer->Add(fhEMCLambda0[i]) ;
1070
1071           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
1072                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
1073                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1074           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
1075           fhEMCLambda1[i]->SetXTitle("E (GeV)");
1076           outputContainer->Add(fhEMCLambda1[i]) ;
1077
1078           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
1079                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
1080                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1081           fhEMCDispersion[i]->SetYTitle("D^{2}");
1082           fhEMCDispersion[i]->SetXTitle("E (GeV)");
1083           outputContainer->Add(fhEMCDispersion[i]) ;
1084
1085           if(fCalorimeter=="EMCAL")
1086           {
1087             fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
1088                                              Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
1089                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1090             fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
1091             fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
1092             outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
1093
1094
1095             fhMCEDispEta[i]  = new TH2F (Form("hEDispEtaE_MC%s",pname[i].Data()),
1096                                          Form("cluster from %s : #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptype[i].Data()),
1097                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1098             fhMCEDispEta[i]->SetXTitle("E (GeV)");
1099             fhMCEDispEta[i]->SetYTitle("#sigma^{2}_{#eta #eta}");
1100             outputContainer->Add(fhMCEDispEta[i]);
1101
1102             fhMCEDispPhi[i]  = new TH2F (Form("hEDispPhiE_MC%s",pname[i].Data()),
1103                                          Form("cluster from %s : #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptype[i].Data()),
1104                                          nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
1105             fhMCEDispPhi[i]->SetXTitle("E (GeV)");
1106             fhMCEDispPhi[i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1107             outputContainer->Add(fhMCEDispPhi[i]);
1108
1109             fhMCESumEtaPhi[i]  = new TH2F (Form("hESumEtaPhiE_MC%s",pname[i].Data()),
1110                                            Form("cluster from %s : #sigma'^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptype[i].Data()),
1111                                            nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
1112             fhMCESumEtaPhi[i]->SetXTitle("E (GeV)");
1113             fhMCESumEtaPhi[i]->SetYTitle("#sigma'^{2}_{#eta #phi}");
1114             outputContainer->Add(fhMCESumEtaPhi[i]);
1115
1116             fhMCEDispEtaPhiDiff[i]  = new TH2F (Form("hEDispEtaPhiDiffE_MC%s",pname[i].Data()),
1117                                                 Form("cluster from %s : #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptype[i].Data()),
1118                                                 nptbins,ptmin,ptmax,200,-10,10);
1119             fhMCEDispEtaPhiDiff[i]->SetXTitle("E (GeV)");
1120             fhMCEDispEtaPhiDiff[i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
1121             outputContainer->Add(fhMCEDispEtaPhiDiff[i]);
1122
1123             fhMCESphericity[i]  = new TH2F (Form("hESphericity_MC%s",pname[i].Data()),
1124                                             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()),
1125                                             nptbins,ptmin,ptmax, 200,-1,1);
1126             fhMCESphericity[i]->SetXTitle("E (GeV)");
1127             fhMCESphericity[i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
1128             outputContainer->Add(fhMCESphericity[i]);
1129
1130             for(Int_t ie = 0; ie < 7; ie++)
1131             {
1132               fhMCDispEtaDispPhi[ie][i] = new TH2F (Form("hMCDispEtaDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1133                                                     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]),
1134                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1135               fhMCDispEtaDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1136               fhMCDispEtaDispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1137               outputContainer->Add(fhMCDispEtaDispPhi[ie][i]);
1138
1139               fhMCLambda0DispEta[ie][i] = new TH2F (Form("hMCLambda0DispEta_EBin%d_MC%s",ie,pname[i].Data()),
1140                                                     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]),
1141                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1142               fhMCLambda0DispEta[ie][i]->SetXTitle("#lambda^{2}_{0}");
1143               fhMCLambda0DispEta[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1144               outputContainer->Add(fhMCLambda0DispEta[ie][i]);
1145
1146               fhMCLambda0DispPhi[ie][i] = new TH2F (Form("hMCLambda0DispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1147                                                     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]),
1148                                                     ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
1149               fhMCLambda0DispPhi[ie][i]->SetXTitle("#lambda^{2}_{0}");
1150               fhMCLambda0DispPhi[ie][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
1151               outputContainer->Add(fhMCLambda0DispPhi[ie][i]);
1152
1153             }
1154           }
1155
1156           fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
1157                                                     Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
1158                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1159           fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
1160           fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
1161           outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
1162
1163           fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
1164                                           Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
1165                                           nptbins,ptmin,ptmax,100,0,1);
1166           fhEMCFracMaxCell[i]->SetYTitle("Fraction");
1167           fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
1168           outputContainer->Add(fhEMCFracMaxCell[i]) ;
1169
1170         }//
1171       } // shower shape histo
1172
1173     } //Not MC reader
1174   }//Histos with MC
1175
1176
1177   if(fAnaType==kSSCalo && fFillSelectClHisto )
1178   {
1179
1180     fhAsymmetryE  = new TH2F ("hAsymmetryE","A = ( E1 - E2 ) / ( E1 + E2 ) vs E",
1181                                nptbins,ptmin,ptmax, 200, -1,1);
1182     fhAsymmetryE->SetXTitle("E (GeV)");
1183     fhAsymmetryE->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1184     outputContainer->Add(fhAsymmetryE);
1185
1186     for(Int_t i = 0; i< 3; i++)
1187     {
1188       fhEAsymmetryLocMax[i]  = new TH2F(Form("hEAsymmetryLocMax%d",i+1),
1189                                         Form("Selected #pi^{0} (#eta) pairs: E vs A = ( E1 - E2 ) / ( E1 + E2 ), %s",nlm[i].Data()),
1190                                         nptbins,ptmin,ptmax,200, -1,1);
1191       fhEAsymmetryLocMax[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1192       fhEAsymmetryLocMax[i]->SetXTitle("E (GeV)");
1193       outputContainer->Add(fhEAsymmetryLocMax[i]) ;
1194     }
1195
1196     for(Int_t ie = 0; ie< 7; ie++)
1197     {
1198
1199       fhAsymmetryLambda0[ie] = new TH2F (Form("hAsymmetryLambda0_EBin%d",ie),
1200                                          Form("#lambda_{0}^{2} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1201                                          ssbins,ssmin,ssmax , 200,-1,1);
1202       fhAsymmetryLambda0[ie]->SetXTitle("#lambda_{0}^{2}");
1203       fhAsymmetryLambda0[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1204       outputContainer->Add(fhAsymmetryLambda0[ie]);
1205
1206       fhAsymmetryDispEta[ie] = new TH2F (Form("hAsymmetryDispEta_EBin%d",ie),
1207                                          Form("#sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1208                                          ssbins,ssmin,ssmax , 200,-1,1);
1209       fhAsymmetryDispEta[ie]->SetXTitle("#sigma^{2}_{#eta #eta}");
1210       fhAsymmetryDispEta[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1211       outputContainer->Add(fhAsymmetryDispEta[ie]);
1212
1213       fhAsymmetryDispPhi[ie] = new TH2F (Form("hAsymmetryDispPhi_EBin%d",ie),
1214                                          Form("#sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",bin[ie],bin[ie+1]),
1215                                          ssbins,ssmin,ssmax , 200,-1,1);
1216       fhAsymmetryDispPhi[ie]->SetXTitle("#sigma^{2}_{#phi #phi}");
1217       fhAsymmetryDispPhi[ie]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1218       outputContainer->Add(fhAsymmetryDispPhi[ie]);
1219     }
1220
1221
1222     if(IsDataMC())
1223     {
1224       for(Int_t i = 0; i< 6; i++)
1225       {
1226         fhMCEAsymmetry[i]  = new TH2F (Form("hEAsymmetry_MC%s",pname[i].Data()),
1227                                        Form("cluster from %s : A = ( E1 - E2 ) / ( E1 + E2 ) vs E",ptype[i].Data()),
1228                                        nptbins,ptmin,ptmax, 200,-1,1);
1229         fhMCEAsymmetry[i]->SetXTitle("E (GeV)");
1230         fhMCEAsymmetry[i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1231         outputContainer->Add(fhMCEAsymmetry[i]);
1232
1233         for(Int_t ie = 0; ie < 7; ie++)
1234         {
1235           fhMCAsymmetryLambda0[ie][i] = new TH2F (Form("hMCAsymmetryLambda0_EBin%d_MC%s",ie,pname[i].Data()),
1236                                                   Form("cluster from %s : #lambda_{0}^{2} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1237                                                   ssbins,ssmin,ssmax , 200,-1,1);
1238           fhMCAsymmetryLambda0[ie][i]->SetXTitle("#lambda_{0}^{2}");
1239           fhMCAsymmetryLambda0[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1240           outputContainer->Add(fhMCAsymmetryLambda0[ie][i]);
1241
1242           fhMCAsymmetryDispEta[ie][i] = new TH2F (Form("hMCAsymmetryDispEta_EBin%d_MC%s",ie,pname[i].Data()),
1243                                                   Form("cluster from %s : #sigma^{2}_{#eta #eta} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1244                                                   ssbins,ssmin,ssmax , 200,-1,1);
1245           fhMCAsymmetryDispEta[ie][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
1246           fhMCAsymmetryDispEta[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1247           outputContainer->Add(fhMCAsymmetryDispEta[ie][i]);
1248
1249           fhMCAsymmetryDispPhi[ie][i] = new TH2F (Form("hMCAsymmetryDispPhi_EBin%d_MC%s",ie,pname[i].Data()),
1250                                                   Form("cluster from %s : #sigma^{2}_{#phi #phi} vs A for %d < E < %d GeV",pname[i].Data(),bin[ie],bin[ie+1]),
1251                                                   ssbins,ssmin,ssmax , 200,-1,1);
1252           fhMCAsymmetryDispPhi[ie][i]->SetXTitle("#sigma^{2}_{#phi #phi}");
1253           fhMCAsymmetryDispPhi[ie][i]->SetYTitle("A = ( E1 - E2 ) / ( E1 + E2 )");
1254           outputContainer->Add(fhMCAsymmetryDispPhi[ie][i]);
1255         }
1256       }
1257     }
1258   }
1259
1260   //Keep neutral meson selection histograms if requiered
1261   //Setting done in AliNeutralMesonSelection
1262
1263   if(fAnaType!=kSSCalo && GetNeutralMesonSelection())
1264   {
1265     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
1266
1267     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
1268       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
1269
1270     delete nmsHistos;
1271   }
1272
1273   return outputContainer ;
1274
1275 }
1276
1277 //____________________________________________________________________________
1278 void AliAnaPi0EbE::Init()
1279 {
1280   //Init
1281   //Do some checks
1282   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1283     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1284     abort();
1285   }
1286   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1287     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1288     abort();
1289   }
1290
1291 }
1292
1293 //____________________________________________________________________________
1294 void AliAnaPi0EbE::InitParameters()
1295 {
1296   //Initialize the parameters of the analysis.
1297   AddToHistogramsName("AnaPi0EbE_");
1298
1299   fInputAODGammaConvName = "PhotonsCTS" ;
1300   fAnaType = kIMCalo ;
1301   fCalorimeter = "EMCAL" ;
1302   fMinDist  = 2.;
1303   fMinDist2 = 4.;
1304   fMinDist3 = 5.;
1305
1306 }
1307
1308 //__________________________________________________________________
1309 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
1310 {
1311   //Do analysis and fill aods
1312
1313   switch(fAnaType)
1314   {
1315     case kIMCalo:
1316       MakeInvMassInCalorimeter();
1317       break;
1318
1319     case kSSCalo:
1320       MakeShowerShapeIdentification();
1321       break;
1322
1323     case kIMCaloTracks:
1324       MakeInvMassInCalorimeterAndCTS();
1325       break;
1326
1327   }
1328 }
1329
1330 //____________________________________________
1331 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
1332 {
1333   //Do analysis and fill aods
1334   //Search for the photon decay in calorimeters
1335   //Read photon list from AOD, produced in class AliAnaPhoton
1336   //Check if 2 photons have the mass of the pi0.
1337
1338   TLorentzVector mom1;
1339   TLorentzVector mom2;
1340   TLorentzVector mom ;
1341
1342   Int_t tag   = 0;
1343   Int_t label = 0;
1344
1345   if(!GetInputAODBranch()){
1346     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1347     abort();
1348   }
1349
1350   //Get shower shape information of clusters
1351   TObjArray *clusters = 0;
1352   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1353   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
1354
1355   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
1356     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1357
1358     //Vertex cut in case of mixed events
1359     Int_t evtIndex1 = 0 ;
1360     if(GetMixedEvent())
1361       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
1362     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
1363     mom1 = *(photon1->Momentum());
1364
1365     //Get original cluster, to recover some information
1366     Int_t iclus = -1;
1367     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1368
1369     if(!cluster1){
1370       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1371       return;
1372     }
1373
1374     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1375     {
1376       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
1377
1378       Int_t evtIndex2 = 0 ;
1379       if(GetMixedEvent())
1380         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1381
1382       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1383         continue ;
1384
1385       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
1386
1387       mom2 = *(photon2->Momentum());
1388
1389       //Get original cluster, to recover some information
1390       Int_t iclus2;
1391       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
1392
1393       if(!cluster2)
1394       {
1395         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1396         return;
1397       }
1398
1399       Float_t e1    = photon1->E();
1400       Float_t e2    = photon2->E();
1401
1402       //Select clusters with good time window difference
1403       Float_t tof1  = cluster1->GetTOF()*1e9;;
1404       Float_t tof2  = cluster2->GetTOF()*1e9;;
1405       Double_t t12diff = tof1-tof2;
1406       fhEPairDiffTime->Fill(e1+e2,    t12diff);
1407       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1408
1409       //Play with the MC stack if available
1410       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1411
1412       // Check the invariant mass for different selection on the local maxima
1413       // Name of AOD method TO BE FIXED
1414       Int_t nMaxima1 = photon1->GetFiducialArea();
1415       Int_t nMaxima2 = photon2->GetFiducialArea();
1416
1417       Double_t mass  = (mom1+mom2).M();
1418       Double_t epair = (mom1+mom2).E();
1419
1420       if(nMaxima1==nMaxima2)
1421       {
1422         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1423         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1424         else                 fhMassPairLocMax[2]->Fill(epair,mass);
1425       }
1426       else if(nMaxima1==1 || nMaxima2==1)
1427       {
1428         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
1429         else                             fhMassPairLocMax[4]->Fill(epair,mass);
1430       }
1431       else
1432         fhMassPairLocMax[5]->Fill(epair,mass);
1433
1434       // combinations with SS axis cut and NLM cut
1435       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1436       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1437       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1438       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1439
1440       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1441       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1442       {
1443         if(GetDebug()>1)
1444           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());
1445
1446         //Fill some histograms about shower shape
1447         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1448         {
1449           FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1450           FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
1451         }
1452
1453         // Tag both photons as decay
1454         photon1->SetTagged(kTRUE);
1455         photon2->SetTagged(kTRUE);
1456
1457         fhPtDecay->Fill(photon1->Pt());
1458         fhEDecay ->Fill(photon1->E() );
1459
1460         fhPtDecay->Fill(photon2->Pt());
1461         fhEDecay ->Fill(photon2->E() );
1462
1463         //Create AOD for analysis
1464         mom = mom1+mom2;
1465
1466         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1467
1468         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1469         pi0.SetDetector(photon1->GetDetector());
1470
1471         // MC
1472         pi0.SetLabel(label);
1473         pi0.SetTag(tag);
1474
1475         //Set the indeces of the original caloclusters
1476         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
1477         //pi0.SetInputFileIndex(input);
1478
1479         AddAODParticle(pi0);
1480
1481       }//pi0
1482
1483     }//2n photon loop
1484
1485   }//1st photon loop
1486
1487   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
1488
1489 }
1490
1491 //__________________________________________________
1492 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1493 {
1494   //Do analysis and fill aods
1495   //Search for the photon decay in calorimeters
1496   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1497   //Check if 2 photons have the mass of the pi0.
1498
1499   TLorentzVector mom1;
1500   TLorentzVector mom2;
1501   TLorentzVector mom ;
1502   Int_t tag   = 0;
1503   Int_t label = 0;
1504   Int_t evtIndex = 0;
1505
1506   // Check calorimeter input
1507   if(!GetInputAODBranch()){
1508     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
1509     abort();
1510   }
1511
1512   // Get the array with conversion photons
1513   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1514   if(!inputAODGammaConv) {
1515
1516     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1517
1518     if(!inputAODGammaConv) {
1519       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1520
1521       return;
1522     }
1523   }
1524
1525   //Get shower shape information of clusters
1526   TObjArray *clusters = 0;
1527   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1528   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
1529
1530   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
1531   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
1532   if(nCTS<=0 || nCalo <=0) {
1533     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1534     return;
1535   }
1536
1537   if(GetDebug() > 1)
1538     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1539
1540   // Do the loop, first calo, second CTS
1541   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1542     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1543     mom1 = *(photon1->Momentum());
1544
1545     //Get original cluster, to recover some information
1546     Int_t iclus = -1;
1547     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1548
1549     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1550       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
1551       if(GetMixedEvent())
1552         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1553       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
1554
1555       mom2 = *(photon2->Momentum());
1556
1557       Double_t mass  = (mom1+mom2).M();
1558       Double_t epair = (mom1+mom2).E();
1559
1560       Int_t nMaxima = photon1->GetFiducialArea();
1561       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1562       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1563       else                fhMassPairLocMax[2]->Fill(epair,mass);
1564
1565       //Play with the MC stack if available
1566       if(IsDataMC())
1567       {
1568         Int_t   label2 = photon2->GetLabel();
1569         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1570
1571         HasPairSameMCMother(photon1, photon2, label, tag) ;
1572       }
1573
1574       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1575       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1576       {
1577         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());
1578
1579         //Fill some histograms about shower shape
1580         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1581         {
1582           FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1583         }
1584
1585         // Tag both photons as decay
1586         photon1->SetTagged(kTRUE);
1587         photon2->SetTagged(kTRUE);
1588
1589         fhPtDecay->Fill(photon1->Pt());
1590         fhEDecay ->Fill(photon1->E() );
1591
1592         //Create AOD for analysis
1593
1594         mom = mom1+mom2;
1595
1596         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1597
1598         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1599         pi0.SetDetector(photon1->GetDetector());
1600
1601         // MC
1602         pi0.SetLabel(label);
1603         pi0.SetTag(tag);
1604
1605         //Set the indeces of the original tracks or caloclusters
1606         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1607         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
1608         //pi0.SetInputFileIndex(input);
1609
1610         AddAODParticle(pi0);
1611
1612       }//pi0
1613     }//2n photon loop
1614
1615   }//1st photon loop
1616
1617   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
1618
1619 }
1620
1621
1622 //_________________________________________________
1623 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
1624 {
1625   //Search for pi0 in fCalorimeter with shower shape analysis
1626
1627   TObjArray * pl        = 0x0;
1628   AliVCaloCells * cells = 0x0;
1629   //Select the Calorimeter of the photon
1630   if      (fCalorimeter == "PHOS" )
1631   {
1632     pl    = GetPHOSClusters();
1633     cells = GetPHOSCells();
1634   }
1635   else if (fCalorimeter == "EMCAL")
1636   {
1637     pl    = GetEMCALClusters();
1638     cells = GetEMCALCells();
1639   }
1640
1641   if(!pl)
1642   {
1643     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1644     return;
1645   }
1646
1647   TLorentzVector mom ;
1648   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
1649   {
1650     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1651
1652     Int_t evtIndex = 0 ;
1653     if (GetMixedEvent())
1654     {
1655       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1656     }
1657
1658     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
1659
1660     //Get Momentum vector,
1661     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1662     {
1663       calo->GetMomentum(mom,GetVertex(evtIndex)) ;
1664     }//Assume that come from vertex in straight line
1665     else
1666     {
1667       Double_t vertex[]={0,0,0};
1668       calo->GetMomentum(mom,vertex) ;
1669     }
1670
1671     //If too small or big pt, skip it
1672     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
1673
1674     //Check acceptance selection
1675     if(IsFiducialCutOn())
1676     {
1677       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
1678       if(! in ) continue ;
1679     }
1680
1681     //Create AOD for analysis
1682     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
1683     aodpi0.SetLabel(calo->GetLabel());
1684
1685     //Set the indeces of the original caloclusters
1686     aodpi0.SetCaloLabel(calo->GetID(),-1);
1687     aodpi0.SetDetector(fCalorimeter);
1688     if(GetDebug() > 1)
1689       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());
1690
1691     //Check Distance to Bad channel, set bit.
1692     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1693     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
1694     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
1695       continue ;
1696
1697     //.......................................
1698     // TOF cut, BE CAREFUL WITH THIS CUT
1699     Double_t tof = calo->GetTOF()*1e9;
1700     if(tof < fTimeCutMin || tof > fTimeCutMax) continue ;
1701
1702
1703     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
1704
1705     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
1706     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
1707     else                         aodpi0.SetDistToBad(0) ;
1708
1709     //Check PID
1710     //PID selection or bit setting
1711     Int_t    nMaxima = 0 ;
1712     Double_t mass    = 0 , angle = 0;
1713     Double_t e1      = 0 , e2    = 0;
1714     //Skip matched clusters with tracks
1715     if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
1716
1717     // Check if cluster is pi0 via cluster splitting
1718     aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleTypeFromClusterSplitting(calo,cells,GetCaloUtils(),
1719                                                                                                  GetVertex(evtIndex),nMaxima,
1720                                                                                                  mass,angle,e1,e2));
1721
1722     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
1723
1724     // If cluster does not pass pid, not pi0, skip it.
1725     // TO DO, add option for Eta ... or conversions
1726     if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
1727
1728     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",
1729                               aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
1730
1731     //Play with the MC stack if available
1732     //Check origin of the candidates
1733     Int_t tag   = 0 ;
1734     if(IsDataMC())
1735     {
1736       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1737           GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1738       {
1739         //aodpi0.SetInputFileIndex(input);
1740         tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodpi0.GetInputFileIndex());
1741         //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
1742         aodpi0.SetTag(tag);
1743         if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
1744       }
1745     }//Work with stack also
1746
1747     //Fill some histograms about shower shape
1748     if(fFillSelectClHisto && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1749     {
1750       Float_t asy =-10;
1751       if(e1+e2 > 0 ) asy = (e1-e2) / (e1+e2);
1752       FillSelectedClusterHistograms(calo, nMaxima, tag, asy);
1753     }
1754
1755     //Add AOD with pi0 object to aod branch
1756     AddAODParticle(aodpi0);
1757
1758   }//loop
1759
1760   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
1761
1762 }
1763 //______________________________________________
1764 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
1765 {
1766   //Do analysis and fill histograms
1767
1768   if(!GetOutputAODBranch())
1769   {
1770     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
1771     abort();
1772   }
1773   //Loop on stored AOD pi0
1774   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1775   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1776
1777   for(Int_t iaod = 0; iaod < naod ; iaod++)
1778   {
1779
1780     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1781     Int_t pdg = pi0->GetIdentifiedParticleType();
1782
1783     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
1784
1785     //Fill pi0 histograms
1786     Float_t ener  = pi0->E();
1787     Float_t pt    = pi0->Pt();
1788     Float_t phi   = pi0->Phi();
1789     if(phi < 0) phi+=TMath::TwoPi();
1790     Float_t eta = pi0->Eta();
1791
1792     fhPt     ->Fill(pt);
1793     fhE      ->Fill(ener);
1794
1795     fhEEta   ->Fill(ener,eta);
1796     fhEPhi   ->Fill(ener,phi);
1797     fhEtaPhi ->Fill(eta,phi);
1798
1799     if(IsDataMC())
1800     {
1801       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1802          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1803         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0))
1804         {
1805           fhPtMC  ->Fill(pt);
1806           fhPhiMC ->Fill(pt,phi);
1807           fhEtaMC ->Fill(pt,eta);
1808         }
1809         else
1810         {
1811           fhPtMCNo  ->Fill(pt);
1812           fhPhiMCNo ->Fill(pt,phi);
1813           fhEtaMCNo ->Fill(pt,eta);
1814         }
1815       }
1816     }//Histograms with MC
1817
1818   }// aod loop
1819
1820 }
1821
1822 //__________________________________________________________________
1823 void AliAnaPi0EbE::Print(const Option_t * opt) const
1824 {
1825   //Print some relevant parameters set for the analysis
1826   if(! opt)
1827     return;
1828
1829   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1830   AliAnaCaloTrackCorrBaseClass::Print("");
1831   printf("Analysis Type = %d \n",  fAnaType) ;
1832   if(fAnaType == kSSCalo){
1833     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1834     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1835     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1836     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1837   }
1838   printf("    \n") ;
1839
1840 }
1841
1842