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