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