PWGGA/CaloTrackCorrelations/AliAnaPi0EbE.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 // Class for the analysis of high pT pi0 event by event
  18 // Pi0/Eta identified by one of the following:
  19 //  -Invariant mass of 2 cluster in calorimeter
  20 //  -Shower shape analysis in calorimeter
  21 //  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
  22 //
  23 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TList.h>
  29 #include <TClonesArray.h>
  30 #include <TObjString.h>
  31
  32 // --- Analysis system ---
  33 #include "AliAnaPi0EbE.h"
  34 #include "AliCaloTrackReader.h"
  35 #include "AliIsolationCut.h"
  36 #include "AliNeutralMesonSelection.h"
  37 #include "AliCaloPID.h"
  38 #include "AliMCAnalysisUtils.h"
  39 #include "AliStack.h"
  40 #include "AliFiducialCut.h"
  41 #include "TParticle.h"
  42 #include "AliVCluster.h"
  43 #include "AliAODEvent.h"
  44 #include "AliAODMCParticle.h"
  45
  46 ClassImp(AliAnaPi0EbE)
  47
  48 //____________________________
  49 AliAnaPi0EbE::AliAnaPi0EbE() :
  50     AliAnaCaloTrackCorrBaseClass(),fAnaType(kIMCalo),            fCalorimeter(""),
  51     fMinDist(0.),fMinDist2(0.),    fMinDist3(0.),
  52     fFillWeightHistograms(kFALSE), fFillTMHisto(0),              fFillSelectClHisto(0),
  53     fInputAODGammaConvName(""),
  54     // Histograms
  55     fhPt(0),                       fhE(0),
  56     fhEEta(0),                     fhEPhi(0),                    fhEtaPhi(0),
  57     fhPtDecay(0),                  fhEDecay(0),
  58     // Shower shape histos
  59     fhEDispersion(0),              fhELambda0(0),                fhELambda1(0),
  60     fhELambda0NoTRD(0),            fhELambda0FracMaxCellCut(0),
  61     fhEFracMaxCell(0),             fhEFracMaxCellNoTRD(0),
  62     fhENCells(0),                  fhETime(0),                   fhEPairDiffTime(0),
  63     // MC histos
  64     fhPtMCNo(0),                   fhPhiMCNo(0),                 fhEtaMCNo(0),
  65     fhPtMC(0),                     fhPhiMC(0),                   fhEtaMC(0),
  66     fhMassPairMCPi0(0),            fhMassPairMCEta(0),
  67     fhAnglePairMCPi0(0),           fhAnglePairMCEta(0),
  68     // Weight studies
  69     fhECellClusterRatio(0),        fhECellClusterLogRatio(0),
  70     fhEMaxCellClusterRatio(0),     fhEMaxCellClusterLogRatio(0),
  71     fhTrackMatchedDEta(0),         fhTrackMatchedDPhi(0),        fhTrackMatchedDEtaDPhi(0),
  72     fhTrackMatchedMCParticle(0),   fhdEdx(0),
  73     fhEOverP(0),                   fhEOverPNoTRD(0),
  74     // Number of local maxima in cluster
  75     fhNLocMax(0),
  76     fhELambda0LocMax1(0),          fhELambda1LocMax1(0),
  77     fhELambda0LocMax2(0),          fhELambda1LocMax2(0),
  78     fhELambda0LocMaxN(0),          fhELambda1LocMaxN(0)
  79 {
  80   //default ctor
  81
  82   for(Int_t i = 0; i < 6; i++){
  83     fhEMCLambda0[i]     = 0;
  84     fhEMCLambda0NoTRD[i]= 0;
  85     fhEMCLambda0FracMaxCellCut[i]= 0;
  86     fhEMCFracMaxCell[i] = 0;
  87     fhEMCLambda1[i]     = 0;
  88     fhEMCDispersion[i]  = 0;
  89   }
  90
  91   //Weight studies
  92   for(Int_t i =0; i < 14; i++){
  93     fhLambda0ForW0[i] = 0;
  94     //fhLambda1ForW0[i] = 0;
  95     if(i<8)fhMassPairLocMax[i] = 0;
  96   }
  97
  98   //Initialize parameters
  99   InitParameters();
 100
 101 }
 102
 103 //_____________________________________________________________________________________
 104 void AliAnaPi0EbE::FillSelectedClusterHistograms(AliVCluster* cluster,
 105                                                  const Int_t nMaxima,
 106                                                  const Int_t tag)
 107 {
 108   // Fill shower shape, timing and other histograms for selected clusters from decay
 109
 110   Float_t e    = cluster->E();
 111   Float_t disp = cluster->GetDispersion()*cluster->GetDispersion();
 112   Float_t l0   = cluster->GetM02();
 113   Float_t l1   = cluster->GetM20();
 114   Int_t   nSM  = GetModuleNumber(cluster);
 115
 116   AliVCaloCells * cell = 0x0;
 117   if(fCalorimeter == "PHOS")
 118     cell = GetPHOSCells();
 119   else
 120     cell = GetEMCALCells();
 121
 122   Float_t maxCellFraction = 0;
 123   GetCaloUtils()->GetMaxEnergyCell(cell, cluster, maxCellFraction);
 124   fhEFracMaxCell->Fill(e,maxCellFraction);
 125
 126   FillWeightHistograms(cluster);
 127
 128   fhEDispersion->Fill(e, disp);
 129   fhELambda0   ->Fill(e, l0  );
 130   fhELambda1   ->Fill(e, l1  );
 131
 132   fhNLocMax->Fill(e,nMaxima);
 133   if     (nMaxima==1) { fhELambda0LocMax1->Fill(e,l0); fhELambda1LocMax1->Fill(e,l1); }
 134   else if(nMaxima==2) { fhELambda0LocMax2->Fill(e,l0); fhELambda1LocMax2->Fill(e,l1); }
 135   else                { fhELambda0LocMaxN->Fill(e,l0); fhELambda1LocMaxN->Fill(e,l1); }
 136
 137   if(fCalorimeter=="EMCAL" && nSM < 6)
 138   {
 139     fhELambda0NoTRD->Fill(e, l0  );
 140     fhEFracMaxCellNoTRD->Fill(e,maxCellFraction);
 141   }
 142
 143   if(maxCellFraction < 0.5)
 144     fhELambda0FracMaxCellCut->Fill(e, l0  );
 145
 146   fhETime  ->Fill(e, cluster->GetTOF()*1.e9);
 147   fhENCells->Fill(e, cluster->GetNCells());
 148
 149   // Fill Track matching control histograms
 150   if(fFillTMHisto)
 151   {
 152     Float_t dZ  = cluster->GetTrackDz();
 153     Float_t dR  = cluster->GetTrackDx();
 154
 155     if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
 156     {
 157       dR = 2000., dZ = 2000.;
 158       GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
 159     }
 160     //printf("Pi0EbE: dPhi %f, dEta %f\n",dR,dZ);
 161
 162     if(fhTrackMatchedDEta && TMath::Abs(dR) < 999)
 163     {
 164       fhTrackMatchedDEta->Fill(e,dZ);
 165       fhTrackMatchedDPhi->Fill(e,dR);
 166       if(e > 0.5) fhTrackMatchedDEtaDPhi->Fill(dZ,dR);
 167     }
 168
 169     // Check dEdx and E/p of matched clusters
 170
 171     if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
 172     {
 173       AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
 174
 175       if(track) {
 176
 177         Float_t dEdx = track->GetTPCsignal();
 178         fhdEdx->Fill(e, dEdx);
 179
 180         Float_t eOverp = e/track->P();
 181         fhEOverP->Fill(e,  eOverp);
 182
 183         if(fCalorimeter=="EMCAL" && nSM < 6) fhEOverPNoTRD->Fill(e,  eOverp);
 184
 185       }
 186       //else
 187       //  printf("AliAnaPi0EbE::FillSelectedClusterHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
 188
 189
 190
 191       if(IsDataMC())
 192       {
 193         if  ( !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion)  )
 194         {
 195           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 196                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 2.5 );
 197           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 0.5 );
 198           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 1.5 );
 199           else                                                                                 fhTrackMatchedMCParticle->Fill(e, 3.5 );
 200
 201         }
 202         else
 203         {
 204           if       ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)      ||
 205                      GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle->Fill(e, 6.5 );
 206           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle->Fill(e, 4.5 );
 207           else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle->Fill(e, 5.5 );
 208           else                                                                                 fhTrackMatchedMCParticle->Fill(e, 7.5 );
 209         }
 210       }  // MC
 211     }
 212   }// Track matching histograms
 213
 214   if(IsDataMC())
 215   {
 216     //Photon1
 217     if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0)  )
 218     {
 219       fhEMCLambda0[kmcPi0]    ->Fill(e, l0);
 220       fhEMCLambda1[kmcPi0]    ->Fill(e, l1);
 221       fhEMCDispersion[kmcPi0] ->Fill(e, disp);
 222
 223       fhEMCFracMaxCell[kmcPi0]->Fill(e,maxCellFraction);
 224       if(fCalorimeter=="EMCAL" && nSM < 6)
 225         fhEMCLambda0NoTRD[kmcPi0]->Fill(e, l0  );
 226       if(maxCellFraction < 0.5)
 227         fhEMCLambda0FracMaxCellCut[kmcPi0]->Fill(e, l0  );
 228
 229     }//pi0
 230     else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta)  )
 231     {
 232       fhEMCLambda0[kmcEta]    ->Fill(e, l0);
 233       fhEMCLambda1[kmcEta]    ->Fill(e, l1);
 234       fhEMCDispersion[kmcEta] ->Fill(e, disp);
 235       fhEMCFracMaxCell[kmcEta]->Fill(e,maxCellFraction);
 236       if(fCalorimeter=="EMCAL" && nSM < 6)
 237         fhEMCLambda0NoTRD[kmcEta]->Fill(e, l0  );
 238       if(maxCellFraction < 0.5)
 239         fhEMCLambda0FracMaxCellCut[kmcEta]->Fill(e, l0  );
 240     }//eta
 241     else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) &&
 242               GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) )
 243     {
 244       fhEMCLambda0[kmcConversion]    ->Fill(e, l0);
 245       fhEMCLambda1[kmcConversion]    ->Fill(e, l1);
 246       fhEMCDispersion[kmcConversion] ->Fill(e, disp);
 247       fhEMCFracMaxCell[kmcConversion]->Fill(e,maxCellFraction);
 248       if(fCalorimeter=="EMCAL" && nSM < 6)
 249         fhEMCLambda0NoTRD[kmcConversion]->Fill(e, l0  );
 250       if(maxCellFraction < 0.5)
 251         fhEMCLambda0FracMaxCellCut[kmcConversion]->Fill(e, l0  );
 252     }//conversion photon
 253     else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) )
 254     {
 255       fhEMCLambda0[kmcPhoton]    ->Fill(e, l0);
 256       fhEMCLambda1[kmcPhoton]    ->Fill(e, l1);
 257       fhEMCDispersion[kmcPhoton] ->Fill(e, disp);
 258       fhEMCFracMaxCell[kmcPhoton]->Fill(e,maxCellFraction);
 259       if(fCalorimeter=="EMCAL" && nSM < 6)
 260         fhEMCLambda0NoTRD[kmcPhoton]->Fill(e, l0  );
 261       if(maxCellFraction < 0.5)
 262         fhEMCLambda0FracMaxCellCut[kmcPhoton]->Fill(e, l0  );
 263     }//photon   no conversion
 264     else if  ( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
 265     {
 266       fhEMCLambda0[kmcElectron]    ->Fill(e, l0);
 267       fhEMCLambda1[kmcElectron]    ->Fill(e, l1);
 268       fhEMCDispersion[kmcElectron] ->Fill(e, disp);
 269       fhEMCFracMaxCell[kmcElectron]->Fill(e,maxCellFraction);
 270       if(fCalorimeter=="EMCAL" && nSM < 6)
 271         fhEMCLambda0NoTRD[kmcElectron]->Fill(e, l0  );
 272       if(maxCellFraction < 0.5)
 273         fhEMCLambda0FracMaxCellCut[kmcElectron]->Fill(e, l0  );
 274     }//electron
 275     else
 276     {
 277       fhEMCLambda0[kmcHadron]    ->Fill(e, l0);
 278       fhEMCLambda1[kmcHadron]    ->Fill(e, l1);
 279       fhEMCDispersion[kmcHadron] ->Fill(e, disp);
 280       fhEMCFracMaxCell[kmcHadron]->Fill(e,maxCellFraction);
 281       if(fCalorimeter=="EMCAL" && nSM < 6)
 282         fhEMCLambda0NoTRD[kmcHadron]->Fill(e, l0  );
 283       if(maxCellFraction < 0.5)
 284         fhEMCLambda0FracMaxCellCut[kmcHadron]->Fill(e, l0  );
 285     }//other particles
 286   }//MC
 287 }
 288
 289 //________________________________________________________
 290 void AliAnaPi0EbE::FillWeightHistograms(AliVCluster *clus)
 291 {
 292   // Calculate weights and fill histograms
 293
 294   if(!fFillWeightHistograms || GetMixedEvent()) return;
 295
 296   AliVCaloCells* cells = 0;
 297   if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
 298   else                        cells = GetPHOSCells();
 299
 300   // First recalculate energy in case non linearity was applied
 301   Float_t  energy = 0;
 302   Float_t  ampMax = 0;
 303   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 304   {
 305
 306     Int_t id       = clus->GetCellsAbsId()[ipos];
 307
 308     //Recalibrate cell energy if needed
 309     Float_t amp = cells->GetCellAmplitude(id);
 310     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 311
 312     energy    += amp;
 313
 314     if(amp> ampMax)
 315       ampMax = amp;
 316
 317   } // energy loop
 318
 319   if(energy <=0 )
 320   {
 321     printf("AliAnaPi0EbE::WeightHistograms()- Wrong calculated energy %f\n",energy);
 322     return;
 323   }
 324
 325   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
 326   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
 327
 328   //Get the ratio and log ratio to all cells in cluster
 329   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++)
 330   {
 331     Int_t id       = clus->GetCellsAbsId()[ipos];
 332
 333     //Recalibrate cell energy if needed
 334     Float_t amp = cells->GetCellAmplitude(id);
 335     GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
 336
 337     fhECellClusterRatio   ->Fill(energy,amp/energy);
 338     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
 339   }
 340
 341   //Recalculate shower shape for different W0
 342   if(fCalorimeter=="EMCAL"){
 343
 344     Float_t l0org = clus->GetM02();
 345     Float_t l1org = clus->GetM20();
 346     Float_t dorg  = clus->GetDispersion();
 347
 348     for(Int_t iw = 0; iw < 14; iw++)
 349     {
 350       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
 351       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
 352
 353       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
 354       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
 355
 356     } // w0 loop
 357
 358     // Set the original values back
 359     clus->SetM02(l0org);
 360     clus->SetM20(l1org);
 361     clus->SetDispersion(dorg);
 362
 363   }// EMCAL
 364 }
 365
 366 //__________________________________________
 367 TObjString * AliAnaPi0EbE::GetAnalysisCuts()
 368 {
 369         //Save parameters used for analysis
 370   TString parList ; //this will be list of parameters used for this analysis.
 371   const Int_t buffersize = 255;
 372   char onePar[buffersize] ;
 373
 374   snprintf(onePar,buffersize,"--- AliAnaPi0EbE ---\n") ;
 375   parList+=onePar ;
 376   snprintf(onePar,buffersize,"fAnaType=%d (Pi0 selection type) \n",fAnaType) ;
 377   parList+=onePar ;
 378
 379   if(fAnaType == kSSCalo)
 380   {
 381     snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 382     parList+=onePar ;
 383     snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 384     parList+=onePar ;
 385     snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 386     parList+=onePar ;
 387     snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 388     parList+=onePar ;
 389   }
 390
 391   //Get parameters set in base class.
 392   parList += GetBaseParametersList() ;
 393
 394   //Get parameters set in PID class.
 395   if(fAnaType == kSSCalo) parList += GetCaloPID()->GetPIDParametersList() ;
 396
 397   return new TObjString(parList) ;
 398 }
 399
 400 //__________________________________________________________________
 401 void AliAnaPi0EbE::HasPairSameMCMother(AliAODPWG4Particle * photon1,
 402                                        AliAODPWG4Particle * photon2,
 403                                        Int_t & label, Int_t & tag)
 404 {
 405   // Check the labels of pare in case mother was same pi0 or eta
 406   // Set the new AOD accordingly
 407
 408   Int_t  label1 = photon1->GetLabel();
 409   Int_t  label2 = photon2->GetLabel();
 410
 411   if(label1 < 0 || label2 < 0 ) return ;
 412
 413   //Int_t tag1 = GetMCAnalysisUtils()->CheckOrigin(label1, GetReader(), photon1->GetInputFileIndex());
 414   //Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex());
 415   Int_t tag1 = photon1->GetTag();
 416   Int_t tag2 = photon2->GetTag();
 417
 418   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Origin of: photon1 %d; photon2 %d \n",tag1, tag2);
 419   if( (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay) &&
 420        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCPi0Decay)    ) ||
 421       (GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCEtaDecay) &&
 422        GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCEtaDecay)    )
 423      )
 424   {
 425
 426     //Check if pi0/eta mother is the same
 427     if(GetReader()->ReadStack())
 428     {
 429       if(label1>=0)
 430       {
 431         TParticle * mother1 = GetMCStack()->Particle(label1);//photon in kine tree
 432         label1 = mother1->GetFirstMother();
 433         //mother1 = GetMCStack()->Particle(label1);//pi0
 434       }
 435       if(label2>=0)
 436       {
 437         TParticle * mother2 = GetMCStack()->Particle(label2);//photon in kine tree
 438         label2 = mother2->GetFirstMother();
 439         //mother2 = GetMCStack()->Particle(label2);//pi0
 440       }
 441     } // STACK
 442     else if(GetReader()->ReadAODMCParticles())
 443     {//&& (input > -1)){
 444       if(label1>=0)
 445       {
 446         AliAODMCParticle * mother1 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon1->GetInputFileIndex()))->At(label1);//photon in kine tree
 447         label1 = mother1->GetMother();
 448         //mother1 = GetMCStack()->Particle(label1);//pi0
 449       }
 450       if(label2>=0)
 451       {
 452         AliAODMCParticle * mother2 = (AliAODMCParticle *) (GetReader()->GetAODMCParticles(photon2->GetInputFileIndex()))->At(label2);//photon in kine tree
 453         label2 = mother2->GetMother();
 454         //mother2 = GetMCStack()->Particle(label2);//pi0
 455       }
 456     }// AOD
 457
 458     //printf("mother1 %d, mother2 %d\n",label1,label2);
 459     if( label1 == label2 && label1>=0 )
 460     {
 461
 462       label = label1;
 463
 464       TLorentzVector mom1 = *(photon1->Momentum());
 465       TLorentzVector mom2 = *(photon2->Momentum());
 466
 467       Double_t angle = mom2.Angle(mom1.Vect());
 468       Double_t mass  = (mom1+mom2).M();
 469       Double_t epair = (mom1+mom2).E();
 470
 471       if(GetMCAnalysisUtils()->CheckTagBit(tag1,AliMCAnalysisUtils::kMCPi0Decay))
 472       {
 473         fhMassPairMCPi0 ->Fill(epair,mass);
 474         fhAnglePairMCPi0->Fill(epair,angle);
 475         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCPi0);
 476       }
 477       else
 478       {
 479         fhMassPairMCEta ->Fill(epair,mass);
 480         fhAnglePairMCEta->Fill(epair,angle);
 481         GetMCAnalysisUtils()->SetTagBit(tag,AliMCAnalysisUtils::kMCEta);
 482       }
 483
 484     } // same label
 485   } // both from eta or pi0 decay
 486
 487 }
 488
 489 //_____________________________________________
 490 TList *  AliAnaPi0EbE::GetCreateOutputObjects()
 491 {
 492   // Create histograms to be saved in output file and
 493   // store them in outputContainer
 494   TList * outputContainer = new TList() ;
 495   outputContainer->SetName("Pi0EbEHistos") ;
 496
 497   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();
 498   Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();
 499   Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 500   Int_t ssbins   = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax  = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin  = GetHistogramRanges()->GetHistoShowerShapeMin();
 501   Int_t tdbins   = GetHistogramRanges()->GetHistoDiffTimeBins() ;    Float_t tdmax  = GetHistogramRanges()->GetHistoDiffTimeMax();     Float_t tdmin  = GetHistogramRanges()->GetHistoDiffTimeMin();
 502   Int_t tbins    = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax   = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin   = GetHistogramRanges()->GetHistoTimeMin();
 503   Int_t nbins    = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax   = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin   = GetHistogramRanges()->GetHistoNClusterCellMin();
 504
 505   Int_t   nmassbins   = GetHistogramRanges()->GetHistoMassBins();
 506   Float_t massmin     = GetHistogramRanges()->GetHistoMassMin();
 507   Float_t massmax     = GetHistogramRanges()->GetHistoMassMax();
 508
 509   Int_t   nresetabins = GetHistogramRanges()->GetHistoTrackResidualEtaBins();
 510   Float_t resetamax   = GetHistogramRanges()->GetHistoTrackResidualEtaMax();
 511   Float_t resetamin   = GetHistogramRanges()->GetHistoTrackResidualEtaMin();
 512   Int_t   nresphibins = GetHistogramRanges()->GetHistoTrackResidualPhiBins();
 513   Float_t resphimax   = GetHistogramRanges()->GetHistoTrackResidualPhiMax();
 514   Float_t resphimin   = GetHistogramRanges()->GetHistoTrackResidualPhiMin();
 515
 516   Int_t   ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();
 517   Float_t dedxmax     = GetHistogramRanges()->GetHistodEdxMax();
 518   Float_t dedxmin     = GetHistogramRanges()->GetHistodEdxMin();
 519   Int_t   nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();
 520   Float_t pOverEmax   = GetHistogramRanges()->GetHistoPOverEMax();
 521   Float_t pOverEmin   = GetHistogramRanges()->GetHistoPOverEMin();
 522
 523
 524   fhPt  = new TH1F("hPt","Number of identified  #pi^{0} (#eta) decay",nptbins,ptmin,ptmax);
 525   fhPt->SetYTitle("N");
 526   fhPt->SetXTitle("p_{T} (GeV/c)");
 527   outputContainer->Add(fhPt) ;
 528
 529   fhE  = new TH1F("hE","Number of identified  #pi^{0} (#eta) decay pairs",nptbins,ptmin,ptmax);
 530   fhE->SetYTitle("N");
 531   fhE->SetXTitle("E (GeV)");
 532   outputContainer->Add(fhE) ;
 533
 534   fhEPhi  = new TH2F
 535   ("hEPhi","Selected #pi^{0} (#eta) pairs: E vs #phi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 536   fhEPhi->SetYTitle("#phi (rad)");
 537   fhEPhi->SetXTitle("E (GeV)");
 538   outputContainer->Add(fhEPhi) ;
 539
 540   fhEEta  = new TH2F
 541   ("hEEta","Selected #pi^{0} (#eta) pairs: E vs #eta",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 542   fhEEta->SetYTitle("#eta");
 543   fhEEta->SetXTitle("E (GeV)");
 544   outputContainer->Add(fhEEta) ;
 545
 546   fhEtaPhi  = new TH2F
 547   ("hEtaPhi","Selected #pi^{0} (#eta) pairs: #eta vs #phi",netabins,etamin,etamax, nphibins,phimin,phimax);
 548   fhEtaPhi->SetYTitle("#phi (rad)");
 549   fhEtaPhi->SetXTitle("#eta");
 550   outputContainer->Add(fhEtaPhi) ;
 551
 552   fhPtDecay  = new TH1F("hPtDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 553   fhPtDecay->SetYTitle("N");
 554   fhPtDecay->SetXTitle("p_{T} (GeV/c)");
 555   outputContainer->Add(fhPtDecay) ;
 556
 557   fhEDecay  = new TH1F("hEDecay","Number of identified  #pi^{0} (#eta) decay photons",nptbins,ptmin,ptmax);
 558   fhEDecay->SetYTitle("N");
 559   fhEDecay->SetXTitle("E (GeV)");
 560   outputContainer->Add(fhEDecay) ;
 561
 562   ////////
 563
 564   if( fFillSelectClHisto && (fAnaType == kIMCalo || fAnaType == kIMCaloTracks) )
 565   {
 566
 567     fhEDispersion  = new TH2F
 568     ("hEDispersion","Selected #pi^{0} (#eta) pairs: E vs dispersion",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 569     fhEDispersion->SetYTitle("D^{2}");
 570     fhEDispersion->SetXTitle("E (GeV)");
 571     outputContainer->Add(fhEDispersion) ;
 572
 573     fhELambda0  = new TH2F
 574     ("hELambda0","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 575     fhELambda0->SetYTitle("#lambda_{0}^{2}");
 576     fhELambda0->SetXTitle("E (GeV)");
 577     outputContainer->Add(fhELambda0) ;
 578
 579     fhELambda1  = new TH2F
 580     ("hELambda1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 581     fhELambda1->SetYTitle("#lambda_{1}^{2}");
 582     fhELambda1->SetXTitle("E (GeV)");
 583     outputContainer->Add(fhELambda1) ;
 584
 585     fhELambda0FracMaxCellCut  = new TH2F
 586     ("hELambda0FracMaxCellCut","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy < 0.5",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 587     fhELambda0FracMaxCellCut->SetYTitle("#lambda_{0}^{2}");
 588     fhELambda0FracMaxCellCut->SetXTitle("E (GeV)");
 589     outputContainer->Add(fhELambda0FracMaxCellCut) ;
 590
 591     fhEFracMaxCell  = new TH2F
 592     ("hEFracMaxCell","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy",nptbins,ptmin,ptmax,100,0,1);
 593     fhEFracMaxCell->SetYTitle("Fraction");
 594     fhEFracMaxCell->SetXTitle("E (GeV)");
 595     outputContainer->Add(fhEFracMaxCell) ;
 596
 597     fhNLocMax = new TH2F("hNLocMax","Number of local maxima in cluster",
 598                          nptbins,ptmin,ptmax,10,0,10);
 599     fhNLocMax ->SetYTitle("N maxima");
 600     fhNLocMax ->SetXTitle("E (GeV)");
 601     outputContainer->Add(fhNLocMax) ;
 602
 603     fhELambda0LocMax1  = new TH2F
 604     ("hELambda0LocMax1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, 1 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 605     fhELambda0LocMax1->SetYTitle("#lambda_{0}^{2}");
 606     fhELambda0LocMax1->SetXTitle("E (GeV)");
 607     outputContainer->Add(fhELambda0LocMax1) ;
 608
 609     fhELambda1LocMax1  = new TH2F
 610     ("hELambda1LocMax1","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, 1 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 611     fhELambda1LocMax1->SetYTitle("#lambda_{1}^{2}");
 612     fhELambda1LocMax1->SetXTitle("E (GeV)");
 613     outputContainer->Add(fhELambda1LocMax1) ;
 614
 615     fhELambda0LocMax2  = new TH2F
 616     ("hELambda0LocMax2","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, 2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 617     fhELambda0LocMax2->SetYTitle("#lambda_{0}^{2}");
 618     fhELambda0LocMax2->SetXTitle("E (GeV)");
 619     outputContainer->Add(fhELambda0LocMax2) ;
 620
 621     fhELambda1LocMax2  = new TH2F
 622     ("hELambda1LocMax2","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, 2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 623     fhELambda1LocMax2->SetYTitle("#lambda_{1}^{2}");
 624     fhELambda1LocMax2->SetXTitle("E (GeV)");
 625     outputContainer->Add(fhELambda1LocMax2) ;
 626
 627     fhELambda0LocMaxN  = new TH2F
 628     ("hELambda0LocMaxN","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, N>2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 629     fhELambda0LocMaxN->SetYTitle("#lambda_{0}^{2}");
 630     fhELambda0LocMaxN->SetXTitle("E (GeV)");
 631     outputContainer->Add(fhELambda0LocMaxN) ;
 632
 633     fhELambda1LocMaxN  = new TH2F
 634     ("hELambda1LocMaxN","Selected #pi^{0} (#eta) pairs: E vs #lambda_{1}, N>2 Local maxima",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 635     fhELambda1LocMaxN->SetYTitle("#lambda_{1}^{2}");
 636     fhELambda1LocMaxN->SetXTitle("E (GeV)");
 637     outputContainer->Add(fhELambda1LocMaxN) ;
 638
 639
 640     if(fCalorimeter=="EMCAL")
 641     {
 642       fhELambda0NoTRD  = new TH2F
 643       ("hELambda0NoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, not behind TRD",nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 644       fhELambda0NoTRD->SetYTitle("#lambda_{0}^{2}");
 645       fhELambda0NoTRD->SetXTitle("E (GeV)");
 646       outputContainer->Add(fhELambda0NoTRD) ;
 647
 648       fhEFracMaxCellNoTRD  = new TH2F
 649       ("hEFracMaxCellNoTRD","Selected #pi^{0} (#eta) pairs: E vs #lambda_{0}, Max cell fraction of energy, not behind TRD",nptbins,ptmin,ptmax,100,0,1);
 650       fhEFracMaxCellNoTRD->SetYTitle("Fraction");
 651       fhEFracMaxCellNoTRD->SetXTitle("E (GeV)");
 652       outputContainer->Add(fhEFracMaxCellNoTRD) ;
 653     }
 654
 655     fhENCells  = new TH2F ("hENCells","N cells in cluster vs E ", nptbins,ptmin,ptmax, nbins,nmin,nmax);
 656     fhENCells->SetXTitle("E (GeV)");
 657     fhENCells->SetYTitle("# of cells in cluster");
 658     outputContainer->Add(fhENCells);
 659
 660     fhETime = new TH2F("hETime","cluster time vs pair E",nptbins,ptmin,ptmax, tbins,tmin,tmax);
 661     fhETime->SetXTitle("E (GeV)");
 662     fhETime->SetYTitle("t (ns)");
 663     outputContainer->Add(fhETime);
 664
 665   }// Invariant mass analysis in calorimeters and calorimeter + conversion photons
 666
 667   if(fAnaType == kIMCalo)
 668   {
 669     fhEPairDiffTime = new TH2F("hEPairDiffTime","cluster pair time difference vs E",nptbins,ptmin,ptmax, tdbins,tdmin,tdmax);
 670     fhEPairDiffTime->SetXTitle("E_{pair} (GeV)");
 671     fhEPairDiffTime->SetYTitle("#Delta t (ns)");
 672     outputContainer->Add(fhEPairDiffTime);
 673
 674     TString combiName [] = {"1LocMax","2LocMax","NLocMax","1LocMax2LocMax","1LocMaxNLocMax","2LocMaxNLocMax","1LocMaxSSBad","NLocMaxSSGood"};
 675     TString combiTitle[] = {"1 Local Maxima in both clusters","2 Local Maxima in both clusters","more than 2 Local Maxima in both clusters",
 676       "1 Local Maxima paired with 2 Local Maxima","1 Local Maxima paired with more than 2 Local Maxima",
 677       "2 Local Maxima paired with more than 2 Local Maxima",
 678       "1 Local Maxima paired with #lambda_{0}^{2}>0.3","N Local Maxima paired with 0.1<#lambda_{0}^{2}<0.3"};
 679
 680     for (Int_t i = 0; i < 8 ; i++)
 681     {
 682
 683       if (fAnaType == kIMCaloTracks && i > 2 ) continue ;
 684
 685       fhMassPairLocMax[i]  = new TH2F
 686       (Form("MassPairLocMax%s",combiName[i].Data()),
 687        Form("Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}, %s", combiTitle[i].Data()),
 688        nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 689       fhMassPairLocMax[i]->SetYTitle("Mass (MeV/c^{2})");
 690       fhMassPairLocMax[i]->SetXTitle("E_{pair} (GeV)");
 691       outputContainer->Add(fhMassPairLocMax[i]) ;
 692     }
 693   }
 694
 695   if(fFillTMHisto)
 696   {
 697     fhTrackMatchedDEta  = new TH2F
 698     ("hTrackMatchedDEta",
 699      "d#eta of cluster-track vs cluster energy",
 700      nptbins,ptmin,ptmax,nresetabins,resetamin,resetamax);
 701     fhTrackMatchedDEta->SetYTitle("d#eta");
 702     fhTrackMatchedDEta->SetXTitle("E_{cluster} (GeV)");
 703
 704     fhTrackMatchedDPhi  = new TH2F
 705     ("hTrackMatchedDPhi",
 706      "d#phi of cluster-track vs cluster energy",
 707      nptbins,ptmin,ptmax,nresphibins,resphimin,resphimax);
 708     fhTrackMatchedDPhi->SetYTitle("d#phi (rad)");
 709     fhTrackMatchedDPhi->SetXTitle("E_{cluster} (GeV)");
 710
 711     fhTrackMatchedDEtaDPhi  = new TH2F
 712     ("hTrackMatchedDEtaDPhi",
 713      "d#eta vs d#phi of cluster-track vs cluster energy",
 714      nresetabins,resetamin,resetamax,nresphibins,resphimin,resphimax);
 715     fhTrackMatchedDEtaDPhi->SetYTitle("d#phi (rad)");
 716     fhTrackMatchedDEtaDPhi->SetXTitle("d#eta");
 717
 718     outputContainer->Add(fhTrackMatchedDEta) ;
 719     outputContainer->Add(fhTrackMatchedDPhi) ;
 720     outputContainer->Add(fhTrackMatchedDEtaDPhi) ;
 721
 722     fhdEdx  = new TH2F ("hdEdx","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 723     fhdEdx->SetXTitle("E (GeV)");
 724     fhdEdx->SetYTitle("<dE/dx>");
 725     outputContainer->Add(fhdEdx);
 726
 727     fhEOverP  = new TH2F ("hEOverP","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 728     fhEOverP->SetXTitle("E (GeV)");
 729     fhEOverP->SetYTitle("E/p");
 730     outputContainer->Add(fhEOverP);
 731
 732     if(fCalorimeter=="EMCAL")
 733     {
 734       fhEOverPNoTRD  = new TH2F ("hEOverPNoTRD","matched track E/p vs cluster E, SM not behind TRD ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 735       fhEOverPNoTRD->SetXTitle("E (GeV)");
 736       fhEOverPNoTRD->SetYTitle("E/p");
 737       outputContainer->Add(fhEOverPNoTRD);
 738     }
 739
 740     if(IsDataMC())
 741     {
 742       fhTrackMatchedMCParticle  = new TH2F
 743       ("hTrackMatchedMCParticle",
 744        "Origin of particle vs energy",
 745        nptbins,ptmin,ptmax,8,0,8);
 746       fhTrackMatchedMCParticle->SetXTitle("E (GeV)");
 747       //fhTrackMatchedMCParticle->SetYTitle("Particle type");
 748
 749       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(1 ,"Photon");
 750       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(2 ,"Electron");
 751       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(3 ,"Meson Merged");
 752       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(4 ,"Rest");
 753       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(5 ,"Conv. Photon");
 754       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(6 ,"Conv. Electron");
 755       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(7 ,"Conv. Merged");
 756       fhTrackMatchedMCParticle->GetYaxis()->SetBinLabel(8 ,"Conv. Rest");
 757
 758       outputContainer->Add(fhTrackMatchedMCParticle);
 759     }
 760   }
 761
 762   if(fFillWeightHistograms)
 763   {
 764
 765     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
 766                                      nptbins,ptmin,ptmax, 100,0,1.);
 767     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
 768     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
 769     outputContainer->Add(fhECellClusterRatio);
 770
 771     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
 772                                         nptbins,ptmin,ptmax, 100,-10,0);
 773     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
 774     fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
 775     outputContainer->Add(fhECellClusterLogRatio);
 776
 777     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected decay photons from neutral meson",
 778                                         nptbins,ptmin,ptmax, 100,0,1.);
 779     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
 780     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
 781     outputContainer->Add(fhEMaxCellClusterRatio);
 782
 783     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected decay photons from neutral meson",
 784                                            nptbins,ptmin,ptmax, 100,-10,0);
 785     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
 786     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
 787     outputContainer->Add(fhEMaxCellClusterLogRatio);
 788
 789     for(Int_t iw = 0; iw < 14; iw++)
 790     {
 791       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),
 792                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 793       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
 794       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
 795       outputContainer->Add(fhLambda0ForW0[iw]);
 796
 797 //      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),
 798 //                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 799 //      fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
 800 //      fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
 801 //      outputContainer->Add(fhLambda1ForW0[iw]);
 802
 803     }
 804   }
 805
 806   if(IsDataMC())
 807   {
 808     if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
 809        GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 810     {
 811
 812       fhPtMC  = new TH1F("hPtMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax);
 813       fhPtMC->SetYTitle("N");
 814       fhPtMC->SetXTitle("p_{T} (GeV/c)");
 815       outputContainer->Add(fhPtMC) ;
 816
 817       fhPhiMC  = new TH2F
 818       ("hPhiMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 819       fhPhiMC->SetYTitle("#phi");
 820       fhPhiMC->SetXTitle("p_{T} (GeV/c)");
 821       outputContainer->Add(fhPhiMC) ;
 822
 823       fhEtaMC  = new TH2F
 824       ("hEtaMC","Identified #pi^{0} (#eta) from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 825       fhEtaMC->SetYTitle("#eta");
 826       fhEtaMC->SetXTitle("p_{T} (GeV/c)");
 827       outputContainer->Add(fhEtaMC) ;
 828
 829       fhPtMCNo  = new TH1F("hPtMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax);
 830       fhPtMCNo->SetYTitle("N");
 831       fhPtMCNo->SetXTitle("p_{T} (GeV/c)");
 832       outputContainer->Add(fhPtMCNo) ;
 833
 834       fhPhiMCNo  = new TH2F
 835       ("hPhiMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 836       fhPhiMCNo->SetYTitle("#phi");
 837       fhPhiMCNo->SetXTitle("p_{T} (GeV/c)");
 838       outputContainer->Add(fhPhiMCNo) ;
 839
 840       fhEtaMCNo  = new TH2F
 841       ("hEtaMCNo","Identified #pi^{0} (#eta) not from #pi^{0} (#eta)",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 842       fhEtaMCNo->SetYTitle("#eta");
 843       fhEtaMCNo->SetXTitle("p_{T} (GeV/c)");
 844       outputContainer->Add(fhEtaMCNo) ;
 845
 846       fhAnglePairMCPi0  = new TH2F
 847       ("AnglePairMCPi0",
 848        "Angle between decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,250,0,0.5);
 849       fhAnglePairMCPi0->SetYTitle("#alpha (rad)");
 850       fhAnglePairMCPi0->SetXTitle("E_{pair} (GeV)");
 851       outputContainer->Add(fhAnglePairMCPi0) ;
 852
 853       fhAnglePairMCEta  = new TH2F
 854       ("AnglePairMCEta",
 855        "Angle between decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,250,0,0.5);
 856       fhAnglePairMCEta->SetYTitle("#alpha (rad)");
 857       fhAnglePairMCEta->SetXTitle("E_{pair} (GeV)");
 858       outputContainer->Add(fhAnglePairMCEta) ;
 859
 860       fhMassPairMCPi0  = new TH2F
 861       ("MassPairMCPi0",
 862        "Mass for decay #gamma pair vs E_{pair}, origin #pi^{0}",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 863       fhMassPairMCPi0->SetYTitle("Mass (MeV/c^{2})");
 864       fhMassPairMCPi0->SetXTitle("E_{pair} (GeV)");
 865       outputContainer->Add(fhMassPairMCPi0) ;
 866
 867       fhMassPairMCEta  = new TH2F
 868       ("MassPairMCEta",
 869        "Mass for decay #gamma pair vs E_{pair}, origin #eta",nptbins,ptmin,ptmax,nmassbins,massmin,massmax);
 870       fhMassPairMCEta->SetYTitle("Mass (MeV/c^{2})");
 871       fhMassPairMCEta->SetXTitle("E_{pair} (GeV)");
 872       outputContainer->Add(fhMassPairMCEta) ;
 873
 874       if(fAnaType == kIMCalo){
 875         TString ptype[] ={"#gamma","#gamma->e^{#pm}","#pi^{0}","#eta","e^{#pm}", "hadron"};
 876         TString pname[] ={"Photon","Conversion",     "Pi0",    "Eta", "Electron","Hadron"};
 877         for(Int_t i = 0; i < 6; i++){
 878
 879           fhEMCLambda0[i]  = new TH2F(Form("hELambda0_MC%s",pname[i].Data()),
 880                                       Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}",ptype[i].Data()),
 881                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 882           fhEMCLambda0[i]->SetYTitle("#lambda_{0}^{2}");
 883           fhEMCLambda0[i]->SetXTitle("E (GeV)");
 884           outputContainer->Add(fhEMCLambda0[i]) ;
 885
 886           if(fCalorimeter=="EMCAL"){
 887             fhEMCLambda0NoTRD[i]  = new TH2F(Form("hELambda0NoTRD_MC%s",pname[i].Data()),
 888                                              Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, NoTRD",ptype[i].Data()),
 889                                              nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 890             fhEMCLambda0NoTRD[i]->SetYTitle("#lambda_{0}^{2}");
 891             fhEMCLambda0NoTRD[i]->SetXTitle("E (GeV)");
 892             outputContainer->Add(fhEMCLambda0NoTRD[i]) ;
 893           }
 894
 895           fhEMCLambda0FracMaxCellCut[i]  = new TH2F(Form("hELambda0FracMaxCellCut_MC%s",pname[i].Data()),
 896                                                     Form("Selected pair, cluster from %s : E vs #lambda_{0}^{2}, Max cell fraction of energy < 0.5 ",ptype[i].Data()),
 897                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 898           fhEMCLambda0FracMaxCellCut[i]->SetYTitle("#lambda_{0}^{2}");
 899           fhEMCLambda0FracMaxCellCut[i]->SetXTitle("E (GeV)");
 900           outputContainer->Add(fhEMCLambda0FracMaxCellCut[i]) ;
 901
 902           fhEMCFracMaxCell[i]  = new TH2F(Form("hEFracMaxCell_MC%s",pname[i].Data()),
 903                                           Form("Selected pair, cluster from %s : E vs Max cell fraction of energy",ptype[i].Data()),
 904                                           nptbins,ptmin,ptmax,100,0,1);
 905           fhEMCFracMaxCell[i]->SetYTitle("Fraction");
 906           fhEMCFracMaxCell[i]->SetXTitle("E (GeV)");
 907           outputContainer->Add(fhEMCFracMaxCell[i]) ;
 908
 909           fhEMCLambda1[i]  = new TH2F(Form("hELambda1_MC%s",pname[i].Data()),
 910                                       Form("Selected pair, cluster from %s : E vs #lambda_{1}^{2}",ptype[i].Data()),
 911                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 912           fhEMCLambda1[i]->SetYTitle("#lambda_{1}^{2}");
 913           fhEMCLambda1[i]->SetXTitle("E (GeV)");
 914           outputContainer->Add(fhEMCLambda1[i]) ;
 915
 916           fhEMCDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pname[i].Data()),
 917                                          Form("Selected pair, cluster from %s : E vs dispersion^{2}",ptype[i].Data()),
 918                                          nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 919           fhEMCDispersion[i]->SetYTitle("D^{2}");
 920           fhEMCDispersion[i]->SetXTitle("E (GeV)");
 921           outputContainer->Add(fhEMCDispersion[i]) ;
 922
 923         }//
 924
 925       }//kIMCalo
 926     } //Not MC reader
 927   }//Histos with MC
 928
 929
 930   //Keep neutral meson selection histograms if requiered
 931   //Setting done in AliNeutralMesonSelection
 932
 933   if(fAnaType!=kSSCalo && GetNeutralMesonSelection()){
 934
 935     TList * nmsHistos = GetNeutralMesonSelection()->GetCreateOutputObjects() ;
 936     if(GetNeutralMesonSelection()->AreNeutralMesonSelectionHistosKept())
 937       for(Int_t i = 0; i < nmsHistos->GetEntries(); i++) outputContainer->Add(nmsHistos->At(i)) ;
 938     delete nmsHistos;
 939
 940   }
 941
 942   return outputContainer ;
 943
 944 }
 945
 946 //____________________________________________________________________________
 947 void AliAnaPi0EbE::Init()
 948 {
 949   //Init
 950   //Do some checks
 951   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
 952     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
 953     abort();
 954   }
 955   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
 956     printf("AliAnaPi0EbE::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
 957     abort();
 958   }
 959
 960 }
 961
 962 //____________________________________________________________________________
 963 void AliAnaPi0EbE::InitParameters()
 964 {
 965   //Initialize the parameters of the analysis.
 966   AddToHistogramsName("AnaPi0EbE_");
 967
 968   fInputAODGammaConvName = "PhotonsCTS" ;
 969   fAnaType = kIMCalo ;
 970   fCalorimeter = "EMCAL" ;
 971   fMinDist  = 2.;
 972   fMinDist2 = 4.;
 973   fMinDist3 = 5.;
 974
 975 }
 976
 977 //__________________________________________________________________
 978 void  AliAnaPi0EbE::MakeAnalysisFillAOD()
 979 {
 980   //Do analysis and fill aods
 981
 982   switch(fAnaType)
 983   {
 984     case kIMCalo:
 985       MakeInvMassInCalorimeter();
 986       break;
 987
 988     case kSSCalo:
 989       MakeShowerShapeIdentification();
 990       break;
 991
 992     case kIMCaloTracks:
 993       MakeInvMassInCalorimeterAndCTS();
 994       break;
 995
 996   }
 997 }
 998
 999 //____________________________________________
1000 void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
1001 {
1002   //Do analysis and fill aods
1003   //Search for the photon decay in calorimeters
1004   //Read photon list from AOD, produced in class AliAnaPhoton
1005   //Check if 2 photons have the mass of the pi0.
1006
1007   TLorentzVector mom1;
1008   TLorentzVector mom2;
1009   TLorentzVector mom ;
1010
1011   Int_t tag   = 0;
1012   Int_t label = 0;
1013
1014   if(!GetInputAODBranch()){
1015     printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - No input calo photons in AOD with name branch < %s >, STOP \n",GetInputAODName().Data());
1016     abort();
1017   }
1018
1019   //Get shower shape information of clusters
1020   TObjArray *clusters = 0;
1021   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1022   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
1023
1024   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast()-1; iphoton++){
1025     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1026
1027     //Vertex cut in case of mixed events
1028     Int_t evtIndex1 = 0 ;
1029     if(GetMixedEvent())
1030       evtIndex1 = GetMixedEvent()->EventIndexForCaloCluster(photon1->GetCaloLabel(0)) ;
1031     if(TMath::Abs(GetVertex(evtIndex1)[2]) > GetZvertexCut()) continue ;  //vertex cut
1032     mom1 = *(photon1->Momentum());
1033
1034     //Get original cluster, to recover some information
1035     Int_t iclus = -1;
1036     AliVCluster *cluster1 = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1037
1038     if(!cluster1){
1039       printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - First cluster not found\n");
1040       return;
1041     }
1042
1043     for(Int_t jphoton = iphoton+1; jphoton < GetInputAODBranch()->GetEntriesFast(); jphoton++)
1044     {
1045       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(jphoton));
1046
1047       Int_t evtIndex2 = 0 ;
1048       if(GetMixedEvent())
1049         evtIndex2 = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1050
1051       if(GetMixedEvent() && (evtIndex1 == evtIndex2))
1052         continue ;
1053
1054       if(TMath::Abs(GetVertex(evtIndex2)[2]) > GetZvertexCut()) continue ;  //vertex cut
1055
1056       mom2 = *(photon2->Momentum());
1057
1058       //Get original cluster, to recover some information
1059       Int_t iclus2;
1060       AliVCluster *cluster2 = FindCluster(clusters,photon2->GetCaloLabel(0),iclus2,iclus+1);
1061
1062       if(!cluster2)
1063       {
1064         printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - Second cluster not found\n");
1065         return;
1066       }
1067
1068       Float_t e1    = photon1->E();
1069       Float_t e2    = photon2->E();
1070
1071       //Select clusters with good time window difference
1072       Float_t tof1  = cluster1->GetTOF()*1e9;;
1073       Float_t tof2  = cluster2->GetTOF()*1e9;;
1074       Double_t t12diff = tof1-tof2;
1075       fhEPairDiffTime->Fill(e1+e2,    t12diff);
1076       if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
1077
1078       //Play with the MC stack if available
1079       if(IsDataMC()) HasPairSameMCMother(photon1, photon2, label, tag) ;
1080
1081       // Check the invariant mass for different selection on the local maxima
1082       // Name of AOD method TO BE FIXED
1083       Int_t nMaxima1 = photon1->GetFiducialArea();
1084       Int_t nMaxima2 = photon2->GetFiducialArea();
1085
1086       Double_t mass  = (mom1+mom2).M();
1087       Double_t epair = (mom1+mom2).E();
1088
1089       if(nMaxima1==nMaxima2)
1090       {
1091         if     (nMaxima1==1) fhMassPairLocMax[0]->Fill(epair,mass);
1092         else if(nMaxima1==2) fhMassPairLocMax[1]->Fill(epair,mass);
1093         else                 fhMassPairLocMax[2]->Fill(epair,mass);
1094       }
1095       else if(nMaxima1==1 || nMaxima2==1)
1096       {
1097         if  (nMaxima1==2 || nMaxima2==2) fhMassPairLocMax[3]->Fill(epair,mass);
1098         else                             fhMassPairLocMax[4]->Fill(epair,mass);
1099       }
1100       else
1101         fhMassPairLocMax[5]->Fill(epair,mass);
1102
1103       // combinations with SS axis cut and NLM cut
1104       if(nMaxima1 == 1 && cluster2->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1105       if(nMaxima2 == 1 && cluster1->GetM02() > 0.3) fhMassPairLocMax[6]->Fill(epair,mass);
1106       if(nMaxima1 >  1 && cluster2->GetM02() < 0.3 && cluster2->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1107       if(nMaxima2 >  1 && cluster1->GetM02() < 0.3 && cluster1->GetM02()> 0.1 ) fhMassPairLocMax[7]->Fill(epair,mass);
1108
1109
1110
1111       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1112       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1113       {
1114         if(GetDebug()>1)
1115           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());
1116
1117         //Fill some histograms about shower shape
1118         if(fFillSelectClHisto && clusters && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1119         {
1120           FillSelectedClusterHistograms(cluster1, nMaxima1, photon1->GetTag());
1121           FillSelectedClusterHistograms(cluster2, nMaxima2, photon2->GetTag());
1122         }
1123
1124         // Tag both photons as decay
1125         photon1->SetTagged(kTRUE);
1126         photon2->SetTagged(kTRUE);
1127
1128         fhPtDecay->Fill(photon1->Pt());
1129         fhEDecay ->Fill(photon1->E() );
1130
1131         fhPtDecay->Fill(photon2->Pt());
1132         fhEDecay ->Fill(photon2->E() );
1133
1134         //Create AOD for analysis
1135         mom = mom1+mom2;
1136
1137         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1138
1139         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1140         pi0.SetDetector(photon1->GetDetector());
1141
1142         // MC
1143         pi0.SetLabel(label);
1144         pi0.SetTag(tag);
1145
1146         //Set the indeces of the original caloclusters
1147         pi0.SetCaloLabel(photon1->GetCaloLabel(0), photon2->GetCaloLabel(0));
1148         //pi0.SetInputFileIndex(input);
1149
1150         AddAODParticle(pi0);
1151
1152       }//pi0
1153
1154     }//2n photon loop
1155
1156   }//1st photon loop
1157
1158   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeter() - End fill AODs \n");
1159
1160 }
1161
1162 //__________________________________________________
1163 void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
1164 {
1165   //Do analysis and fill aods
1166   //Search for the photon decay in calorimeters
1167   //Read photon list from AOD, produced in class AliAnaPhoton and AliGammaConversion
1168   //Check if 2 photons have the mass of the pi0.
1169
1170   TLorentzVector mom1;
1171   TLorentzVector mom2;
1172   TLorentzVector mom ;
1173   Int_t tag   = 0;
1174   Int_t label = 0;
1175   Int_t evtIndex = 0;
1176
1177   // Check calorimeter input
1178   if(!GetInputAODBranch()){
1179     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input calo photons in AOD branch with name < %s > , STOP\n",GetInputAODName().Data());
1180     abort();
1181   }
1182
1183   // Get the array with conversion photons
1184   TClonesArray * inputAODGammaConv = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaConvName);
1185   if(!inputAODGammaConv) {
1186
1187     inputAODGammaConv = (TClonesArray *) GetReader()->GetInputEvent()->FindListObject(fInputAODGammaConvName);
1188
1189     if(!inputAODGammaConv) {
1190       printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - No input gamma conversions in AOD branch with name < %s >\n",fInputAODGammaConvName.Data());
1191
1192       return;
1193     }
1194   }
1195
1196   //Get shower shape information of clusters
1197   TObjArray *clusters = 0;
1198   if     (fCalorimeter=="EMCAL") clusters = GetEMCALClusters();
1199   else if(fCalorimeter=="PHOS")  clusters = GetPHOSClusters() ;
1200
1201   Int_t nCTS  = inputAODGammaConv->GetEntriesFast();
1202   Int_t nCalo = GetInputAODBranch()->GetEntriesFast();
1203   if(nCTS<=0 || nCalo <=0) {
1204     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - nCalo %d, nCTS %d, cannot loop\n",nCalo,nCTS);
1205     return;
1206   }
1207
1208   if(GetDebug() > 1)
1209     printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - Number of conversion photons %d\n",nCTS);
1210
1211   // Do the loop, first calo, second CTS
1212   for(Int_t iphoton = 0; iphoton < GetInputAODBranch()->GetEntriesFast(); iphoton++){
1213     AliAODPWG4Particle * photon1 =  (AliAODPWG4Particle*) (GetInputAODBranch()->At(iphoton));
1214     mom1 = *(photon1->Momentum());
1215
1216     //Get original cluster, to recover some information
1217     Int_t iclus = -1;
1218     AliVCluster *cluster = FindCluster(clusters,photon1->GetCaloLabel(0),iclus);
1219
1220     for(Int_t jphoton = 0; jphoton < nCTS; jphoton++){
1221       AliAODPWG4Particle * photon2 =  (AliAODPWG4Particle*) (inputAODGammaConv->At(jphoton));
1222       if(GetMixedEvent())
1223         evtIndex = GetMixedEvent()->EventIndexForCaloCluster(photon2->GetCaloLabel(0)) ;
1224       if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
1225
1226       mom2 = *(photon2->Momentum());
1227
1228       Double_t mass  = (mom1+mom2).M();
1229       Double_t epair = (mom1+mom2).E();
1230
1231       Int_t nMaxima = photon1->GetFiducialArea();
1232       if     (nMaxima==1) fhMassPairLocMax[0]->Fill(epair,mass);
1233       else if(nMaxima==2) fhMassPairLocMax[1]->Fill(epair,mass);
1234       else                fhMassPairLocMax[2]->Fill(epair,mass);
1235
1236       //Play with the MC stack if available
1237       if(IsDataMC())
1238       {
1239         Int_t   label2 = photon2->GetLabel();
1240         if(label2 >= 0 )photon2->SetTag(GetMCAnalysisUtils()->CheckOrigin(label2, GetReader(), photon2->GetInputFileIndex()));
1241
1242         HasPairSameMCMother(photon1, photon2, label, tag) ;
1243       }
1244
1245       //Select good pair (good phi, pt cuts, aperture and invariant mass)
1246       if(GetNeutralMesonSelection()->SelectPair(mom1, mom2,fCalorimeter))
1247       {
1248         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());
1249
1250         //Fill some histograms about shower shape
1251         if(fFillSelectClHisto && cluster && GetReader()->GetDataType()!=AliCaloTrackReader::kMC)
1252         {
1253           FillSelectedClusterHistograms(cluster, nMaxima, photon1->GetTag());
1254         }
1255
1256         // Tag both photons as decay
1257         photon1->SetTagged(kTRUE);
1258         photon2->SetTagged(kTRUE);
1259
1260         fhPtDecay->Fill(photon1->Pt());
1261         fhEDecay ->Fill(photon1->E() );
1262
1263         //fhPtDecay->Fill(photon2->Pt());
1264         //fhEDecay ->Fill(photon2->E() );
1265
1266         //Create AOD for analysis
1267
1268         mom = mom1+mom2;
1269
1270         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
1271
1272         pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
1273         pi0.SetDetector(photon1->GetDetector());
1274
1275         // MC
1276         pi0.SetLabel(label);
1277         pi0.SetTag(tag);
1278
1279         //Set the indeces of the original tracks or caloclusters
1280         pi0.SetCaloLabel(photon1->GetCaloLabel(0), -1);
1281         pi0.SetTrackLabel(photon2->GetTrackLabel(0), photon2->GetTrackLabel(1));
1282         //pi0.SetInputFileIndex(input);
1283
1284         AddAODParticle(pi0);
1285
1286       }//pi0
1287     }//2n photon loop
1288
1289   }//1st photon loop
1290
1291   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS() - End fill AODs \n");
1292
1293 }
1294
1295
1296 //_________________________________________________
1297 void  AliAnaPi0EbE::MakeShowerShapeIdentification()
1298 {
1299   //Search for pi0 in fCalorimeter with shower shape analysis
1300
1301   TObjArray * pl = 0x0;
1302   //Select the Calorimeter of the photon
1303   if      (fCalorimeter == "PHOS" )
1304     pl = GetPHOSClusters();
1305   else if (fCalorimeter == "EMCAL")
1306     pl = GetEMCALClusters();
1307
1308   if(!pl) {
1309     Info("MakeShowerShapeIdentification","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1310     return;
1311   }
1312
1313   TLorentzVector mom ;
1314   for(Int_t icalo = 0; icalo < pl->GetEntriesFast(); icalo++)
1315   {
1316     AliVCluster * calo = (AliVCluster*) (pl->At(icalo));
1317
1318     Int_t evtIndex = 0 ;
1319     if (GetMixedEvent())
1320     {
1321       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1322     }
1323     if(TMath::Abs(GetVertex(evtIndex)[2]) > GetZvertexCut()) continue ;  //vertex cut
1324
1325     //Get Momentum vector,
1326     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1327       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1328     else{
1329       Double_t vertex[]={0,0,0};
1330       calo->GetMomentum(mom,vertex) ;
1331     }
1332
1333     //If too small or big pt, skip it
1334     if(mom.Pt() < GetMinPt() || mom.Pt() > GetMaxPt() ) continue ;
1335     //Check acceptance selection
1336     if(IsFiducialCutOn())
1337     {
1338       Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
1339       if(! in ) continue ;
1340     }
1341
1342     //Create AOD for analysis
1343     AliAODPWG4Particle aodpi0 = AliAODPWG4Particle(mom);
1344     aodpi0.SetLabel(calo->GetLabel());
1345
1346     //Set the indeces of the original caloclusters
1347     aodpi0.SetCaloLabel(calo->GetID(),-1);
1348     aodpi0.SetDetector(fCalorimeter);
1349     if(GetDebug() > 1)
1350       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());
1351
1352     //Check Distance to Bad channel, set bit.
1353     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1354     if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
1355     if(distBad < fMinDist) //In bad channel (PHOS cristal size 2.2x2.2 cm)
1356       continue ;
1357
1358     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: Bad channel cut passed %4.2f\n",distBad);
1359
1360     if     (distBad > fMinDist3) aodpi0.SetDistToBad(2) ;
1361     else if(distBad > fMinDist2) aodpi0.SetDistToBad(1) ;
1362     else                         aodpi0.SetDistToBad(0) ;
1363
1364     //Check PID
1365     //PID selection or bit setting
1366     if(IsCaloPIDOn())
1367     {
1368       //Skip matched clusters with tracks
1369       if(IsTrackMatched(calo, GetReader()->GetInputEvent())) continue ;
1370
1371       // Get most probable PID, 2 options check bayesian PID weights or redo PID
1372       // By default, redo PID
1373
1374       aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(calo));//PID recalculated
1375
1376       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
1377
1378       //If cluster does not pass pid, not pi0, skip it.
1379       if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
1380
1381     }
1382     else
1383     {
1384       //Set PID bits for later selection
1385       //GetPDG already called in SetPIDBits.
1386       GetCaloPID()->SetPIDBits(calo,&aodpi0, GetCaloUtils(), GetReader()->GetInputEvent());
1387       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");
1388     }
1389
1390     if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
1391
1392     //Play with the MC stack if available
1393     //Check origin of the candidates
1394     if(IsDataMC())
1395     {
1396       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1397          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1398         //aodpi0.SetInputFileIndex(input);
1399         Int_t tag       =0;
1400         tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabel(),GetReader(), aodpi0.GetInputFileIndex());
1401         //GetMCAnalysisUtils()->CheckMultipleOrigin(calo->GetLabels(),calo->GetNLabels(), GetReader(), aodpi0.GetInputFileIndex(), tag);
1402         aodpi0.SetTag(tag);
1403         if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Origin of candidate %d\n",aodpi0.GetTag());
1404       }
1405     }//Work with stack also
1406
1407     //Add AOD with pi0 object to aod branch
1408     AddAODParticle(aodpi0);
1409
1410   }//loop
1411
1412   if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - End fill AODs \n");
1413
1414 }
1415 //______________________________________________
1416 void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
1417 {
1418   //Do analysis and fill histograms
1419
1420   if(!GetOutputAODBranch())
1421   {
1422     printf("AliAnaPi0EbE::MakeAnalysisFillHistograms()  - No output pi0 in AOD branch with name < %s >,STOP \n",GetOutputAODName().Data());
1423     abort();
1424   }
1425   //Loop on stored AOD pi0
1426   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1427   if(GetDebug() > 0) printf("AliAnaPi0EbE::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1428
1429   for(Int_t iaod = 0; iaod < naod ; iaod++)
1430   {
1431
1432     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1433     Int_t pdg = pi0->GetIdentifiedParticleType();
1434
1435     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;
1436
1437     //Fill pi0 histograms
1438     Float_t ener  = pi0->E();
1439     Float_t pt    = pi0->Pt();
1440     Float_t phi   = pi0->Phi();
1441     if(phi < 0) phi+=TMath::TwoPi();
1442     Float_t eta = pi0->Eta();
1443
1444     fhPt     ->Fill(pt);
1445     fhE      ->Fill(ener);
1446
1447     fhEEta   ->Fill(ener,eta);
1448     fhEPhi   ->Fill(ener,phi);
1449     fhEtaPhi ->Fill(eta,phi);
1450
1451     if(IsDataMC())
1452     {
1453       if((GetReader()->GetDataType() == AliCaloTrackReader::kMC && fAnaType!=kSSCalo) ||
1454          GetReader()->GetDataType() != AliCaloTrackReader::kMC){
1455         if(GetMCAnalysisUtils()->CheckTagBit(pi0->GetTag(), AliMCAnalysisUtils::kMCPi0))
1456         {
1457           fhPtMC  ->Fill(pt);
1458           fhPhiMC ->Fill(pt,phi);
1459           fhEtaMC ->Fill(pt,eta);
1460         }
1461         else
1462         {
1463           fhPtMCNo  ->Fill(pt);
1464           fhPhiMCNo ->Fill(pt,phi);
1465           fhEtaMCNo ->Fill(pt,eta);
1466         }
1467       }
1468     }//Histograms with MC
1469
1470   }// aod loop
1471
1472 }
1473
1474 //__________________________________________________________________
1475 void AliAnaPi0EbE::Print(const Option_t * opt) const
1476 {
1477   //Print some relevant parameters set for the analysis
1478   if(! opt)
1479     return;
1480
1481   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1482   AliAnaCaloTrackCorrBaseClass::Print("");
1483   printf("Analysis Type = %d \n",  fAnaType) ;
1484   if(fAnaType == kSSCalo){
1485     printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
1486     printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1487     printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1488     printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1489   }
1490   printf("    \n") ;
1491
1492 }
1493
1494