PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

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