PWG4/PartCorrDep/AliAnaPhoton.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 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 /* $Id: AliAnaPhoton.cxx 28688 2008-09-11 15:04:07Z gconesab $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for the photon identification.
  20 // Clusters from calorimeters are identified as photons
  21 // and kept in the AOD. Few histograms produced.
  22 // Produces input for other analysis classes like AliAnaPi0,
  23 // AliAnaParticleHadronCorrelation ...
  24 //
  25 // -- Author: Gustavo Conesa (LNF-INFN)
  26 //////////////////////////////////////////////////////////////////////////////
  27
  28
  29 // --- ROOT system ---
  30 #include <TH2F.h>
  31 #include <TH3D.h>
  32 #include <TClonesArray.h>
  33 #include <TObjString.h>
  34 //#include <Riostream.h>
  35 #include "TParticle.h"
  36 #include "TDatabasePDG.h"
  37
  38 // --- Analysis system ---
  39 #include "AliAnaPhoton.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliStack.h"
  42 #include "AliCaloPID.h"
  43 #include "AliMCAnalysisUtils.h"
  44 #include "AliFiducialCut.h"
  45 #include "AliVCluster.h"
  46 #include "AliAODMCParticle.h"
  47 #include "AliMixedEvent.h"
  48
  49
  50 ClassImp(AliAnaPhoton)
  51
  52 //____________________________________________________________________________
  53 AliAnaPhoton::AliAnaPhoton() :
  54     AliAnaPartCorrBaseClass(),    fCalorimeter(""),
  55     fMinDist(0.),                 fMinDist2(0.),                fMinDist3(0.),
  56     fRejectTrackMatch(0),         fTimeCutMin(-1),              fTimeCutMax(999999),
  57     fNCellsCut(0),                fFillSSHistograms(kFALSE),
  58     fCheckConversion(kFALSE),     fRemoveConvertedPair(kFALSE),
  59     fAddConvertedPairsToAOD(kFALSE),
  60     fMassCut(0),                  fConvAsymCut(1.),             fConvDEtaCut(2.),
  61     fConvDPhiMinCut(-1.),         fConvDPhiMaxCut(7.),
  62
  63     // Histograms
  64     fhNCellsE(0),                 fhEPhoton(0),                 fhPtPhoton(0),
  65     fhPhiPhoton(0),               fhEtaPhoton(0),
  66     fhEtaPhiPhoton(0),            fhEtaPhi05Photon(0),
  67
  68     // Conversion histograms
  69     fhPtPhotonConv(0),            fhEtaPhiPhotonConv(0),        fhEtaPhi05PhotonConv(0),
  70     fhConvDeltaEta(0),            fhConvDeltaPhi(0),            fhConvDeltaEtaPhi(0),
  71     fhConvAsym(0),                fhConvPt(0),
  72     fhConvDistEta(0),             fhConvDistEn(0),              fhConvDistMass(0),
  73     fhConvDistEtaCutEta(0),       fhConvDistEnCutEta(0),        fhConvDistMassCutEta(0),
  74     fhConvDistEtaCutMass(0),      fhConvDistEnCutMass(0),
  75     fhConvDistEtaCutAsy(0),       fhConvDistEnCutAsy(0),
  76
  77     //Shower shape histograms
  78     fhDispE(0),                   fhLam0E(0),                   fhLam1E(0),
  79     fhdDispE(0),                  fhdLam0E(0),                  fhdLam1E(0),
  80     fhDispETRD(0),                fhLam0ETRD(0),                fhLam1ETRD(0),
  81     fhdDispETRD(0),               fhdLam0ETRD(0),               fhdLam1ETRD(0),
  82
  83     fhNCellsLam0LowE(0),          fhNCellsLam1LowE(0),          fhNCellsDispLowE(0),
  84     fhNCellsLam0HighE(0),         fhNCellsLam1HighE(0),         fhNCellsDispHighE(0),
  85     fhNCellsdLam0LowE(0),         fhNCellsdLam1LowE(0),         fhNCellsdDispLowE(0),
  86     fhNCellsdLam0HighE(0),        fhNCellsdLam1HighE(0),        fhNCellsdDispHighE(0),
  87
  88     fhEtaLam0LowE(0),             fhPhiLam0LowE(0),
  89     fhEtaLam0HighE(0),            fhPhiLam0HighE(0),
  90     fhLam0DispLowE(0),            fhLam0DispHighE(0),
  91     fhLam1Lam0LowE(0),            fhLam1Lam0HighE(0),
  92     fhDispLam1LowE(0),            fhDispLam1HighE(0),
  93     fhEtadLam0LowE(0),            fhPhidLam0LowE(0),
  94     fhEtadLam0HighE(0),           fhPhidLam0HighE(0),
  95     fhdLam0dDispLowE(0),          fhdLam0dDispHighE(0),
  96     fhdLam1dLam0LowE(0),          fhdLam1dLam0HighE(0),
  97     fhdDispdLam1LowE(0),          fhdDispdLam1HighE(0),
  98
  99     //MC histograms
 100     fhDeltaE(0),                  fhDeltaPt(0),
 101     fhRatioE(0),                  fhRatioPt(0),
 102     fh2E(0),                      fh2Pt(0),
 103
 104     // Conversion MC histograms
 105     fhPtConversionTagged(0),           fhPtAntiNeutronTagged(0),
 106     fhPtAntiProtonTagged(0),           fhPtUnknownTagged(0),
 107     fhEtaPhiConversion(0),             fhEtaPhi05Conversion(0),
 108
 109     fhConvDeltaEtaMCConversion(0),     fhConvDeltaPhiMCConversion(0),  fhConvDeltaEtaPhiMCConversion(0),
 110     fhConvAsymMCConversion(0),         fhConvPtMCConversion(0),
 111     fhConvDispersionMCConversion(0),   fhConvM02MCConversion(0),
 112
 113     fhConvDeltaEtaMCAntiNeutron(0),    fhConvDeltaPhiMCAntiNeutron(0), fhConvDeltaEtaPhiMCAntiNeutron(0),
 114     fhConvAsymMCAntiNeutron(0),        fhConvPtMCAntiNeutron(0),
 115     fhConvDispersionMCAntiNeutron(0),  fhConvM02MCAntiNeutron(0),
 116     fhConvDeltaEtaMCAntiProton(0),     fhConvDeltaPhiMCAntiProton(0),  fhConvDeltaEtaPhiMCAntiProton(0),
 117     fhConvAsymMCAntiProton(0),         fhConvPtMCAntiProton(0),
 118     fhConvDispersionMCAntiProton(0),   fhConvM02MCAntiProton(0),
 119     fhConvDeltaEtaMCString(0),         fhConvDeltaPhiMCString(0),      fhConvDeltaEtaPhiMCString(0),
 120     fhConvAsymMCString(0),             fhConvPtMCString(0),
 121     fhConvDispersionMCString(0),       fhConvM02MCString(0),
 122     fhConvDistMCConversion(0),         fhConvDistMCConversionCuts(0),
 123
 124     // Photon SS MC histograms
 125     fhMCPhotonELambda0NoOverlap(0),    fhMCPhotonELambda0TwoOverlap(0),  fhMCPhotonELambda0NOverlap(0),
 126     fhMCPhotonEdLambda0NoOverlap(0),   fhMCPhotonEdLambda0TwoOverlap(0), fhMCPhotonEdLambda0NOverlap(0),
 127
 128     //Embedding
 129     fhEmbeddedSignalFractionEnergy(0),
 130     fhEmbedPhotonELambda0FullSignal(0),    fhEmbedPhotonEdLambda0FullSignal(0),
 131     fhEmbedPhotonELambda0MostlySignal(0),  fhEmbedPhotonEdLambda0MostlySignal(0),
 132     fhEmbedPhotonELambda0MostlyBkg(0),     fhEmbedPhotonEdLambda0MostlyBkg(0),
 133     fhEmbedPhotonELambda0FullBkg(0),       fhEmbedPhotonEdLambda0FullBkg(0),
 134     fhEmbedPi0ELambda0FullSignal(0),       fhEmbedPi0EdLambda0FullSignal(0),
 135     fhEmbedPi0ELambda0MostlySignal(0),     fhEmbedPi0EdLambda0MostlySignal(0),
 136     fhEmbedPi0ELambda0MostlyBkg(0),        fhEmbedPi0EdLambda0MostlyBkg(0),
 137     fhEmbedPi0ELambda0FullBkg(0),          fhEmbedPi0EdLambda0FullBkg(0)
 138 {
 139   //default ctor
 140
 141   for(Int_t i = 0; i < 12; i++){
 142     fhPtMC [i] = 0;
 143     fhMCE  [i] = 0;
 144     fhPhiMC[i] = 0;
 145     fhEtaMC[i] = 0;
 146   }
 147
 148   for(Int_t i = 0; i < 7; i++){
 149     fhPtPrimMC [i] = 0;
 150     fhEPrimMC  [i] = 0;
 151     fhPhiPrimMC[i] = 0;
 152     fhYPrimMC  [i] = 0;
 153
 154     fhPtPrimMCAcc [i] = 0;
 155     fhEPrimMCAcc  [i] = 0;
 156     fhPhiPrimMCAcc[i] = 0;
 157     fhYPrimMCAcc  [i] = 0;
 158   }
 159
 160   for(Int_t i = 0; i < 6; i++){
 161     fhMCELambda0    [i] = 0;
 162     fhMCELambda1    [i] = 0;
 163     fhMCEDispersion [i] = 0;
 164     fhMCEdLambda0   [i] = 0;
 165     fhMCEdLambda1   [i] = 0;
 166     fhMCEdDispersion[i] = 0;
 167   }
 168
 169   //Initialize parameters
 170   InitParameters();
 171
 172 }//____________________________________________________________________________
 173 AliAnaPhoton::~AliAnaPhoton()
 174 {
 175   //dtor
 176
 177 }
 178
 179 //__________________________________________________________________
 180 Bool_t  AliAnaPhoton::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
 181 {
 182   //Select clusters if they pass different cuts
 183   if(GetDebug() > 2)
 184     printf("AliAnaPhoton::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 185            GetReader()->GetEventNumber(),
 186            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 187
 188   //.......................................
 189   //If too small or big energy, skip it
 190   if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
 191   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
 192
 193   //.......................................
 194   // TOF cut, BE CAREFUL WITH THIS CUT
 195   Double_t tof = calo->GetTOF()*1e9;
 196   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
 197   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
 198
 199   //.......................................
 200   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
 201   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
 202
 203   //.......................................
 204   //Check acceptance selection
 205   if(IsFiducialCutOn()){
 206     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,fCalorimeter) ;
 207     if(! in ) return kFALSE ;
 208   }
 209   if(GetDebug() > 2) printf("Fiducial cut passed \n");
 210
 211   //.......................................
 212   //Skip matched clusters with tracks
 213   if(fRejectTrackMatch){
 214     if(IsTrackMatched(calo)) {
 215       if(GetDebug() > 2) printf("\t Reject track-matched clusters\n");
 216       return kFALSE ;
 217     }
 218     else
 219       if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 220   }// reject matched clusters
 221
 222   //.......................................
 223   //Check Distance to Bad channel, set bit.
 224   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 225   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 226   if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 227     return kFALSE ;
 228   }
 229   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 230   //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
 231
 232   if(GetDebug() > 0)
 233     printf("AliAnaPhoton::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 234            GetReader()->GetEventNumber(),
 235            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 236
 237   //All checks passed, cluster selected
 238   return kTRUE;
 239
 240 }
 241
 242 //_____________________________________________________________
 243 void AliAnaPhoton::FillAcceptanceHistograms(){
 244   //Fill acceptance histograms if MC data is available
 245
 246   if(GetReader()->ReadStack()){
 247     AliStack * stack = GetMCStack();
 248     if(stack){
 249       for(Int_t i=0 ; i<stack->GetNtrack(); i++){
 250         TParticle * prim = stack->Particle(i) ;
 251         Int_t pdg = prim->GetPdgCode();
 252         //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
 253         //                             prim->GetName(), prim->GetPdgCode());
 254
 255         if(pdg == 22){
 256
 257           // Get tag of this particle photon from fragmentation, decay, prompt ...
 258           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 259           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
 260             //A conversion photon from a hadron, skip this kind of photon
 261             //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
 262             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
 263             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
 264             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
 265             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
 266             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
 267             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
 268             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
 269             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
 270             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 271             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
 272             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
 273             //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
 274
 275             return;
 276           }
 277
 278           //Get photon kinematics
 279           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 280
 281           Double_t photonY   = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 282           Double_t photonE   = prim->Energy() ;
 283           Double_t photonPt  = prim->Pt() ;
 284           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
 285           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 286           Double_t photonEta = prim->Eta() ;
 287
 288
 289           //Check if photons hit the Calorimeter
 290           TLorentzVector lv;
 291           prim->Momentum(lv);
 292           Bool_t inacceptance = kFALSE;
 293           if     (fCalorimeter == "PHOS"){
 294             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 295               Int_t mod ;
 296               Double_t x,z ;
 297               if(GetPHOSGeometry()->ImpactOnEmc(prim,mod,z,x))
 298                 inacceptance = kTRUE;
 299               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 300             }
 301             else{
 302               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 303                 inacceptance = kTRUE ;
 304               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 305             }
 306           }
 307           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 308             if(GetEMCALGeometry()){
 309
 310               Int_t absID=0;
 311
 312               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 313
 314               if( absID >= 0)
 315                 inacceptance = kTRUE;
 316
 317               //                  if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
 318               //                    inacceptance = kTRUE;
 319               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 320             }
 321             else{
 322               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 323                 inacceptance = kTRUE ;
 324               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 325             }
 326           }       //In EMCAL
 327
 328           //Fill histograms
 329
 330           fhYPrimMC[mcPhoton]->Fill(photonPt, photonY) ;
 331           if(TMath::Abs(photonY) < 1.0)
 332           {
 333             fhEPrimMC  [mcPhoton]->Fill(photonE ) ;
 334             fhPtPrimMC [mcPhoton]->Fill(photonPt) ;
 335             fhPhiPrimMC[mcPhoton]->Fill(photonE , photonPhi) ;
 336             fhYPrimMC[mcPhoton]  ->Fill(photonE , photonEta) ;
 337           }
 338           if(inacceptance){
 339             fhEPrimMCAcc[mcPhoton]  ->Fill(photonE ) ;
 340             fhPtPrimMCAcc[mcPhoton] ->Fill(photonPt) ;
 341             fhPhiPrimMCAcc[mcPhoton]->Fill(photonE , photonPhi) ;
 342             fhYPrimMCAcc[mcPhoton]  ->Fill(photonE , photonY) ;
 343           }//Accepted
 344
 345           //Origin of photon
 346           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
 347           {
 348             fhYPrimMC[mcPrompt]->Fill(photonPt, photonY) ;
 349             if(TMath::Abs(photonY) < 1.0){
 350               fhEPrimMC  [mcPrompt]->Fill(photonE ) ;
 351               fhPtPrimMC [mcPrompt]->Fill(photonPt) ;
 352               fhPhiPrimMC[mcPrompt]->Fill(photonE , photonPhi) ;
 353               fhYPrimMC[mcPrompt]  ->Fill(photonE , photonEta) ;
 354             }
 355             if(inacceptance){
 356               fhEPrimMCAcc[mcPrompt]  ->Fill(photonE ) ;
 357               fhPtPrimMCAcc[mcPrompt] ->Fill(photonPt) ;
 358               fhPhiPrimMCAcc[mcPrompt]->Fill(photonE , photonPhi) ;
 359               fhYPrimMCAcc[mcPrompt]  ->Fill(photonE , photonY) ;
 360             }//Accepted
 361           }
 362           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
 363           {
 364             fhYPrimMC[mcFragmentation]->Fill(photonPt, photonY) ;
 365             if(TMath::Abs(photonY) < 1.0){
 366               fhEPrimMC  [mcFragmentation]->Fill(photonE ) ;
 367               fhPtPrimMC [mcFragmentation]->Fill(photonPt) ;
 368               fhPhiPrimMC[mcFragmentation]->Fill(photonE , photonPhi) ;
 369               fhYPrimMC[mcFragmentation]  ->Fill(photonE , photonEta) ;
 370             }
 371             if(inacceptance){
 372               fhEPrimMCAcc[mcFragmentation]  ->Fill(photonE ) ;
 373               fhPtPrimMCAcc[mcFragmentation] ->Fill(photonPt) ;
 374               fhPhiPrimMCAcc[mcFragmentation]->Fill(photonE , photonPhi) ;
 375               fhYPrimMCAcc[mcFragmentation]  ->Fill(photonE , photonY) ;
 376             }//Accepted
 377           }
 378           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
 379           {
 380             fhYPrimMC[mcISR]->Fill(photonPt, photonY) ;
 381             if(TMath::Abs(photonY) < 1.0){
 382               fhEPrimMC  [mcISR]->Fill(photonE ) ;
 383               fhPtPrimMC [mcISR]->Fill(photonPt) ;
 384               fhPhiPrimMC[mcISR]->Fill(photonE , photonPhi) ;
 385               fhYPrimMC[mcISR]->Fill(photonE , photonEta) ;
 386             }
 387             if(inacceptance){
 388               fhEPrimMCAcc[mcISR]  ->Fill(photonE ) ;
 389               fhPtPrimMCAcc[mcISR] ->Fill(photonPt) ;
 390               fhPhiPrimMCAcc[mcISR]->Fill(photonE , photonPhi) ;
 391               fhYPrimMCAcc[mcISR]  ->Fill(photonE , photonY) ;
 392             }//Accepted
 393           }
 394           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 395           {
 396             fhYPrimMC[mcPi0Decay]->Fill(photonPt, photonY) ;
 397             if(TMath::Abs(photonY) < 1.0){
 398               fhEPrimMC  [mcPi0Decay]->Fill(photonE ) ;
 399               fhPtPrimMC [mcPi0Decay]->Fill(photonPt) ;
 400               fhPhiPrimMC[mcPi0Decay]->Fill(photonE , photonPhi) ;
 401               fhYPrimMC[mcPi0Decay]  ->Fill(photonE , photonEta) ;
 402             }
 403             if(inacceptance){
 404               fhEPrimMCAcc[mcPi0Decay]  ->Fill(photonE ) ;
 405               fhPtPrimMCAcc[mcPi0Decay] ->Fill(photonPt) ;
 406               fhPhiPrimMCAcc[mcPi0Decay]->Fill(photonE , photonPhi) ;
 407               fhYPrimMCAcc[mcPi0Decay]  ->Fill(photonE , photonY) ;
 408             }//Accepted
 409           }
 410           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 411                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 412           {
 413             fhYPrimMC[mcOtherDecay]->Fill(photonPt, photonY) ;
 414             if(TMath::Abs(photonY) < 1.0){
 415               fhEPrimMC  [mcOtherDecay]->Fill(photonE ) ;
 416               fhPtPrimMC [mcOtherDecay]->Fill(photonPt) ;
 417               fhPhiPrimMC[mcOtherDecay]->Fill(photonE , photonPhi) ;
 418               fhYPrimMC[mcOtherDecay]  ->Fill(photonE , photonEta) ;
 419             }
 420             if(inacceptance){
 421               fhEPrimMCAcc[mcOtherDecay]  ->Fill(photonE ) ;
 422               fhPtPrimMCAcc[mcOtherDecay] ->Fill(photonPt) ;
 423               fhPhiPrimMCAcc[mcOtherDecay]->Fill(photonE , photonPhi) ;
 424               fhYPrimMCAcc[mcOtherDecay]  ->Fill(photonE , photonY) ;
 425             }//Accepted
 426           }
 427           else
 428           {
 429             fhYPrimMC[mcOther]->Fill(photonPt, photonY) ;
 430             if(TMath::Abs(photonY) < 1.0){
 431               fhEPrimMC  [mcOther]->Fill(photonE ) ;
 432               fhPtPrimMC [mcOther]->Fill(photonPt) ;
 433               fhPhiPrimMC[mcOther]->Fill(photonE , photonPhi) ;
 434               fhYPrimMC[mcOther]  ->Fill(photonE , photonEta) ;
 435             }
 436             if(inacceptance){
 437               fhEPrimMCAcc[mcOther]  ->Fill(photonE ) ;
 438               fhPtPrimMCAcc[mcOther] ->Fill(photonPt) ;
 439               fhPhiPrimMCAcc[mcOther]->Fill(photonE , photonPhi) ;
 440               fhYPrimMCAcc[mcOther]  ->Fill(photonE , photonY) ;
 441             }//Accepted
 442           }//Other origin
 443         }// Primary photon
 444       }//loop on primaries
 445     }//stack exists and data is MC
 446   }//read stack
 447   else if(GetReader()->ReadAODMCParticles()){
 448     TClonesArray * mcparticles = GetReader()->GetAODMCParticles(0);
 449     if(mcparticles){
 450       Int_t nprim = mcparticles->GetEntriesFast();
 451
 452       for(Int_t i=0; i < nprim; i++)
 453       {
 454         AliAODMCParticle * prim = (AliAODMCParticle *) mcparticles->At(i);
 455
 456         Int_t pdg = prim->GetPdgCode();
 457
 458         if(pdg == 22){
 459
 460           // Get tag of this particle photon from fragmentation, decay, prompt ...
 461           Int_t tag = GetMCAnalysisUtils()->CheckOrigin(i,GetReader(), 0);
 462           if(!GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton)){
 463             //A conversion photon from a hadron, skip this kind of photon
 464 //            printf("AliAnaPhoton::FillAcceptanceHistograms() - not a photon, weird!: tag %d, conv %d, pi0 %d, hadron %d, electron %d, unk %d, muon %d,pion %d, proton %d, neutron %d, kaon %d, antiproton %d, antineutron %d\n",tag,
 465 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion),
 466 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0),
 467 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOther),
 468 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron),
 469 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCUnknown),
 470 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCMuon),
 471 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPion),
 472 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCProton),
 473 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron),
 474 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCKaon),
 475 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton),
 476 //                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron));
 477
 478             return;
 479           }
 480
 481           //Get photon kinematics
 482           if(prim->E() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 483
 484           Double_t photonY  = 0.5*TMath::Log((prim->E()-prim->Pz())/(prim->E()+prim->Pz())) ;
 485           Double_t photonE   = prim->E() ;
 486           Double_t photonPt  = prim->Pt() ;
 487           Double_t photonPhi = TMath::RadToDeg()*prim->Phi() ;
 488           if(photonPhi < 0) photonPhi+=TMath::TwoPi();
 489           Double_t photonEta = prim->Eta() ;
 490
 491           //Check if photons hit the Calorimeter
 492           TLorentzVector lv;
 493           lv.SetPxPyPzE(prim->Px(),prim->Py(),prim->Pz(),prim->E());
 494           Bool_t inacceptance = kFALSE;
 495           if     (fCalorimeter == "PHOS"){
 496             if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 497               Int_t mod ;
 498               Double_t x,z ;
 499               Double_t vtx[]={prim->Xv(),prim->Yv(),prim->Zv()};
 500               if(GetPHOSGeometry()->ImpactOnEmc(vtx, prim->Theta(),prim->Phi(),mod,z,x))
 501                 inacceptance = kTRUE;
 502               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 503             }
 504             else{
 505               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 506                 inacceptance = kTRUE ;
 507               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 508             }
 509           }
 510           else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 511             if(GetEMCALGeometry()){
 512
 513               Int_t absID=0;
 514
 515               GetEMCALGeometry()->GetAbsCellIdFromEtaPhi(prim->Eta(),prim->Phi(),absID);
 516
 517               if( absID >= 0)
 518                 inacceptance = kTRUE;
 519
 520               if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 521             }
 522             else{
 523               if(GetFiducialCut()->IsInFiducialCut(lv,fCalorimeter))
 524                 inacceptance = kTRUE ;
 525               if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 526             }
 527           }       //In EMCAL
 528
 529           //Fill histograms
 530
 531           fhYPrimMC[mcPhoton]->Fill(photonPt, photonY) ;
 532           if(TMath::Abs(photonY) < 1.0)
 533           {
 534             fhEPrimMC  [mcPhoton]->Fill(photonE ) ;
 535             fhPtPrimMC [mcPhoton]->Fill(photonPt) ;
 536             fhPhiPrimMC[mcPhoton]->Fill(photonE , photonPhi) ;
 537             fhYPrimMC[mcPhoton]->Fill(photonE , photonEta) ;
 538           }
 539           if(inacceptance){
 540             fhEPrimMCAcc[mcPhoton]  ->Fill(photonE ) ;
 541             fhPtPrimMCAcc[mcPhoton] ->Fill(photonPt) ;
 542             fhPhiPrimMCAcc[mcPhoton]->Fill(photonE , photonPhi) ;
 543             fhYPrimMCAcc[mcPhoton]  ->Fill(photonE , photonY) ;
 544           }//Accepted
 545
 546
 547           //Origin of photon
 548           if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt))
 549           {
 550             fhYPrimMC[mcPrompt]->Fill(photonPt, photonY) ;
 551             if(TMath::Abs(photonY) < 1.0){
 552               fhEPrimMC  [mcPrompt]->Fill(photonE ) ;
 553               fhPtPrimMC [mcPrompt]->Fill(photonPt) ;
 554               fhPhiPrimMC[mcPrompt]->Fill(photonE , photonPhi) ;
 555               fhYPrimMC[mcPrompt]->Fill(photonE , photonEta) ;
 556             }
 557             if(inacceptance){
 558               fhEPrimMCAcc[mcPrompt]  ->Fill(photonE ) ;
 559               fhPtPrimMCAcc[mcPrompt] ->Fill(photonPt) ;
 560               fhPhiPrimMCAcc[mcPrompt]->Fill(photonE , photonPhi) ;
 561               fhYPrimMCAcc[mcPrompt]  ->Fill(photonE , photonY) ;
 562             }//Accepted
 563           }
 564           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
 565           {
 566             fhYPrimMC[mcFragmentation]->Fill(photonPt, photonY) ;
 567             if(TMath::Abs(photonY) < 1.0){
 568               fhEPrimMC  [mcFragmentation]->Fill(photonE ) ;
 569               fhPtPrimMC [mcFragmentation]->Fill(photonPt) ;
 570               fhPhiPrimMC[mcFragmentation]->Fill(photonE , photonPhi) ;
 571               fhYPrimMC[mcFragmentation]->Fill(photonE , photonEta) ;
 572             }
 573             if(inacceptance){
 574               fhEPrimMCAcc[mcFragmentation]  ->Fill(photonE ) ;
 575               fhPtPrimMCAcc[mcFragmentation] ->Fill(photonPt) ;
 576               fhPhiPrimMCAcc[mcFragmentation]->Fill(photonE , photonPhi) ;
 577               fhYPrimMCAcc[mcFragmentation]  ->Fill(photonE , photonY) ;
 578             }//Accepted
 579           }
 580           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
 581           {
 582             fhYPrimMC[mcISR]->Fill(photonPt, photonY) ;
 583             if(TMath::Abs(photonY) < 1.0){
 584               fhEPrimMC  [mcISR]->Fill(photonE ) ;
 585               fhPtPrimMC [mcISR]->Fill(photonPt) ;
 586               fhPhiPrimMC[mcISR]->Fill(photonE , photonPhi) ;
 587               fhYPrimMC[mcISR]->Fill(photonE , photonEta) ;
 588             }
 589             if(inacceptance){
 590               fhEPrimMCAcc[mcISR]  ->Fill(photonE ) ;
 591               fhPtPrimMCAcc[mcISR] ->Fill(photonPt) ;
 592               fhPhiPrimMCAcc[mcISR]->Fill(photonE , photonPhi) ;
 593               fhYPrimMCAcc[mcISR]  ->Fill(photonE , photonY) ;
 594             }//Accepted
 595           }
 596           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay))
 597           {
 598             fhYPrimMC[mcPi0Decay]->Fill(photonPt, photonY) ;
 599             if(TMath::Abs(photonY) < 1.0){
 600               fhEPrimMC  [mcPi0Decay]->Fill(photonE ) ;
 601               fhPtPrimMC [mcPi0Decay]->Fill(photonPt) ;
 602               fhPhiPrimMC[mcPi0Decay]->Fill(photonE , photonPhi) ;
 603               fhYPrimMC[mcPi0Decay]->Fill(photonE , photonEta) ;
 604             }
 605             if(inacceptance){
 606               fhEPrimMCAcc[mcPi0Decay]  ->Fill(photonE ) ;
 607               fhPtPrimMCAcc[mcPi0Decay] ->Fill(photonPt) ;
 608               fhPhiPrimMCAcc[mcPi0Decay]->Fill(photonE , photonPhi) ;
 609               fhYPrimMCAcc[mcPi0Decay]  ->Fill(photonE , photonY) ;
 610             }//Accepted
 611           }
 612           else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
 613                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
 614           {
 615             fhYPrimMC[mcOtherDecay]->Fill(photonPt, photonY) ;
 616             if(TMath::Abs(photonY) < 1.0){
 617               fhEPrimMC  [mcOtherDecay]->Fill(photonE ) ;
 618               fhPtPrimMC [mcOtherDecay]->Fill(photonPt) ;
 619               fhPhiPrimMC[mcOtherDecay]->Fill(photonE , photonPhi) ;
 620               fhYPrimMC[mcOtherDecay]->Fill(photonE , photonEta) ;
 621             }
 622             if(inacceptance){
 623               fhEPrimMCAcc[mcOtherDecay]  ->Fill(photonE ) ;
 624               fhPtPrimMCAcc[mcOtherDecay] ->Fill(photonPt) ;
 625               fhPhiPrimMCAcc[mcOtherDecay]->Fill(photonE , photonPhi) ;
 626               fhYPrimMCAcc[mcOtherDecay]  ->Fill(photonE , photonY) ;
 627             }//Accepted
 628           }
 629           else
 630           {
 631             fhYPrimMC[mcOther]->Fill(photonPt, photonY) ;
 632             if(TMath::Abs(photonY) < 1.0){
 633               fhEPrimMC  [mcOther]->Fill(photonE ) ;
 634               fhPtPrimMC [mcOther]->Fill(photonPt) ;
 635               fhPhiPrimMC[mcOther]->Fill(photonE , photonPhi) ;
 636               fhYPrimMC[mcOther]->Fill(photonE , photonEta) ;
 637             }
 638             if(inacceptance){
 639               fhEPrimMCAcc[mcOther]  ->Fill(photonE ) ;
 640               fhPtPrimMCAcc[mcOther] ->Fill(photonPt) ;
 641               fhPhiPrimMCAcc[mcOther]->Fill(photonE , photonPhi) ;
 642               fhYPrimMCAcc[mcOther]  ->Fill(photonE , photonY) ;
 643             }//Accepted
 644           }//Other origin
 645         }// Primary photon
 646       }//loop on primaries
 647
 648     }//mc array exists and data is MC
 649   }     // read AOD MC
 650 }
 651
 652 //__________________________________________________________________
 653 void  AliAnaPhoton::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag){
 654
 655   //Fill cluster Shower Shape histograms
 656
 657   if(!fFillSSHistograms || GetMixedEvent()) return;
 658
 659   Float_t energy  = cluster->E();
 660   Int_t   ncells  = cluster->GetNCells();
 661   Int_t   ncells2 = ncells*ncells;
 662   Float_t lambda0 = cluster->GetM02();
 663   Float_t lambda1 = cluster->GetM20();
 664   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
 665
 666   TLorentzVector mom;
 667   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 668     cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
 669   else{
 670     Double_t vertex[]={0,0,0};
 671     cluster->GetMomentum(mom,vertex) ;
 672   }
 673
 674   Float_t eta = mom.Eta();
 675   Float_t phi = mom.Phi();
 676   if(phi < 0) phi+=TMath::TwoPi();
 677
 678   fhLam0E ->Fill(energy,lambda0);
 679   fhLam1E ->Fill(energy,lambda1);
 680   fhDispE ->Fill(energy,disp);
 681   fhdLam0E->Fill(energy,lambda0/ncells2);
 682   fhdLam1E->Fill(energy,lambda1/ncells2);
 683   fhdDispE->Fill(energy,disp/ncells2);
 684
 685   if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
 686     fhLam0ETRD->Fill(energy,lambda0);
 687     fhLam1ETRD->Fill(energy,lambda1);
 688     fhDispETRD->Fill(energy,disp);
 689     fhdLam0ETRD->Fill(energy,lambda0/ncells2);
 690     fhdLam1ETRD->Fill(energy,lambda1/ncells2);
 691     fhdDispETRD->Fill(energy,disp/ncells2);
 692   }
 693
 694   if(energy < 2){
 695     fhNCellsLam0LowE ->Fill(ncells,lambda0);
 696     fhNCellsLam1LowE ->Fill(ncells,lambda1);
 697     fhNCellsDispLowE ->Fill(ncells,disp);
 698     fhNCellsdLam0LowE->Fill(ncells,lambda0/ncells2);
 699     fhNCellsdLam1LowE->Fill(ncells,lambda1/ncells2);
 700     fhNCellsdDispLowE->Fill(ncells,disp/ncells2);
 701
 702     fhLam1Lam0LowE  ->Fill(lambda1,lambda0);
 703     fhLam0DispLowE  ->Fill(lambda0,disp);
 704     fhDispLam1LowE  ->Fill(disp,lambda1);
 705     fhEtaLam0LowE   ->Fill(eta,lambda0);
 706     fhPhiLam0LowE   ->Fill(phi,lambda0);
 707
 708     fhdLam1dLam0LowE->Fill(lambda1/ncells2,lambda0/ncells2);
 709     fhdLam0dDispLowE->Fill(lambda0/ncells2,disp/ncells2);
 710     fhdDispdLam1LowE->Fill(disp/ncells2,lambda1/ncells2);
 711     fhEtadLam0LowE  ->Fill(eta,lambda0/ncells2);
 712     fhPhidLam0LowE  ->Fill(phi,lambda0/ncells2);
 713   }
 714   else {
 715     fhNCellsLam0HighE ->Fill(ncells,lambda0);
 716     fhNCellsLam1HighE ->Fill(ncells,lambda1);
 717     fhNCellsDispHighE ->Fill(ncells,disp);
 718     fhNCellsdLam0HighE->Fill(ncells,lambda0/ncells2);
 719     fhNCellsdLam1HighE->Fill(ncells,lambda1/ncells2);
 720     fhNCellsdDispHighE->Fill(ncells,disp/ncells2);
 721
 722     fhLam1Lam0HighE  ->Fill(lambda1,lambda0);
 723     fhLam0DispHighE  ->Fill(lambda0,disp);
 724     fhDispLam1HighE  ->Fill(disp,lambda1);
 725     fhEtaLam0HighE   ->Fill(eta, lambda0);
 726     fhPhiLam0HighE   ->Fill(phi, lambda0);
 727
 728     fhdLam1dLam0HighE->Fill(lambda1/ncells2,lambda0/ncells2);
 729     fhdLam0dDispHighE->Fill(lambda0/ncells2,disp/ncells2);
 730     fhdDispdLam1HighE->Fill(disp/ncells2,lambda1/ncells2);
 731     fhEtadLam0HighE  ->Fill(eta, lambda0/ncells2);
 732     fhPhidLam0HighE  ->Fill(phi, lambda0/ncells2);
 733   }
 734
 735   if(IsDataMC()){
 736
 737
 738     //Fill histograms to check shape of embedded clusters
 739     Float_t fraction = 0;
 740     if(GetReader()->IsEmbeddedClusterSelectionOn()){
 741       AliVCaloCells * cells = 0;
 742       if(fCalorimeter == "EMCAL") cells = GetReader()->GetEMCALCells();
 743       else                        cells = GetReader()->GetPHOSCells(); //only working for EMCAL, but in case in future it works
 744
 745       Float_t clusterE = 0; // recalculate in case corrections applied.
 746       Float_t cellE    = 0;
 747       for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
 748         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
 749         clusterE+=cellE;
 750         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
 751       }
 752
 753       //Fraction of total energy due to the embedded signal
 754       fraction/=clusterE;
 755
 756       if(GetDebug() > 1 ) printf("AliAnaPhoton::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
 757
 758       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
 759
 760     }  // embedded fraction
 761
 762     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 763        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
 764        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
 765        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
 766       fhMCELambda0[mcssPhoton]    ->Fill(energy, lambda0);
 767       fhMCEdLambda0[mcssPhoton]   ->Fill(energy, lambda0/ncells2);
 768       fhMCELambda1[mcssPhoton]    ->Fill(energy, lambda1);
 769       fhMCEdLambda1[mcssPhoton]   ->Fill(energy, lambda1/ncells2);
 770       fhMCEDispersion[mcssPhoton] ->Fill(energy, disp);
 771       fhMCEdDispersion[mcssPhoton]->Fill(energy, disp/ncells2);
 772
 773
 774       if(!GetReader()->IsEmbeddedClusterSelectionOn()){
 775         //Check particle overlaps in cluster
 776
 777         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
 778         Int_t ancPDG = 0, ancStatus = -1;
 779         TLorentzVector momentum; TVector3 prodVertex;
 780         Int_t ancLabel = 0;
 781         Int_t noverlaps = 1;
 782         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
 783           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
 784           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
 785         }
 786
 787         if(noverlaps == 1){
 788           fhMCPhotonELambda0NoOverlap  ->Fill(energy, lambda0);
 789           fhMCPhotonEdLambda0NoOverlap ->Fill(energy, lambda0/ncells2);
 790         }
 791         else if(noverlaps == 2){
 792           fhMCPhotonELambda0TwoOverlap ->Fill(energy, lambda0);
 793           fhMCPhotonEdLambda0TwoOverlap->Fill(energy, lambda0/ncells2);
 794         }
 795         else if(noverlaps > 2){
 796           fhMCPhotonELambda0NOverlap   ->Fill(energy, lambda0);
 797           fhMCPhotonEdLambda0NOverlap  ->Fill(energy, lambda0/ncells2);
 798         }
 799         else {
 800           printf("AliAnaPhoton::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
 801         }
 802       }//No embedding
 803
 804       //Fill histograms to check shape of embedded clusters
 805       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 806
 807         if     (fraction > 0.9)
 808         {
 809           fhEmbedPhotonELambda0FullSignal   ->Fill(energy, lambda0);
 810           fhEmbedPhotonEdLambda0FullSignal  ->Fill(energy, lambda0/ncells2);
 811         }
 812         else if(fraction > 0.5)
 813         {
 814           fhEmbedPhotonELambda0MostlySignal ->Fill(energy, lambda0);
 815           fhEmbedPhotonEdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
 816         }
 817         else if(fraction > 0.1)
 818         {
 819           fhEmbedPhotonELambda0MostlyBkg    ->Fill(energy, lambda0);
 820           fhEmbedPhotonEdLambda0MostlyBkg   ->Fill(energy, lambda0/ncells2);
 821         }
 822         else
 823         {
 824           fhEmbedPhotonELambda0FullBkg      ->Fill(energy, lambda0);
 825           fhEmbedPhotonEdLambda0FullBkg     ->Fill(energy, lambda0/ncells2);
 826         }
 827       } // embedded
 828
 829     }//photon   no conversion
 830     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)){
 831       fhMCELambda0[mcssElectron]    ->Fill(energy, lambda0);
 832       fhMCEdLambda0[mcssElectron]   ->Fill(energy, lambda0/ncells2);
 833       fhMCELambda1[mcssElectron]    ->Fill(energy, lambda1);
 834       fhMCEdLambda1[mcssElectron]   ->Fill(energy, lambda1/ncells2);
 835       fhMCEDispersion[mcssElectron] ->Fill(energy, disp);
 836       fhMCEdDispersion[mcssElectron]->Fill(energy, disp/ncells2);
 837     }//electron
 838     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 839               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
 840       fhMCELambda0[mcssConversion]    ->Fill(energy, lambda0);
 841       fhMCEdLambda0[mcssConversion]   ->Fill(energy, lambda0/ncells2);
 842       fhMCELambda1[mcssConversion]    ->Fill(energy, lambda1);
 843       fhMCEdLambda1[mcssConversion]   ->Fill(energy, lambda1/ncells2);
 844       fhMCEDispersion[mcssConversion] ->Fill(energy, disp);
 845       fhMCEdDispersion[mcssConversion]->Fill(energy, disp/ncells2);
 846
 847     }//conversion photon
 848     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  ){
 849       fhMCELambda0[mcssPi0]    ->Fill(energy, lambda0);
 850       fhMCEdLambda0[mcssPi0]   ->Fill(energy, lambda0/ncells2);
 851       fhMCELambda1[mcssPi0]    ->Fill(energy, lambda1);
 852       fhMCEdLambda1[mcssPi0]   ->Fill(energy, lambda1/ncells2);
 853       fhMCEDispersion[mcssPi0] ->Fill(energy, disp);
 854       fhMCEdDispersion[mcssPi0]->Fill(energy, disp/ncells2);
 855
 856       //Fill histograms to check shape of embedded clusters
 857       if(GetReader()->IsEmbeddedClusterSelectionOn()){
 858
 859         if     (fraction > 0.9)
 860         {
 861           fhEmbedPi0ELambda0FullSignal   ->Fill(energy, lambda0);
 862           fhEmbedPi0EdLambda0FullSignal  ->Fill(energy, lambda0/ncells2);
 863         }
 864         else if(fraction > 0.5)
 865         {
 866           fhEmbedPi0ELambda0MostlySignal ->Fill(energy, lambda0);
 867           fhEmbedPi0EdLambda0MostlySignal->Fill(energy, lambda0/ncells2);
 868         }
 869         else if(fraction > 0.1)
 870         {
 871           fhEmbedPi0ELambda0MostlyBkg    ->Fill(energy, lambda0);
 872           fhEmbedPi0EdLambda0MostlyBkg   ->Fill(energy, lambda0/ncells2);
 873         }
 874         else
 875         {
 876           fhEmbedPi0ELambda0FullBkg      ->Fill(energy, lambda0);
 877           fhEmbedPi0EdLambda0FullBkg     ->Fill(energy, lambda0/ncells2);
 878         }
 879       } // embedded
 880
 881     }//pi0
 882     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ){
 883       fhMCELambda0[mcssEta]    ->Fill(energy, lambda0);
 884       fhMCEdLambda0[mcssEta]   ->Fill(energy, lambda0/ncells2);
 885       fhMCELambda1[mcssEta]    ->Fill(energy, lambda1);
 886       fhMCEdLambda1[mcssEta]   ->Fill(energy, lambda1/ncells2);
 887       fhMCEDispersion[mcssEta] ->Fill(energy, disp);
 888       fhMCEdDispersion[mcssEta]->Fill(energy, disp/ncells2);
 889     }//eta
 890     else {
 891       fhMCELambda0[mcssOther]    ->Fill(energy, lambda0);
 892       fhMCEdLambda0[mcssOther]   ->Fill(energy, lambda0/ncells2);
 893       fhMCELambda1[mcssOther]    ->Fill(energy, lambda1);
 894       fhMCEdLambda1[mcssOther]   ->Fill(energy, lambda1/ncells2);
 895       fhMCEDispersion[mcssOther] ->Fill(energy, disp);
 896       fhMCEdDispersion[mcssOther]->Fill(energy, disp/ncells2);
 897     }//other particles
 898
 899   }//MC data
 900
 901 }
 902
 903 //________________________________________________________________________
 904 TObjString *  AliAnaPhoton::GetAnalysisCuts()
 905 {
 906   //Save parameters used for analysis
 907   TString parList ; //this will be list of parameters used for this analysis.
 908   const Int_t buffersize = 255;
 909   char onePar[buffersize] ;
 910
 911   snprintf(onePar,buffersize,"--- AliAnaPhoton ---\n") ;
 912   parList+=onePar ;
 913   snprintf(onePar,buffersize,"Calorimeter: %s\n",fCalorimeter.Data()) ;
 914   parList+=onePar ;
 915   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 916   parList+=onePar ;
 917   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 918   parList+=onePar ;
 919   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 920   parList+=onePar ;
 921   snprintf(onePar,buffersize,"fRejectTrackMatch: %d\n",fRejectTrackMatch) ;
 922   parList+=onePar ;
 923   snprintf(onePar,buffersize,"Conversion Selection: fConvAsymCut %1.2f, fConvDEtaCut %1.2f fConvDPhiCut (%1.2f,%1.2f)\n",
 924            fConvAsymCut, fConvDEtaCut, fConvDPhiMinCut, fConvDPhiMaxCut) ;
 925   parList+=onePar ;
 926
 927   //Get parameters set in base class.
 928   parList += GetBaseParametersList() ;
 929
 930   //Get parameters set in PID class.
 931   parList += GetCaloPID()->GetPIDParametersList() ;
 932
 933   //Get parameters set in FiducialCut class (not available yet)
 934   //parlist += GetFidCut()->GetFidCutParametersList()
 935
 936   return new TObjString(parList) ;
 937 }
 938
 939 //________________________________________________________________________
 940 TList *  AliAnaPhoton::GetCreateOutputObjects()
 941 {
 942   // Create histograms to be saved in output file and
 943   // store them in outputContainer
 944   TList * outputContainer = new TList() ;
 945   outputContainer->SetName("PhotonHistos") ;
 946
 947   Int_t nptbins  = GetHistoPtBins();  Float_t ptmax  = GetHistoPtMax();  Float_t ptmin  = GetHistoPtMin();
 948   Int_t nphibins = GetHistoPhiBins(); Float_t phimax = GetHistoPhiMax(); Float_t phimin = GetHistoPhiMin();
 949   Int_t netabins = GetHistoEtaBins(); Float_t etamax = GetHistoEtaMax(); Float_t etamin = GetHistoEtaMin();
 950   Int_t ssbins   = GetHistoShowerShapeBins();  Float_t ssmax = GetHistoShowerShapeMax();  Float_t ssmin = GetHistoShowerShapeMin();
 951
 952   fhNCellsE  = new TH2F ("hNCellsE","# of cells in cluster vs E of clusters", 100,0, 20, 20,0,20);
 953   fhNCellsE->SetXTitle("E (GeV)");
 954   fhNCellsE->SetYTitle("# of cells in cluster");
 955   outputContainer->Add(fhNCellsE);
 956
 957   fhEPhoton  = new TH1F("hEPhoton","Number of #gamma over calorimeter vs energy",nptbins,ptmin,ptmax);
 958   fhEPhoton->SetYTitle("N");
 959   fhEPhoton->SetXTitle("E_{#gamma}(GeV)");
 960   outputContainer->Add(fhEPhoton) ;
 961
 962   fhPtPhoton  = new TH1F("hPtPhoton","Number of #gamma over calorimeter vs p_{T}",nptbins,ptmin,ptmax);
 963   fhPtPhoton->SetYTitle("N");
 964   fhPtPhoton->SetXTitle("p_{T #gamma}(GeV/c)");
 965   outputContainer->Add(fhPtPhoton) ;
 966
 967   fhPhiPhoton  = new TH2F
 968     ("hPhiPhoton","#phi_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 969   fhPhiPhoton->SetYTitle("#phi (rad)");
 970   fhPhiPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 971   outputContainer->Add(fhPhiPhoton) ;
 972
 973   fhEtaPhoton  = new TH2F
 974     ("hEtaPhoton","#eta_{#gamma} vs p_{T}",nptbins,ptmin,ptmax,netabins,etamin,etamax);
 975   fhEtaPhoton->SetYTitle("#eta");
 976   fhEtaPhoton->SetXTitle("p_{T #gamma} (GeV/c)");
 977   outputContainer->Add(fhEtaPhoton) ;
 978
 979   fhEtaPhiPhoton  = new TH2F
 980   ("hEtaPhiPhoton","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
 981   fhEtaPhiPhoton->SetYTitle("#phi (rad)");
 982   fhEtaPhiPhoton->SetXTitle("#eta");
 983   outputContainer->Add(fhEtaPhiPhoton) ;
 984   if(GetMinPt() < 0.5){
 985     fhEtaPhi05Photon  = new TH2F
 986     ("hEtaPhi05Photon","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
 987     fhEtaPhi05Photon->SetYTitle("#phi (rad)");
 988     fhEtaPhi05Photon->SetXTitle("#eta");
 989     outputContainer->Add(fhEtaPhi05Photon) ;
 990   }
 991
 992   //Conversion
 993   if(fCheckConversion){
 994
 995     fhEtaPhiConversion  = new TH2F
 996     ("hEtaPhiConversion","#eta vs #phi for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
 997     fhEtaPhiConversion->SetYTitle("#phi (rad)");
 998     fhEtaPhiConversion->SetXTitle("#eta");
 999     outputContainer->Add(fhEtaPhiConversion) ;
1000     if(GetMinPt() < 0.5){
1001       fhEtaPhi05Conversion  = new TH2F
1002       ("hEtaPhi05Conversion","#eta vs #phi, E > 0.5, for converted clusters",netabins,etamin,etamax,nphibins,phimin,phimax);
1003       fhEtaPhi05Conversion->SetYTitle("#phi (rad)");
1004       fhEtaPhi05Conversion->SetXTitle("#eta");
1005       outputContainer->Add(fhEtaPhi05Conversion) ;
1006     }
1007
1008     fhPtPhotonConv  = new TH1F("hPtPhotonConv","Number of #gamma over calorimeter, conversion",nptbins,ptmin,ptmax);
1009     fhPtPhotonConv->SetYTitle("N");
1010     fhPtPhotonConv->SetXTitle("p_{T #gamma}(GeV/c)");
1011     outputContainer->Add(fhPtPhotonConv) ;
1012
1013     fhEtaPhiPhotonConv  = new TH2F
1014     ("hEtaPhiPhotonConv","#eta vs #phi",netabins,etamin,etamax,nphibins,phimin,phimax);
1015     fhEtaPhiPhotonConv->SetYTitle("#phi (rad)");
1016     fhEtaPhiPhotonConv->SetXTitle("#eta");
1017     outputContainer->Add(fhEtaPhiPhotonConv) ;
1018     if(GetMinPt() < 0.5){
1019       fhEtaPhi05PhotonConv  = new TH2F
1020       ("hEtaPhi05PhotonConv","#eta vs #phi, E > 0.5",netabins,etamin,etamax,nphibins,phimin,phimax);
1021       fhEtaPhi05PhotonConv->SetYTitle("#phi (rad)");
1022       fhEtaPhi05PhotonConv->SetXTitle("#eta");
1023       outputContainer->Add(fhEtaPhi05PhotonConv) ;
1024     }
1025
1026     fhConvDeltaEta  = new TH2F
1027     ("hConvDeltaEta","#Delta #eta of selected conversion pairs",100,0,fMassCut,netabins*2,-0.5,0.5);
1028     fhConvDeltaEta->SetYTitle("#Delta #eta");
1029     fhConvDeltaEta->SetXTitle("Pair Mass (GeV/c^2)");
1030     outputContainer->Add(fhConvDeltaEta) ;
1031
1032     fhConvDeltaPhi  = new TH2F
1033     ("hConvDeltaPhi","#Delta #phi of selected conversion pairs",100,0,fMassCut,nphibins*2,-0.5,0.5);
1034     fhConvDeltaPhi->SetYTitle("#Delta #phi");
1035     fhConvDeltaPhi->SetXTitle("Pair Mass (GeV/c^2)");
1036     outputContainer->Add(fhConvDeltaPhi) ;
1037
1038     fhConvDeltaEtaPhi  = new TH2F
1039     ("hConvDeltaEtaPhi","#Delta #eta vs #Delta #phi of selected conversion pairs",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1040     fhConvDeltaEtaPhi->SetYTitle("#Delta #phi");
1041     fhConvDeltaEtaPhi->SetXTitle("#Delta #eta");
1042     outputContainer->Add(fhConvDeltaEtaPhi) ;
1043
1044     fhConvAsym  = new TH2F
1045     ("hConvAsym","Asymmetry of selected conversion pairs",100,0,fMassCut,100,0,1);
1046     fhConvAsym->SetYTitle("Asymmetry");
1047     fhConvAsym->SetXTitle("Pair Mass (GeV/c^2)");
1048     outputContainer->Add(fhConvAsym) ;
1049
1050     fhConvPt  = new TH2F
1051     ("hConvPt","p_{T} of selected conversion pairs",100,0,fMassCut,100,0.,10.);
1052     fhConvPt->SetYTitle("Pair p_{T} (GeV/c)");
1053     fhConvPt->SetXTitle("Pair Mass (GeV/c^2)");
1054     outputContainer->Add(fhConvPt) ;
1055
1056     fhConvDistEta  = new TH2F
1057     ("hConvDistEta","distance to conversion vertex",100,-0.7,0.7,100,0.,5.);
1058     fhConvDistEta->SetXTitle("#eta");
1059     fhConvDistEta->SetYTitle(" distance (m)");
1060     outputContainer->Add(fhConvDistEta) ;
1061
1062     fhConvDistEn  = new TH2F
1063     ("hConvDistEn","distance to conversion vertex",nptbins,ptmin,ptmax,100,0.,5.);
1064     fhConvDistEn->SetXTitle("E (GeV)");
1065     fhConvDistEn->SetYTitle(" distance (m)");
1066     outputContainer->Add(fhConvDistEn) ;
1067
1068     fhConvDistMass  = new TH2F
1069     ("hConvDistMass","distance to conversion vertex",100,0,fMassCut,100,0.,5.);
1070     fhConvDistMass->SetXTitle("m (GeV/c^2)");
1071     fhConvDistMass->SetYTitle(" distance (m)");
1072     outputContainer->Add(fhConvDistMass) ;
1073
1074     fhConvDistEtaCutEta  = new TH2F
1075     ("hConvDistEtaCutEta","distance to conversion vertex, dEta < 0.05",100,-0.7,0.7,100,0.,5.);
1076     fhConvDistEtaCutEta->SetXTitle("#eta");
1077     fhConvDistEtaCutEta->SetYTitle(" distance (m)");
1078     outputContainer->Add(fhConvDistEtaCutEta) ;
1079
1080     fhConvDistEnCutEta  = new TH2F
1081     ("hConvDistEnCutEta","distance to conversion vertex, dEta < 0.05",nptbins,ptmin,ptmax,100,0.,5.);
1082     fhConvDistEnCutEta->SetXTitle("E (GeV)");
1083     fhConvDistEnCutEta->SetYTitle(" distance (m)");
1084     outputContainer->Add(fhConvDistEnCutEta) ;
1085
1086     fhConvDistMassCutEta  = new TH2F
1087     ("hConvDistMassCutEta","distance to conversion vertex, dEta < 0.05",100,0,fMassCut,100,0.,5.);
1088     fhConvDistMassCutEta->SetXTitle("m (GeV/c^2)");
1089     fhConvDistMassCutEta->SetYTitle(" distance (m)");
1090     outputContainer->Add(fhConvDistMassCutEta) ;
1091
1092     fhConvDistEtaCutMass  = new TH2F
1093     ("hConvDistEtaCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",100,-0.7,0.7,100,0.,5.);
1094     fhConvDistEtaCutMass->SetXTitle("#eta");
1095     fhConvDistEtaCutMass->SetYTitle(" distance (m)");
1096     outputContainer->Add(fhConvDistEtaCutMass) ;
1097
1098     fhConvDistEnCutMass  = new TH2F
1099     ("hConvDistEnCutMass","distance to conversion vertex, dEta < 0.05, m < 10 MeV",nptbins,ptmin,ptmax,100,0.,5.);
1100     fhConvDistEnCutMass->SetXTitle("E (GeV)");
1101     fhConvDistEnCutMass->SetYTitle(" distance (m)");
1102     outputContainer->Add(fhConvDistEnCutMass) ;
1103
1104     fhConvDistEtaCutAsy  = new TH2F
1105     ("hConvDistEtaCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",100,-0.7,0.7,100,0.,5.);
1106     fhConvDistEtaCutAsy->SetXTitle("#eta");
1107     fhConvDistEtaCutAsy->SetYTitle(" distance (m)");
1108     outputContainer->Add(fhConvDistEtaCutAsy) ;
1109
1110     fhConvDistEnCutAsy  = new TH2F
1111     ("hConvDistEnCutAsy","distance to conversion vertex, dEta < 0.05, m < 10 MeV, A < 0.1",nptbins,ptmin,ptmax,100,0.,5.);
1112     fhConvDistEnCutAsy->SetXTitle("E (GeV)");
1113     fhConvDistEnCutAsy->SetYTitle(" distance (m)");
1114     outputContainer->Add(fhConvDistEnCutAsy) ;
1115
1116   } // check conversion
1117
1118
1119   //Shower shape
1120   if(fFillSSHistograms){
1121
1122     fhLam0E  = new TH2F ("hLam0E","#lambda_{0}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1123     fhLam0E->SetYTitle("#lambda_{0}^{2}");
1124     fhLam0E->SetXTitle("E (GeV)");
1125     outputContainer->Add(fhLam0E);
1126
1127     fhLam1E  = new TH2F ("hLam1E","#lambda_{1}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1128     fhLam1E->SetYTitle("#lambda_{1}^{2}");
1129     fhLam1E->SetXTitle("E (GeV)");
1130     outputContainer->Add(fhLam1E);
1131
1132     fhDispE  = new TH2F ("hDispE"," dispersion^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1133     fhDispE->SetYTitle("D^{2}");
1134     fhDispE->SetXTitle("E (GeV) ");
1135     outputContainer->Add(fhDispE);
1136
1137     fhdLam0E  = new TH2F ("hdLam0E","#lambda_{0}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1138     fhdLam0E->SetYTitle("d#lambda_{0}^{2}");
1139     fhdLam0E->SetXTitle("E (GeV)");
1140     outputContainer->Add(fhdLam0E);
1141
1142     fhdLam1E  = new TH2F ("hdLam1E","#lambda_{1}^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1143     fhdLam1E->SetYTitle("d#lambda_{1}^{2}");
1144     fhdLam1E->SetXTitle("E (GeV)");
1145     outputContainer->Add(fhdLam1E);
1146
1147     fhdDispE  = new TH2F ("hdDispE"," dispersion^{2}/N_{cells}^{2} vs E", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1148     fhdDispE->SetYTitle("dD^{2}");
1149     fhdDispE->SetXTitle("E (GeV) ");
1150     outputContainer->Add(fhdDispE);
1151
1152
1153     if(fCalorimeter == "EMCAL"){
1154       fhLam0ETRD  = new TH2F ("hLam0ETRD","#lambda_{0}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1155       fhLam0ETRD->SetYTitle("#lambda_{0}^{2}");
1156       fhLam0ETRD->SetXTitle("E (GeV)");
1157       outputContainer->Add(fhLam0ETRD);
1158
1159       fhLam1ETRD  = new TH2F ("hLam1ETRD","#lambda_{1}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1160       fhLam1ETRD->SetYTitle("#lambda_{1}^{2}");
1161       fhLam1ETRD->SetXTitle("E (GeV)");
1162       outputContainer->Add(fhLam1ETRD);
1163
1164       fhDispETRD  = new TH2F ("hDispETRD"," dispersion^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1165       fhDispETRD->SetYTitle("Dispersion^{2}");
1166       fhDispETRD->SetXTitle("E (GeV) ");
1167       outputContainer->Add(fhDispETRD);
1168
1169       fhdLam0ETRD  = new TH2F ("hdLam0ETRD","#lambda_{0}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1170       fhdLam0ETRD->SetYTitle("d#lambda_{0}^{2}");
1171       fhdLam0ETRD->SetXTitle("E (GeV)");
1172       outputContainer->Add(fhdLam0ETRD);
1173
1174       fhdLam1ETRD  = new TH2F ("hdLam1ETRD","#lambda_{1}^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1175       fhdLam1ETRD->SetYTitle("d#lambda_{1}^{2}");
1176       fhdLam1ETRD->SetXTitle("E (GeV)");
1177       outputContainer->Add(fhdLam1ETRD);
1178
1179       fhdDispETRD  = new TH2F ("hdDispETRD"," dispersion^{2}/N_{cells}^{2} vs E, EMCAL SM covered by TRD", nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1180       fhdDispETRD->SetYTitle("dD^{2}");
1181       fhdDispETRD->SetXTitle("E (GeV) ");
1182       outputContainer->Add(fhdDispETRD);
1183
1184     }
1185
1186     fhNCellsLam0LowE  = new TH2F ("hNCellsLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1187     fhNCellsLam0LowE->SetXTitle("N_{cells}");
1188     fhNCellsLam0LowE->SetYTitle("#lambda_{0}^{2}");
1189     outputContainer->Add(fhNCellsLam0LowE);
1190
1191     fhNCellsLam0HighE  = new TH2F ("hNCellsLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1192     fhNCellsLam0HighE->SetXTitle("N_{cells}");
1193     fhNCellsLam0HighE->SetYTitle("#lambda_{0}^{2}");
1194     outputContainer->Add(fhNCellsLam0HighE);
1195
1196     fhNCellsLam1LowE  = new TH2F ("hNCellsLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1197     fhNCellsLam1LowE->SetXTitle("N_{cells}");
1198     fhNCellsLam1LowE->SetYTitle("#lambda_{0}^{2}");
1199     outputContainer->Add(fhNCellsLam1LowE);
1200
1201     fhNCellsLam1HighE  = new TH2F ("hNCellsLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1202     fhNCellsLam1HighE->SetXTitle("N_{cells}");
1203     fhNCellsLam1HighE->SetYTitle("#lambda_{0}^{2}");
1204     outputContainer->Add(fhNCellsLam1HighE);
1205
1206     fhNCellsDispLowE  = new TH2F ("hNCellsDispLowE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1207     fhNCellsDispLowE->SetXTitle("N_{cells}");
1208     fhNCellsDispLowE->SetYTitle("D^{2}");
1209     outputContainer->Add(fhNCellsDispLowE);
1210
1211     fhNCellsDispHighE  = new TH2F ("hNCellsDispHighE","N_{cells} in cluster vs dispersion^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax);
1212     fhNCellsDispHighE->SetXTitle("N_{cells}");
1213     fhNCellsDispHighE->SetYTitle("D^{2}");
1214     outputContainer->Add(fhNCellsDispHighE);
1215
1216
1217     fhNCellsdLam0LowE  = new TH2F ("hNCellsdLam0LowE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1218     fhNCellsdLam0LowE->SetXTitle("N_{cells}");
1219     fhNCellsdLam0LowE->SetYTitle("#lambda_{0}^{2}");
1220     outputContainer->Add(fhNCellsdLam0LowE);
1221
1222     fhNCellsdLam0HighE  = new TH2F ("hNCellsdLam0HighE","N_{cells} in cluster vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1223     fhNCellsdLam0HighE->SetXTitle("N_{cells}");
1224     fhNCellsdLam0HighE->SetYTitle("#lambda_{0}^{2}");
1225     outputContainer->Add(fhNCellsdLam0HighE);
1226
1227     fhNCellsdLam1LowE  = new TH2F ("hNCellsdLam1LowE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1228     fhNCellsdLam1LowE->SetXTitle("N_{cells}");
1229     fhNCellsdLam1LowE->SetYTitle("#lambda_{0}^{2}");
1230     outputContainer->Add(fhNCellsdLam1LowE);
1231
1232     fhNCellsdLam1HighE  = new TH2F ("hNCellsdLam1HighE","N_{cells} in cluster vs #lambda_{1}^{2}/N_{cells}^{2}, E > 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1233     fhNCellsdLam1HighE->SetXTitle("N_{cells}");
1234     fhNCellsdLam1HighE->SetYTitle("#lambda_{0}^{2}");
1235     outputContainer->Add(fhNCellsdLam1HighE);
1236
1237     fhNCellsdDispLowE  = new TH2F ("hNCellsdDispLowE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1238     fhNCellsdDispLowE->SetXTitle("N_{cells}");
1239     fhNCellsdDispLowE->SetYTitle("D^{2}");
1240     outputContainer->Add(fhNCellsdDispLowE);
1241
1242     fhNCellsdDispHighE  = new TH2F ("hNCellsdDispHighE","N_{cells} in cluster vs dispersion^{2}/N_{cells}^{2}, E < 2 GeV", 20,0, 20, ssbins,ssmin,ssmax/50);
1243     fhNCellsdDispHighE->SetXTitle("N_{cells}");
1244     fhNCellsdDispHighE->SetYTitle("D^{2}");
1245     outputContainer->Add(fhNCellsdDispHighE);
1246
1247
1248     fhEtaLam0LowE  = new TH2F ("hEtaLam0LowE","#eta vs #lambda_{0}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1249     fhEtaLam0LowE->SetYTitle("#lambda_{0}^{2}");
1250     fhEtaLam0LowE->SetXTitle("#eta");
1251     outputContainer->Add(fhEtaLam0LowE);
1252
1253     fhPhiLam0LowE  = new TH2F ("hPhiLam0LowE","#phi vs #lambda_{0}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1254     fhPhiLam0LowE->SetYTitle("#lambda_{0}^{2}");
1255     fhPhiLam0LowE->SetXTitle("#phi");
1256     outputContainer->Add(fhPhiLam0LowE);
1257
1258     fhEtaLam0HighE  = new TH2F ("hEtaLam0HighE","#eta vs #lambda_{0}^{2}, E > 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax);
1259     fhEtaLam0HighE->SetYTitle("#lambda_{0}^{2}");
1260     fhEtaLam0HighE->SetXTitle("#eta");
1261     outputContainer->Add(fhEtaLam0HighE);
1262
1263     fhPhiLam0HighE  = new TH2F ("hPhiLam0HighE","#phi vs #lambda_{0}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax);
1264     fhPhiLam0HighE->SetYTitle("#lambda_{0}^{2}");
1265     fhPhiLam0HighE->SetXTitle("#phi");
1266     outputContainer->Add(fhPhiLam0HighE);
1267
1268     fhLam1Lam0LowE  = new TH2F ("hLam1Lam0LowE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1269     fhLam1Lam0LowE->SetYTitle("#lambda_{0}^{2}");
1270     fhLam1Lam0LowE->SetXTitle("#lambda_{1}^{2}");
1271     outputContainer->Add(fhLam1Lam0LowE);
1272
1273     fhLam1Lam0HighE  = new TH2F ("hLam1Lam0HighE","#lambda_{0}^{2} vs #lambda_{1}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1274     fhLam1Lam0HighE->SetYTitle("#lambda_{0}^{2}");
1275     fhLam1Lam0HighE->SetXTitle("#lambda_{1}^{2}");
1276     outputContainer->Add(fhLam1Lam0HighE);
1277
1278     fhLam0DispLowE  = new TH2F ("hLam0DispLowE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1279     fhLam0DispLowE->SetXTitle("#lambda_{0}^{2}");
1280     fhLam0DispLowE->SetYTitle("D^{2}");
1281     outputContainer->Add(fhLam0DispLowE);
1282
1283     fhLam0DispHighE  = new TH2F ("hLam0DispHighE","#lambda_{0}^{2} vs dispersion^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1284     fhLam0DispHighE->SetXTitle("#lambda_{0}^{2}");
1285     fhLam0DispHighE->SetYTitle("D^{2}");
1286     outputContainer->Add(fhLam0DispHighE);
1287
1288     fhDispLam1LowE  = new TH2F ("hDispLam1LowE","Dispersion^{2} vs #lambda_{1}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1289     fhDispLam1LowE->SetXTitle("D^{2}");
1290     fhDispLam1LowE->SetYTitle("#lambda_{1}^{2}");
1291     outputContainer->Add(fhDispLam1LowE);
1292
1293     fhDispLam1HighE  = new TH2F ("hDispLam1HighE","Dispersion^{2} vs #lambda_{1^{2}} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax, ssbins,ssmin,ssmax);
1294     fhDispLam1HighE->SetXTitle("D^{2}");
1295     fhDispLam1HighE->SetYTitle("#lambda_{1}^{2}");
1296     outputContainer->Add(fhDispLam1HighE);
1297
1298
1299
1300     fhEtadLam0LowE  = new TH2F ("hEtadLam0LowE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
1301     fhEtadLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1302     fhEtadLam0LowE->SetXTitle("#eta");
1303     outputContainer->Add(fhEtadLam0LowE);
1304
1305     fhPhidLam0LowE  = new TH2F ("hPhidLam0LowE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E < 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
1306     fhPhidLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1307     fhPhidLam0LowE->SetXTitle("#phi");
1308     outputContainer->Add(fhPhidLam0LowE);
1309
1310     fhEtadLam0HighE  = new TH2F ("hEtadLam0HighE","#eta vs #lambda_{0}^{2}/N_{cells}^{2}", netabins,etamin,etamax, ssbins,ssmin,ssmax/50);
1311     fhEtadLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1312     fhEtadLam0HighE->SetXTitle("#eta");
1313     outputContainer->Add(fhEtadLam0HighE);
1314
1315     fhPhidLam0HighE  = new TH2F ("hPhidLam0HighE","#phi vs #lambda_{0}^{2}/N_{cells}^{2}, E > 2 GeV", nphibins,phimin,phimax, ssbins,ssmin,ssmax/50);
1316     fhPhidLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1317     fhPhidLam0HighE->SetXTitle("#phi");
1318     outputContainer->Add(fhPhidLam0HighE);
1319
1320     fhdLam1dLam0LowE  = new TH2F ("hdLam1dLam0LowE","#lambda_{0}^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1321     fhdLam1dLam0LowE->SetYTitle("d#lambda_{0}^{2}");
1322     fhdLam1dLam0LowE->SetXTitle("d#lambda_{1}^{2}");
1323     outputContainer->Add(fhdLam1dLam0LowE);
1324
1325     fhdLam1dLam0HighE  = new TH2F ("hdLam1dLam0HighE","#lambda_{0}^{2}/N_{cells}^{2}vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1326     fhdLam1dLam0HighE->SetYTitle("d#lambda_{0}^{2}");
1327     fhdLam1dLam0HighE->SetXTitle("d#lambda_{1}^{2}");
1328     outputContainer->Add(fhdLam1dLam0HighE);
1329
1330     fhdLam0dDispLowE  = new TH2F ("hdLam0dDispLowE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1331     fhdLam0dDispLowE->SetXTitle("d#lambda_{0}^{2}");
1332     fhdLam0dDispLowE->SetYTitle("dD^{2}");
1333     outputContainer->Add(fhdLam0dDispLowE);
1334
1335     fhdLam0dDispHighE  = new TH2F ("hdLam0dDispHighE","#lambda_{0}^{2}/N_{cells}^{2} vs dispersion^{2}/N_{cells}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1336     fhdLam0dDispHighE->SetXTitle("d#lambda_{0}^{2}");
1337     fhdLam0dDispHighE->SetYTitle("dD^{2}");
1338     outputContainer->Add(fhdLam0dDispHighE);
1339
1340     fhdDispdLam1LowE  = new TH2F ("hdDispddLam1LowE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1}^{2}/N_{cells}^{2} in cluster of E < 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1341     fhdDispdLam1LowE->SetXTitle("dD^{2}");
1342     fhdDispdLam1LowE->SetYTitle("d#lambda_{1}^{2}");
1343     outputContainer->Add(fhdDispdLam1LowE);
1344
1345     fhdDispdLam1HighE  = new TH2F ("hdDispdLam1HighE","Dispersion^{2}/N_{cells}^{2} vs #lambda_{1^{2}}/N_{cells}^{2} in cluster of E > 2 GeV",  ssbins,ssmin,ssmax/50, ssbins,ssmin,ssmax/50);
1346     fhdDispdLam1HighE->SetXTitle("dD^{2}");
1347     fhdDispdLam1HighE->SetYTitle("d#lambda_{1}^{2}");
1348     outputContainer->Add(fhdDispdLam1HighE);
1349
1350
1351   } // Shower shape
1352
1353
1354   if(IsDataMC()){
1355     fhDeltaE  = new TH1F ("hDeltaE","MC - Reco E ", 200,-50,50);
1356     fhDeltaE->SetXTitle("#Delta E (GeV)");
1357     outputContainer->Add(fhDeltaE);
1358
1359     fhDeltaPt  = new TH1F ("hDeltaPt","MC - Reco p_{T} ", 200,-50,50);
1360     fhDeltaPt->SetXTitle("#Delta p_{T} (GeV/c)");
1361     outputContainer->Add(fhDeltaPt);
1362
1363     fhRatioE  = new TH1F ("hRatioE","Reco/MC E ", 200,0,2);
1364     fhRatioE->SetXTitle("E_{reco}/E_{gen}");
1365     outputContainer->Add(fhRatioE);
1366
1367     fhRatioPt  = new TH1F ("hRatioPt","Reco/MC p_{T} ", 200,0,2);
1368     fhRatioPt->SetXTitle("p_{T, reco}/p_{T, gen}");
1369     outputContainer->Add(fhRatioPt);
1370
1371     fh2E  = new TH2F ("h2E","E distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1372     fh2E->SetXTitle("E_{rec} (GeV)");
1373     fh2E->SetYTitle("E_{gen} (GeV)");
1374     outputContainer->Add(fh2E);
1375
1376     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
1377     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
1378     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
1379     outputContainer->Add(fh2Pt);
1380
1381     TString ptype[] = { "#gamma","#gamma_{prompt}","#gamma_{fragmentation}","#gamma_{ISR}",
1382                         "#gamma_{#pi decay}","#gamma_{decay}","hadron?",
1383                         "#pi^{0}","#gamma->e^{#pm}","e^{#pm}",
1384                         "Anti-N","Anti-P","String"                                } ;
1385
1386     TString pname[] = { "Photon","PhotonPrompt","PhotonFragmentation","PhotonISR",
1387                         "PhotonPi0Decay", "PhotonOtherDecay", "Hadron",
1388                         "Pi0","Conversion","Electron",
1389                         "AntiNeutron","AntiProton","String"                       } ;
1390
1391     for(Int_t i = 0; i < 13; i++){
1392
1393       fhMCE[i]  = new TH1F(Form("hE_MC%s",pname[i].Data()),
1394                                 Form("cluster from %s : E ",ptype[i].Data()),
1395                                 nptbins,ptmin,ptmax);
1396       fhMCE[i]->SetXTitle("E (GeV)");
1397       outputContainer->Add(fhMCE[i]) ;
1398
1399       fhPtMC[i]  = new TH1F(Form("hPt_MC%s",pname[i].Data()),
1400                            Form("cluster from %s : p_{T} ",ptype[i].Data()),
1401                            nptbins,ptmin,ptmax);
1402       fhPtMC[i]->SetXTitle("p_{T} (GeV/c)");
1403       outputContainer->Add(fhPtMC[i]) ;
1404
1405       fhEtaMC[i]  = new TH2F(Form("hEta_MC%s",pname[i].Data()),
1406                            Form("cluster from %s : #eta ",ptype[i].Data()),
1407                            nptbins,ptmin,ptmax,netabins,etamin,etamax);
1408       fhEtaMC[i]->SetYTitle("#eta");
1409       fhEtaMC[i]->SetXTitle("E (GeV)");
1410       outputContainer->Add(fhEtaMC[i]) ;
1411
1412       fhPhiMC[i]  = new TH2F(Form("hPhi_MC%s",pname[i].Data()),
1413                            Form("cluster from %s : #phi ",ptype[i].Data()),
1414                            nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1415       fhPhiMC[i]->SetYTitle("#phi (rad)");
1416       fhPhiMC[i]->SetXTitle("E (GeV)");
1417       outputContainer->Add(fhPhiMC[i]) ;
1418
1419     }
1420
1421     for(Int_t i = 0; i < 7; i++){
1422       fhEPrimMC[i]  = new TH1F(Form("hEPrim_MC%s",pname[i].Data()),
1423                            Form("primary photon %s : E ",ptype[i].Data()),
1424                            nptbins,ptmin,ptmax);
1425       fhEPrimMC[i]->SetXTitle("E (GeV)");
1426       outputContainer->Add(fhEPrimMC[i]) ;
1427
1428       fhPtPrimMC[i]  = new TH1F(Form("hPtPrim_MC%s",pname[i].Data()),
1429                             Form("primary photon %s : p_{T} ",ptype[i].Data()),
1430                             nptbins,ptmin,ptmax);
1431       fhPtPrimMC[i]->SetXTitle("p_{T} (GeV/c)");
1432       outputContainer->Add(fhPtPrimMC[i]) ;
1433
1434       fhYPrimMC[i]  = new TH2F(Form("hYPrim_MC%s",pname[i].Data()),
1435                              Form("primary photon %s : Rapidity ",ptype[i].Data()),
1436                              nptbins,ptmin,ptmax,800,-8,8);
1437       fhYPrimMC[i]->SetYTitle("Rapidity");
1438       fhYPrimMC[i]->SetXTitle("E (GeV)");
1439       outputContainer->Add(fhYPrimMC[i]) ;
1440
1441       fhPhiPrimMC[i]  = new TH2F(Form("hPhiPrim_MC%s",pname[i].Data()),
1442                              Form("primary photon %s : #phi ",ptype[i].Data()),
1443                              nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1444       fhPhiPrimMC[i]->SetYTitle("#phi (rad)");
1445       fhPhiPrimMC[i]->SetXTitle("E (GeV)");
1446       outputContainer->Add(fhPhiPrimMC[i]) ;
1447
1448
1449       fhEPrimMCAcc[i]  = new TH1F(Form("hEPrimAcc_MC%s",pname[i].Data()),
1450                                Form("primary photon %s in acceptance: E ",ptype[i].Data()),
1451                                nptbins,ptmin,ptmax);
1452       fhEPrimMCAcc[i]->SetXTitle("E (GeV)");
1453       outputContainer->Add(fhEPrimMCAcc[i]) ;
1454
1455       fhPtPrimMCAcc[i]  = new TH1F(Form("hPtPrimAcc_MC%s",pname[i].Data()),
1456                                 Form("primary photon %s in acceptance: p_{T} ",ptype[i].Data()),
1457                                 nptbins,ptmin,ptmax);
1458       fhPtPrimMCAcc[i]->SetXTitle("p_{T} (GeV/c)");
1459       outputContainer->Add(fhPtPrimMCAcc[i]) ;
1460
1461       fhYPrimMCAcc[i]  = new TH2F(Form("hYPrimAcc_MC%s",pname[i].Data()),
1462                                  Form("primary photon %s in acceptance: Rapidity ",ptype[i].Data()),
1463                                  nptbins,ptmin,ptmax,100,-1,1);
1464       fhYPrimMCAcc[i]->SetYTitle("Rapidity");
1465       fhYPrimMCAcc[i]->SetXTitle("E (GeV)");
1466       outputContainer->Add(fhYPrimMCAcc[i]) ;
1467
1468       fhPhiPrimMCAcc[i]  = new TH2F(Form("hPhiPrimAcc_MC%s",pname[i].Data()),
1469                                  Form("primary photon %s in acceptance: #phi ",ptype[i].Data()),
1470                                  nptbins,ptmin,ptmax,nphibins,phimin,phimax);
1471       fhPhiPrimMCAcc[i]->SetYTitle("#phi (rad)");
1472       fhPhiPrimMCAcc[i]->SetXTitle("E (GeV)");
1473       outputContainer->Add(fhPhiPrimMCAcc[i]) ;
1474
1475     }
1476
1477     if(fCheckConversion){
1478       fhPtConversionTagged  = new TH1F("hPtMCConversionTagged","Number of converted #gamma over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1479       fhPtConversionTagged->SetYTitle("N");
1480       fhPtConversionTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1481       outputContainer->Add(fhPtConversionTagged) ;
1482
1483
1484       fhPtAntiNeutronTagged  = new TH1F("hPtMCAntiNeutronTagged","Number of AntiNeutron id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1485       fhPtAntiNeutronTagged->SetYTitle("N");
1486       fhPtAntiNeutronTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1487       outputContainer->Add(fhPtAntiNeutronTagged) ;
1488
1489       fhPtAntiProtonTagged  = new TH1F("hPtMCAntiProtonTagged","Number of AntiProton id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1490       fhPtAntiProtonTagged->SetYTitle("N");
1491       fhPtAntiProtonTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1492       outputContainer->Add(fhPtAntiProtonTagged) ;
1493
1494       fhPtUnknownTagged  = new TH1F("hPtMCUnknownTagged","Number of Unknown id as Photon over calorimeter, tagged as converted",nptbins,ptmin,ptmax);
1495       fhPtUnknownTagged->SetYTitle("N");
1496       fhPtUnknownTagged->SetXTitle("p_{T #gamma}(GeV/c)");
1497       outputContainer->Add(fhPtUnknownTagged) ;
1498
1499       fhConvDeltaEtaMCConversion  = new TH2F
1500       ("hConvDeltaEtaMCConversion","#Delta #eta of selected conversion pairs from real conversions",100,0,fMassCut,netabins,-0.5,0.5);
1501       fhConvDeltaEtaMCConversion->SetYTitle("#Delta #eta");
1502       fhConvDeltaEtaMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1503       outputContainer->Add(fhConvDeltaEtaMCConversion) ;
1504
1505       fhConvDeltaPhiMCConversion  = new TH2F
1506       ("hConvDeltaPhiMCConversion","#Delta #phi of selected conversion pairs from real conversions",100,0,fMassCut,nphibins,-0.5,0.5);
1507       fhConvDeltaPhiMCConversion->SetYTitle("#Delta #phi");
1508       fhConvDeltaPhiMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1509       outputContainer->Add(fhConvDeltaPhiMCConversion) ;
1510
1511       fhConvDeltaEtaPhiMCConversion  = new TH2F
1512       ("hConvDeltaEtaPhiMCConversion","#Delta #eta vs #Delta #phi of selected conversion pairs, from real conversions",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1513       fhConvDeltaEtaPhiMCConversion->SetYTitle("#Delta #phi");
1514       fhConvDeltaEtaPhiMCConversion->SetXTitle("#Delta #eta");
1515       outputContainer->Add(fhConvDeltaEtaPhiMCConversion) ;
1516
1517       fhConvAsymMCConversion  = new TH2F
1518       ("hConvAsymMCConversion","Asymmetry of selected conversion pairs from real conversions",100,0,fMassCut,100,0,1);
1519       fhConvAsymMCConversion->SetYTitle("Asymmetry");
1520       fhConvAsymMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1521       outputContainer->Add(fhConvAsymMCConversion) ;
1522
1523       fhConvPtMCConversion  = new TH2F
1524       ("hConvPtMCConversion","p_{T} of selected conversion pairs from real conversions",100,0,fMassCut,100,0.,10.);
1525       fhConvPtMCConversion->SetYTitle("Pair p_{T} (GeV/c)");
1526       fhConvPtMCConversion->SetXTitle("Pair Mass (GeV/c^2)");
1527       outputContainer->Add(fhConvPtMCConversion) ;
1528
1529       fhConvDispersionMCConversion  = new TH2F
1530       ("hConvDispersionMCConversion","p_{T} of selected conversion pairs from real conversions",100,0.,1.,100,0.,1.);
1531       fhConvDispersionMCConversion->SetYTitle("Dispersion cluster 1");
1532       fhConvDispersionMCConversion->SetXTitle("Dispersion cluster 2");
1533       outputContainer->Add(fhConvDispersionMCConversion) ;
1534
1535       fhConvM02MCConversion  = new TH2F
1536       ("hConvM02MCConversion","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1537       fhConvM02MCConversion->SetYTitle("M02 cluster 1");
1538       fhConvM02MCConversion->SetXTitle("M02 cluster 2");
1539       outputContainer->Add(fhConvM02MCConversion) ;
1540
1541       fhConvDeltaEtaMCAntiNeutron  = new TH2F
1542       ("hConvDeltaEtaMCAntiNeutron","#Delta #eta of selected conversion pairs from anti-neutrons",100,0,fMassCut,netabins,-0.5,0.5);
1543       fhConvDeltaEtaMCAntiNeutron->SetYTitle("#Delta #eta");
1544       fhConvDeltaEtaMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1545       outputContainer->Add(fhConvDeltaEtaMCAntiNeutron) ;
1546
1547       fhConvDeltaPhiMCAntiNeutron  = new TH2F
1548       ("hConvDeltaPhiMCAntiNeutron","#Delta #phi of selected conversion pairs from anti-neutrons",100,0,fMassCut,nphibins,-0.5,0.5);
1549       fhConvDeltaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
1550       fhConvDeltaPhiMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1551       outputContainer->Add(fhConvDeltaPhiMCAntiNeutron) ;
1552
1553       fhConvDeltaEtaPhiMCAntiNeutron  = new TH2F
1554       ("hConvDeltaEtaPhiMCAntiNeutron","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-neutrons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1555       fhConvDeltaEtaPhiMCAntiNeutron->SetYTitle("#Delta #phi");
1556       fhConvDeltaEtaPhiMCAntiNeutron->SetXTitle("#Delta #eta");
1557       outputContainer->Add(fhConvDeltaEtaPhiMCAntiNeutron) ;
1558
1559       fhConvAsymMCAntiNeutron  = new TH2F
1560       ("hConvAsymMCAntiNeutron","Asymmetry of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0,1);
1561       fhConvAsymMCAntiNeutron->SetYTitle("Asymmetry");
1562       fhConvAsymMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1563       outputContainer->Add(fhConvAsymMCAntiNeutron) ;
1564
1565       fhConvPtMCAntiNeutron  = new TH2F
1566       ("hConvPtMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0,fMassCut,100,0.,10.);
1567       fhConvPtMCAntiNeutron->SetYTitle("Pair p_{T} (GeV/c)");
1568       fhConvPtMCAntiNeutron->SetXTitle("Pair Mass (GeV/c^2)");
1569       outputContainer->Add(fhConvPtMCAntiNeutron) ;
1570
1571       fhConvDispersionMCAntiNeutron  = new TH2F
1572       ("hConvDispersionMCAntiNeutron","p_{T} of selected conversion pairs from anti-neutrons",100,0.,1.,100,0.,1.);
1573       fhConvDispersionMCAntiNeutron->SetYTitle("Dispersion cluster 1");
1574       fhConvDispersionMCAntiNeutron->SetXTitle("Dispersion cluster 2");
1575       outputContainer->Add(fhConvDispersionMCAntiNeutron) ;
1576
1577       fhConvM02MCAntiNeutron  = new TH2F
1578       ("hConvM02MCAntiNeutron","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1579       fhConvM02MCAntiNeutron->SetYTitle("M02 cluster 1");
1580       fhConvM02MCAntiNeutron->SetXTitle("M02 cluster 2");
1581       outputContainer->Add(fhConvM02MCAntiNeutron) ;
1582
1583       fhConvDeltaEtaMCAntiProton  = new TH2F
1584       ("hConvDeltaEtaMCAntiProton","#Delta #eta of selected conversion pairs from anti-protons",100,0,fMassCut,netabins,-0.5,0.5);
1585       fhConvDeltaEtaMCAntiProton->SetYTitle("#Delta #eta");
1586       fhConvDeltaEtaMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1587       outputContainer->Add(fhConvDeltaEtaMCAntiProton) ;
1588
1589       fhConvDeltaPhiMCAntiProton  = new TH2F
1590       ("hConvDeltaPhiMCAntiProton","#Delta #phi of selected conversion pairs from anti-protons",100,0,fMassCut,nphibins,-0.5,0.5);
1591       fhConvDeltaPhiMCAntiProton->SetYTitle("#Delta #phi");
1592       fhConvDeltaPhiMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1593       outputContainer->Add(fhConvDeltaPhiMCAntiProton) ;
1594
1595       fhConvDeltaEtaPhiMCAntiProton  = new TH2F
1596       ("hConvDeltaEtaPhiMCAntiProton","#Delta #eta vs #Delta #phi of selected conversion pairs from anti-protons",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1597       fhConvDeltaEtaPhiMCAntiProton->SetYTitle("#Delta #phi");
1598       fhConvDeltaEtaPhiMCAntiProton->SetXTitle("#Delta #eta");
1599       outputContainer->Add(fhConvDeltaEtaPhiMCAntiProton) ;
1600
1601       fhConvAsymMCAntiProton  = new TH2F
1602       ("hConvAsymMCAntiProton","Asymmetry of selected conversion pairs from anti-protons",100,0,fMassCut,100,0,1);
1603       fhConvAsymMCAntiProton->SetYTitle("Asymmetry");
1604       fhConvAsymMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1605       outputContainer->Add(fhConvAsymMCAntiProton) ;
1606
1607       fhConvPtMCAntiProton  = new TH2F
1608       ("hConvPtMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0,fMassCut,100,0.,10.);
1609       fhConvPtMCAntiProton->SetYTitle("Pair p_{T} (GeV/c)");
1610       fhConvPtMCAntiProton->SetXTitle("Pair Mass (GeV/c^2)");
1611       outputContainer->Add(fhConvPtMCAntiProton) ;
1612
1613       fhConvDispersionMCAntiProton  = new TH2F
1614       ("hConvDispersionMCAntiProton","p_{T} of selected conversion pairs from anti-protons",100,0.,1.,100,0.,1.);
1615       fhConvDispersionMCAntiProton->SetYTitle("Dispersion cluster 1");
1616       fhConvDispersionMCAntiProton->SetXTitle("Dispersion cluster 2");
1617       outputContainer->Add(fhConvDispersionMCAntiProton) ;
1618
1619       fhConvM02MCAntiProton  = new TH2F
1620       ("hConvM02MCAntiProton","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1621       fhConvM02MCAntiProton->SetYTitle("M02 cluster 1");
1622       fhConvM02MCAntiProton->SetXTitle("M02 cluster 2");
1623       outputContainer->Add(fhConvM02MCAntiProton) ;
1624
1625       fhConvDeltaEtaMCString  = new TH2F
1626       ("hConvDeltaEtaMCString","#Delta #eta of selected conversion pairs from string",100,0,fMassCut,netabins,-0.5,0.5);
1627       fhConvDeltaEtaMCString->SetYTitle("#Delta #eta");
1628       fhConvDeltaEtaMCString->SetXTitle("Pair Mass (GeV/c^2)");
1629       outputContainer->Add(fhConvDeltaEtaMCString) ;
1630
1631       fhConvDeltaPhiMCString  = new TH2F
1632       ("hConvDeltaPhiMCString","#Delta #phi of selected conversion pairs from string",100,0,fMassCut,nphibins,-0.5,0.5);
1633       fhConvDeltaPhiMCString->SetYTitle("#Delta #phi");
1634       fhConvDeltaPhiMCString->SetXTitle("Pair Mass (GeV/c^2)");
1635       outputContainer->Add(fhConvDeltaPhiMCString) ;
1636
1637       fhConvDeltaEtaPhiMCString  = new TH2F
1638       ("hConvDeltaEtaPhiMCString","#Delta #eta vs #Delta #phi of selected conversion pairs from string",netabins,-0.5,0.5,nphibins,-0.5,0.5);
1639       fhConvDeltaEtaPhiMCString->SetYTitle("#Delta #phi");
1640       fhConvDeltaEtaPhiMCString->SetXTitle("#Delta #eta");
1641       outputContainer->Add(fhConvDeltaEtaPhiMCString) ;
1642
1643       fhConvAsymMCString  = new TH2F
1644       ("hConvAsymMCString","Asymmetry of selected conversion pairs from string",100,0,fMassCut,100,0,1);
1645       fhConvAsymMCString->SetYTitle("Asymmetry");
1646       fhConvAsymMCString->SetXTitle("Pair Mass (GeV/c^2)");
1647       outputContainer->Add(fhConvAsymMCString) ;
1648
1649       fhConvPtMCString  = new TH2F
1650       ("hConvPtMCString","p_{T} of selected conversion pairs from string",100,0,fMassCut,100,0.,10.);
1651       fhConvPtMCString->SetYTitle("Pair p_{T} (GeV/c)");
1652       fhConvPtMCString->SetXTitle("Pair Mass (GeV/c^2)");
1653       outputContainer->Add(fhConvPtMCString) ;
1654
1655       fhConvDispersionMCString  = new TH2F
1656       ("hConvDispersionMCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1657       fhConvDispersionMCString->SetYTitle("Dispersion cluster 1");
1658       fhConvDispersionMCString->SetXTitle("Dispersion cluster 2");
1659       outputContainer->Add(fhConvDispersionMCString) ;
1660
1661       fhConvM02MCString  = new TH2F
1662       ("hConvM02MCString","p_{T} of selected conversion pairs from string",100,0.,1.,100,0.,1.);
1663       fhConvM02MCString->SetYTitle("M02 cluster 1");
1664       fhConvM02MCString->SetXTitle("M02 cluster 2");
1665       outputContainer->Add(fhConvM02MCString) ;
1666
1667       fhConvDistMCConversion  = new TH2F
1668       ("hConvDistMCConversion","calculated conversion distance vs real vertes for MC conversion",100,0.,5.,100,0.,5.);
1669       fhConvDistMCConversion->SetYTitle("distance");
1670       fhConvDistMCConversion->SetXTitle("vertex R");
1671       outputContainer->Add(fhConvDistMCConversion) ;
1672
1673       fhConvDistMCConversionCuts  = new TH2F
1674       ("hConvDistMCConversionCuts","calculated conversion distance vs real vertes for MC conversion, deta < 0.05, m < 10 MeV, asym < 0.1",100,0.,5.,100,0.,5.);
1675       fhConvDistMCConversionCuts->SetYTitle("distance");
1676       fhConvDistMCConversionCuts->SetXTitle("vertex R");
1677       outputContainer->Add(fhConvDistMCConversionCuts) ;
1678     }
1679
1680     if(fFillSSHistograms){
1681
1682       TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
1683
1684       TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
1685
1686       for(Int_t i = 0; i < 6; i++){
1687
1688         fhMCELambda0[i]  = new TH2F(Form("hELambda0_MC%s",pnamess[i].Data()),
1689                                     Form("cluster from %s : E vs #lambda_{0}^{2}",ptypess[i].Data()),
1690                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1691         fhMCELambda0[i]->SetYTitle("#lambda_{0}^{2}");
1692         fhMCELambda0[i]->SetXTitle("E (GeV)");
1693         outputContainer->Add(fhMCELambda0[i]) ;
1694
1695         fhMCEdLambda0[i]  = new TH2F(Form("hEdLambda0_MC%s",pnamess[i].Data()),
1696                                      Form("cluster from %s : E vs #lambda_{0}^{2}/N_{cells}^{2}",ptypess[i].Data()),
1697                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1698         fhMCEdLambda0[i]->SetYTitle("d#lambda_{0}^{2}");
1699         fhMCEdLambda0[i]->SetXTitle("E (GeV)");
1700         outputContainer->Add(fhMCEdLambda0[i]) ;
1701
1702         fhMCELambda1[i]  = new TH2F(Form("hELambda1_MC%s",pnamess[i].Data()),
1703                                     Form("cluster from %s : E vs #lambda_{1}^{2}",ptypess[i].Data()),
1704                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1705         fhMCELambda1[i]->SetYTitle("#lambda_{1}^{2}");
1706         fhMCELambda1[i]->SetXTitle("E (GeV)");
1707         outputContainer->Add(fhMCELambda1[i]) ;
1708
1709         fhMCEdLambda1[i]  = new TH2F(Form("hEdLambda1_MC%s",pnamess[i].Data()),
1710                                      Form("cluster from %s : E vs d#lambda_{1}^{2}/N_{cells}^{2}",ptypess[i].Data()),
1711                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1712         fhMCEdLambda1[i]->SetYTitle("d#lambda_{1}^{2}");
1713         fhMCEdLambda1[i]->SetXTitle("E (GeV)");
1714         outputContainer->Add(fhMCEdLambda1[i]) ;
1715
1716         fhMCEDispersion[i]  = new TH2F(Form("hEDispersion_MC%s",pnamess[i].Data()),
1717                                        Form("cluster from %s : E vs dispersion^{2}",ptypess[i].Data()),
1718                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1719         fhMCEDispersion[i]->SetYTitle("D^{2}");
1720         fhMCEDispersion[i]->SetXTitle("E (GeV)");
1721         outputContainer->Add(fhMCEDispersion[i]) ;
1722
1723         fhMCEdDispersion[i]  = new TH2F(Form("hEdDispersion_MC%s",pnamess[i].Data()),
1724                                         Form("cluster from %s : E vs dispersion^{2}/N_{cells}^{2}",ptypess[i].Data()),
1725                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1726         fhMCEdDispersion[i]->SetYTitle("dD^{2}");
1727         fhMCEdDispersion[i]->SetXTitle("E (GeV)");
1728         outputContainer->Add(fhMCEdDispersion[i]) ;
1729
1730        }// loop
1731
1732       if(!GetReader()->IsEmbeddedClusterSelectionOn())
1733       {
1734         fhMCPhotonELambda0NoOverlap  = new TH2F("hELambda0_MCPhoton_NoOverlap",
1735                                                 "cluster from Photon : E vs #lambda_{0}^{2}",
1736                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1737         fhMCPhotonELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
1738         fhMCPhotonELambda0NoOverlap->SetXTitle("E (GeV)");
1739         outputContainer->Add(fhMCPhotonELambda0NoOverlap) ;
1740
1741         fhMCPhotonEdLambda0NoOverlap  = new TH2F("hEdLambda0_MCPhoton_NoOverlap",
1742                                                  "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1743                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1744         fhMCPhotonEdLambda0NoOverlap->SetYTitle("d#lambda_{0}^{2}");
1745         fhMCPhotonEdLambda0NoOverlap->SetXTitle("E (GeV)");
1746         outputContainer->Add(fhMCPhotonEdLambda0NoOverlap) ;
1747
1748         fhMCPhotonELambda0TwoOverlap  = new TH2F("hELambda0_MCPhoton_TwoOverlap",
1749                                                  "cluster from Photon : E vs #lambda_{0}^{2}",
1750                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1751         fhMCPhotonELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
1752         fhMCPhotonELambda0TwoOverlap->SetXTitle("E (GeV)");
1753         outputContainer->Add(fhMCPhotonELambda0TwoOverlap) ;
1754
1755         fhMCPhotonEdLambda0TwoOverlap  = new TH2F("hEdLambda0_MCPhoton_TwoOverlap",
1756                                                   "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1757                                                   nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1758         fhMCPhotonEdLambda0TwoOverlap->SetYTitle("d#lambda_{0}^{2}");
1759         fhMCPhotonEdLambda0TwoOverlap->SetXTitle("E (GeV)");
1760         outputContainer->Add(fhMCPhotonEdLambda0TwoOverlap) ;
1761
1762         fhMCPhotonELambda0NOverlap  = new TH2F("hELambda0_MCPhoton_NOverlap",
1763                                                "cluster from Photon : E vs #lambda_{0}^{2}",
1764                                                nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1765         fhMCPhotonELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1766         fhMCPhotonELambda0NOverlap->SetXTitle("E (GeV)");
1767         outputContainer->Add(fhMCPhotonELambda0NOverlap) ;
1768
1769         fhMCPhotonEdLambda0NOverlap  = new TH2F("hEdLambda0_MCPhoton_NOverlap",
1770                                                 "cluster from Photon : E vs #lambda_{0}^{2}/N_{cells}^{2}",
1771                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax/50);
1772         fhMCPhotonEdLambda0NOverlap->SetYTitle("d#lambda_{0}^{2}");
1773         fhMCPhotonEdLambda0NOverlap->SetXTitle("E (GeV)");
1774         outputContainer->Add(fhMCPhotonEdLambda0NOverlap) ;
1775       } //No embedding
1776
1777       //Fill histograms to check shape of embedded clusters
1778       if(GetReader()->IsEmbeddedClusterSelectionOn())
1779       {
1780
1781         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
1782                                                     "Energy Fraction of embedded signal versus cluster energy",
1783                                                     nptbins,ptmin,ptmax,100,0.,1.);
1784         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1785         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1786         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1787
1788         fhEmbedPhotonELambda0FullSignal  = new TH2F("hELambda0_EmbedPhoton_FullSignal",
1789                                                 "cluster from Photon embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1790                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1791         fhEmbedPhotonELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1792         fhEmbedPhotonELambda0FullSignal->SetXTitle("E (GeV)");
1793         outputContainer->Add(fhEmbedPhotonELambda0FullSignal) ;
1794
1795         fhEmbedPhotonEdLambda0FullSignal  = new TH2F("hEdLambda0_EmbedPhoton_FullSignal",
1796                                           "cluster from Photon with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
1797                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1798         fhEmbedPhotonEdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1799         fhEmbedPhotonEdLambda0FullSignal->SetXTitle("E (GeV)");
1800         outputContainer->Add(fhEmbedPhotonEdLambda0FullSignal) ;
1801
1802
1803         fhEmbedPhotonELambda0MostlySignal  = new TH2F("hELambda0_EmbedPhoton_MostlySignal",
1804                                           "cluster from Photon embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1805                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1806         fhEmbedPhotonELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1807         fhEmbedPhotonELambda0MostlySignal->SetXTitle("E (GeV)");
1808         outputContainer->Add(fhEmbedPhotonELambda0MostlySignal) ;
1809
1810         fhEmbedPhotonEdLambda0MostlySignal  = new TH2F("hEdLambda0_EmbedPhoton_MostlySignal",
1811                                            "cluster from Photon with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
1812                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1813         fhEmbedPhotonEdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1814         fhEmbedPhotonEdLambda0MostlySignal->SetXTitle("E (GeV)");
1815         outputContainer->Add(fhEmbedPhotonEdLambda0MostlySignal) ;
1816
1817
1818         fhEmbedPhotonELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPhoton_MostlyBkg",
1819                                           "cluster from Photon embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1820                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1821         fhEmbedPhotonELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1822         fhEmbedPhotonELambda0MostlyBkg->SetXTitle("E (GeV)");
1823         outputContainer->Add(fhEmbedPhotonELambda0MostlyBkg) ;
1824
1825         fhEmbedPhotonEdLambda0MostlyBkg  = new TH2F("hEdLambda0_EmbedPhoton_MostlyBkg",
1826                                            "cluster from Photon with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
1827                                            nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1828         fhEmbedPhotonEdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1829         fhEmbedPhotonEdLambda0MostlyBkg->SetXTitle("E (GeV)");
1830         outputContainer->Add(fhEmbedPhotonEdLambda0MostlyBkg) ;
1831
1832
1833         fhEmbedPhotonELambda0FullBkg  = new TH2F("hELambda0_EmbedPhoton_FullBkg",
1834                                                    "cluster from Photonm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1835                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1836         fhEmbedPhotonELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1837         fhEmbedPhotonELambda0FullBkg->SetXTitle("E (GeV)");
1838         outputContainer->Add(fhEmbedPhotonELambda0FullBkg) ;
1839
1840         fhEmbedPhotonEdLambda0FullBkg  = new TH2F("hEdLambda0_EmbedPhoton_FullBkg",
1841                                                     "cluster from Photon with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
1842                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1843         fhEmbedPhotonEdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1844         fhEmbedPhotonEdLambda0FullBkg->SetXTitle("E (GeV)");
1845         outputContainer->Add(fhEmbedPhotonEdLambda0FullBkg) ;
1846
1847
1848         fhEmbedPi0ELambda0FullSignal  = new TH2F("hELambda0_EmbedPi0_FullSignal",
1849                                                     "cluster from Pi0 embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1850                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1851         fhEmbedPi0ELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1852         fhEmbedPi0ELambda0FullSignal->SetXTitle("E (GeV)");
1853         outputContainer->Add(fhEmbedPi0ELambda0FullSignal) ;
1854
1855         fhEmbedPi0EdLambda0FullSignal  = new TH2F("hEdLambda0_EmbedPi0_FullSignal",
1856                                                      "cluster from Pi0 with more than 90% energy in cluster: E vs d#lambda_{0}^{2}",
1857                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1858         fhEmbedPi0EdLambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1859         fhEmbedPi0EdLambda0FullSignal->SetXTitle("E (GeV)");
1860         outputContainer->Add(fhEmbedPi0EdLambda0FullSignal) ;
1861
1862
1863         fhEmbedPi0ELambda0MostlySignal  = new TH2F("hELambda0_EmbedPi0_MostlySignal",
1864                                                       "cluster from Pi0 embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1865                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1866         fhEmbedPi0ELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1867         fhEmbedPi0ELambda0MostlySignal->SetXTitle("E (GeV)");
1868         outputContainer->Add(fhEmbedPi0ELambda0MostlySignal) ;
1869
1870         fhEmbedPi0EdLambda0MostlySignal  = new TH2F("hEdLambda0_EmbedPi0_MostlySignal",
1871                                                        "cluster from Pi0 with 50% to 90% energy in cluster: E vs d#lambda_{0}^{2}",
1872                                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1873         fhEmbedPi0EdLambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1874         fhEmbedPi0EdLambda0MostlySignal->SetXTitle("E (GeV)");
1875         outputContainer->Add(fhEmbedPi0EdLambda0MostlySignal) ;
1876
1877
1878         fhEmbedPi0ELambda0MostlyBkg  = new TH2F("hELambda0_EmbedPi0_MostlyBkg",
1879                                                    "cluster from Pi0 embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1880                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1881         fhEmbedPi0ELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1882         fhEmbedPi0ELambda0MostlyBkg->SetXTitle("E (GeV)");
1883         outputContainer->Add(fhEmbedPi0ELambda0MostlyBkg) ;
1884
1885         fhEmbedPi0EdLambda0MostlyBkg  = new TH2F("hEdLambda0_EmbedPi0_MostlyBkg",
1886                                                     "cluster from Pi0 with 10% to 50% energy in cluster: E vs d#lambda_{0}^{2}",
1887                                                     nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1888         fhEmbedPi0EdLambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1889         fhEmbedPi0EdLambda0MostlyBkg->SetXTitle("E (GeV)");
1890         outputContainer->Add(fhEmbedPi0EdLambda0MostlyBkg) ;
1891
1892
1893         fhEmbedPi0ELambda0FullBkg  = new TH2F("hELambda0_EmbedPi0_FullBkg",
1894                                                  "cluster from Pi0 embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1895                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1896         fhEmbedPi0ELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1897         fhEmbedPi0ELambda0FullBkg->SetXTitle("E (GeV)");
1898         outputContainer->Add(fhEmbedPi0ELambda0FullBkg) ;
1899
1900         fhEmbedPi0EdLambda0FullBkg  = new TH2F("hEdLambda0_EmbedPi0_FullBkg",
1901                                                   "cluster from Pi0 with 0% to 10% energy in cluster: E vs d#lambda_{0}^{2}",
1902                                                   nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1903         fhEmbedPi0EdLambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1904         fhEmbedPi0EdLambda0FullBkg->SetXTitle("E (GeV)");
1905         outputContainer->Add(fhEmbedPi0EdLambda0FullBkg) ;
1906
1907       }// embedded histograms
1908
1909
1910     }// Fill SS MC histograms
1911
1912   }//Histos with MC
1913
1914   //Store calo PID histograms
1915   if(fRejectTrackMatch){
1916     TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
1917     for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
1918       outputContainer->Add(caloPIDHistos->At(i)) ;
1919     }
1920     delete caloPIDHistos;
1921   }
1922
1923   return outputContainer ;
1924
1925 }
1926
1927 //____________________________________________________________________________
1928 void AliAnaPhoton::Init()
1929 {
1930
1931   //Init
1932   //Do some checks
1933   if(fCalorimeter == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1934     printf("AliAnaPhoton::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1935     abort();
1936   }
1937   else  if(fCalorimeter == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1938     printf("AliAnaPhoton::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1939     abort();
1940   }
1941
1942 }
1943
1944 //____________________________________________________________________________
1945 void AliAnaPhoton::InitParameters()
1946 {
1947
1948   //Initialize the parameters of the analysis.
1949   AddToHistogramsName("AnaPhoton_");
1950
1951   fCalorimeter = "EMCAL" ;
1952   fMinDist     = 2.;
1953   fMinDist2    = 4.;
1954   fMinDist3    = 5.;
1955   fMassCut     = 0.03; //30 MeV
1956
1957   fTimeCutMin  = -1;
1958   fTimeCutMax  = 9999999;
1959   fNCellsCut   = 0;
1960
1961   fRejectTrackMatch       = kTRUE ;
1962   fCheckConversion        = kFALSE;
1963   fRemoveConvertedPair    = kFALSE;
1964   fAddConvertedPairsToAOD = kFALSE;
1965
1966 }
1967
1968 //__________________________________________________________________
1969 void  AliAnaPhoton::MakeAnalysisFillAOD()
1970 {
1971   //Do photon analysis and fill aods
1972
1973   //Get the vertex
1974   Double_t v[3] = {0,0,0}; //vertex ;
1975   GetReader()->GetVertex(v);
1976
1977   //Select the Calorimeter of the photon
1978   TObjArray * pl = 0x0;
1979   if(fCalorimeter == "PHOS")
1980     pl = GetPHOSClusters();
1981   else if (fCalorimeter == "EMCAL")
1982     pl = GetEMCALClusters();
1983
1984   if(!pl) {
1985     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",fCalorimeter.Data());
1986     return;
1987   }
1988
1989   //Init arrays, variables, get number of clusters
1990   TLorentzVector mom, mom2 ;
1991   Int_t nCaloClusters = pl->GetEntriesFast();
1992   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1993   Bool_t * indexConverted = 0x0;
1994   if(fCheckConversion){
1995     indexConverted = new Bool_t[nCaloClusters];
1996     for (Int_t i = 0; i < nCaloClusters; i++)
1997       indexConverted[i] = kFALSE;
1998         }
1999
2000   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillAOD() - input %s cluster entries %d\n", fCalorimeter.Data(), nCaloClusters);
2001
2002   //----------------------------------------------------
2003   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
2004   //----------------------------------------------------
2005   // Loop on clusters
2006   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++){
2007
2008           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
2009     //printf("calo %d, %f\n",icalo,calo->E());
2010
2011     //Get the index where the cluster comes, to retrieve the corresponding vertex
2012     Int_t evtIndex = 0 ;
2013     if (GetMixedEvent()) {
2014       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
2015       //Get the vertex and check it is not too large in z
2016       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
2017     }
2018
2019     //Cluster selection, not charged, with photon id and in fiducial cut
2020     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
2021       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
2022     else{
2023       Double_t vertex[]={0,0,0};
2024       calo->GetMomentum(mom,vertex) ;
2025     }
2026
2027     //--------------------------------------
2028     // Cluster selection
2029     //--------------------------------------
2030     if(!ClusterSelected(calo,mom)) continue;
2031
2032     //----------------------------
2033     //Create AOD for analysis
2034     //----------------------------
2035     AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
2036
2037     //...............................................
2038     //Set the indeces of the original caloclusters (MC, ID), and calorimeter
2039     Int_t label = calo->GetLabel();
2040     aodph.SetLabel(label);
2041     aodph.SetCaloLabel(calo->GetID(),-1);
2042     aodph.SetDetector(fCalorimeter);
2043     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
2044
2045     //...............................................
2046     //Set bad channel distance bit
2047     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
2048     if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
2049     else if(distBad > fMinDist2) aodph.SetDistToBad(1) ;
2050     else                         aodph.SetDistToBad(0) ;
2051     //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
2052
2053     //--------------------------------------------------------------------------------------
2054     //Play with the MC stack if available
2055     //--------------------------------------------------------------------------------------
2056
2057     //Check origin of the candidates
2058     if(IsDataMC()){
2059       aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
2060
2061       if(GetDebug() > 0)
2062         printf("AliAnaPhoton::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
2063     }//Work with stack also
2064
2065     //--------------------------------------------------------------------------------------
2066     //Fill some shower shape histograms before PID is applied
2067     //--------------------------------------------------------------------------------------
2068
2069     FillShowerShapeHistograms(calo,aodph.GetTag());
2070
2071     //-------------------------------------
2072     //PID selection or bit setting
2073     //-------------------------------------
2074     // MC
2075     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
2076       //Get most probable PID, check PID weights (in MC this option is mandatory)
2077       aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
2078       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
2079       //If primary is not photon, skip it.
2080       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
2081     }
2082     //...............................................
2083     // Data, PID check on
2084     else if(IsCaloPIDOn()){
2085       //Get most probable PID, 2 options check PID weights
2086       //or redo PID, recommended option for MCEal.
2087       if(!IsCaloPIDRecalculationOn())
2088         aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
2089       else
2090         aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
2091
2092       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
2093
2094       //If cluster does not pass pid, not photon, skip it.
2095       if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
2096
2097     }
2098     //...............................................
2099     // Data, PID check off
2100     else{
2101       //Set PID bits for later selection (AliAnaPi0 for example)
2102       //GetPDG already called in SetPIDBits.
2103       GetCaloPID()->SetPIDBits(fCalorimeter,calo,&aodph, GetCaloUtils());
2104       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");
2105     }
2106
2107     if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
2108
2109
2110     //--------------------------------------------------------------------------------------
2111     // Conversions pairs analysis
2112     // Check if cluster comes from a conversion in the material in front of the calorimeter
2113     // Do invariant mass of all pairs, if mass is close to 0, then it is conversion.
2114     //--------------------------------------------------------------------------------------
2115
2116     // Do analysis only if there are more than one cluster
2117     if( nCaloClusters > 1 && fCheckConversion){
2118       Bool_t bConverted = kFALSE;
2119       Int_t id2 = -1;
2120
2121       //Check if set previously as converted couple, if so skip its use.
2122       if (indexConverted[icalo]) continue;
2123
2124       // Second cluster loop
2125       for(Int_t jcalo = icalo + 1 ; jcalo < nCaloClusters ; jcalo++) {
2126         //Check if set previously as converted couple, if so skip its use.
2127         if (indexConverted[jcalo]) continue;
2128         //printf("Check Conversion indeces %d and %d\n",icalo,jcalo);
2129         AliVCluster * calo2 =  (AliVCluster*) (pl->At(jcalo));  //Get cluster kinematics
2130
2131
2132         //Mixed event, get index of event
2133         Int_t evtIndex2 = 0 ;
2134         if (GetMixedEvent()) {
2135           evtIndex2=GetMixedEvent()->EventIndexForCaloCluster(calo2->GetID()) ;
2136
2137         }
2138
2139         //Get kinematics of second cluster
2140         if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
2141           calo2->GetMomentum(mom2,GetVertex(evtIndex2)) ;}//Assume that come from vertex in straight line
2142         else{
2143           Double_t vertex[]={0,0,0};
2144           calo2->GetMomentum(mom2,vertex) ;
2145         }
2146
2147         //--------------------------------------
2148         // Cluster selection
2149         //--------------------------------------
2150
2151         if(!ClusterSelected(calo2,mom2)) continue;
2152
2153         //................................................
2154         // Get TOF of each cluster in pair, calculate difference if small,
2155         // take this pair. Way to reject clusters from hadrons (or pileup?)
2156         Double_t t12diff = calo2->GetTOF()-calo->GetTOF()*1e9;
2157         if(TMath::Abs(t12diff) > GetPairTimeCut()) continue;
2158
2159         //................................................
2160         //Get mass of pair, if small, take this pair.
2161         Float_t pairM     = (mom+mom2).M();
2162         //printf("\t both in calo, mass %f, cut %f\n",pairM,fMassCut);
2163         if(pairM < fMassCut){
2164           aodph.SetTagged(kFALSE);
2165           id2 = calo2->GetID();
2166           indexConverted[icalo]=kTRUE;
2167           indexConverted[jcalo]=kTRUE;
2168           Float_t asymmetry = TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E());
2169           Float_t dPhi      = mom.Phi()-mom2.Phi();
2170           Float_t dEta      = mom.Eta()-mom2.Eta();
2171
2172           //...............................................
2173           //Fill few histograms with kinematics of the pair
2174           //FIXME, move all this to MakeAnalysisFillHistograms ...
2175
2176           fhConvDeltaEta   ->Fill( pairM, dPhi      );
2177           fhConvDeltaPhi   ->Fill( pairM, dEta      );
2178           fhConvAsym       ->Fill( pairM, asymmetry );
2179           fhConvDeltaEtaPhi->Fill( dEta , dPhi      );
2180           fhConvPt         ->Fill( pairM, (mom+mom2).Pt());
2181
2182           //Estimate conversion distance, T. Awes, M. Ivanov
2183           //Under the assumption that the pair has zero mass, and that each electron
2184           //of the pair has the same momentum, they will each have the same bend radius
2185           //given by R=p/(qB) = p / (300 B) with p in [MeV/c], B in [Tesla] and R in [m].
2186           //With nominal ALICE magnet current of 30kA B=0.5T, and so with E_cluster=p,
2187           //R = E/1.5 [cm].  Under these assumptions, the distance from the conversion
2188           //point to the MCEal can be related to the separation distance, L=2y, on the MCEal
2189           //as d = sqrt(R^2 -(R-y)^2) = sqrt(2Ry - y^2). And since R>>y we can write as
2190           //d = sqrt(E*L/1.5) where E is the cluster energy and L is the distance in cm between
2191           //the clusters.
2192           Float_t pos1[3];
2193           calo->GetPosition(pos1);
2194           Float_t pos2[3];
2195           calo2->GetPosition(pos2);
2196           Float_t clustDist = TMath::Sqrt((pos1[0]-pos2[0])*(pos1[0]-pos2[0])+
2197                                           (pos1[1]-pos2[1])*(pos1[1]-pos2[1])+
2198                                           (pos1[2]-pos2[2])*(pos1[2]-pos2[2]));
2199
2200           Float_t convDist  = TMath::Sqrt(mom.E() *clustDist*0.01/0.15);
2201           Float_t convDist2 = TMath::Sqrt(mom2.E()*clustDist*0.01/0.15);
2202           //printf("l = %f, e1 = %f, d1=%f, e2 = %f, d2=%f\n",clustDist,mom.E(),convDist,mom2.E(),convDist2);
2203           if(GetDebug() > 2)
2204             printf("AliAnaPhoton::MakeAnalysisFillAOD(): Pair with mass %2.3f < %2.3f, %1.2f < dPhi %2.2f < %2.2f, dEta %f < %2.2f, asymmetry %2.2f< %2.2f; \n    cluster1 id %d, e %2.3f  SM %d, eta %2.3f, phi %2.3f ; \n    cluster2 id %d, e %2.3f, SM %d,eta %2.3f, phi %2.3f\n",
2205                    pairM,fMassCut,fConvDPhiMinCut, dPhi, fConvDPhiMaxCut, dEta, fConvDEtaCut, asymmetry, fConvAsymCut,
2206                    calo->GetID(),calo->E(),GetCaloUtils()->GetModuleNumber(calo), mom.Eta(), mom.Phi(),
2207                    id2, calo2->E(), GetCaloUtils()->GetModuleNumber(calo2),mom2.Eta(), mom2.Phi());
2208
2209           fhConvDistEta ->Fill(mom .Eta(), convDist );
2210           fhConvDistEta ->Fill(mom2.Eta(), convDist2);
2211           fhConvDistEn  ->Fill(mom .E(),   convDist );
2212           fhConvDistEn  ->Fill(mom2.E(),   convDist2);
2213           fhConvDistMass->Fill((mom+mom2).M(), convDist );
2214           //dEta cut
2215           if(dEta<0.05){
2216             fhConvDistEtaCutEta ->Fill(mom .Eta(), convDist );
2217             fhConvDistEtaCutEta ->Fill(mom2.Eta(), convDist2);
2218             fhConvDistEnCutEta  ->Fill(mom .E(),   convDist );
2219             fhConvDistEnCutEta  ->Fill(mom2.E(),   convDist2);
2220             fhConvDistMassCutEta->Fill((mom+mom2).M(), convDist );
2221             //mass cut
2222             if(pairM<0.01){//10 MeV
2223               fhConvDistEtaCutMass ->Fill(mom .Eta(), convDist );
2224               fhConvDistEtaCutMass ->Fill(mom2.Eta(), convDist2);
2225               fhConvDistEnCutMass  ->Fill(mom .E(),   convDist );
2226               fhConvDistEnCutMass  ->Fill(mom2.E(),   convDist2);
2227               // asymmetry cut
2228               if(asymmetry<0.1){
2229                 fhConvDistEtaCutAsy ->Fill(mom .Eta(), convDist );
2230                 fhConvDistEtaCutAsy ->Fill(mom2.Eta(), convDist2);
2231                 fhConvDistEnCutAsy  ->Fill(mom .E(),   convDist );
2232                 fhConvDistEnCutAsy  ->Fill(mom2.E(),   convDist2);
2233               }//asymmetry cut
2234             }//mass cut
2235           }//dEta cut
2236
2237           //...............................................
2238           //Select pairs in a eta-phi window
2239           if(TMath::Abs(dEta) < fConvDEtaCut    &&
2240              TMath::Abs(dPhi) < fConvDPhiMaxCut &&
2241              TMath::Abs(dPhi) > fConvDPhiMinCut &&
2242              asymmetry        < fConvAsymCut       ){
2243             bConverted = kTRUE;
2244           }
2245           //printf("Accepted? %d\n",bConverted);
2246           //...........................................
2247           //Fill more histograms, simulated data
2248           //FIXME, move all this to MakeAnalysisFillHistograms ...
2249           if(IsDataMC()){
2250
2251             //Check the origin of the pair, look for conversion, antinucleons or jet correlations (strings)
2252             Int_t ancPDG    = 0;
2253             Int_t ancStatus = 0;
2254             TLorentzVector momentum;
2255             TVector3 prodVertex;
2256             Int_t ancLabel  = GetMCAnalysisUtils()->CheckCommonAncestor(calo->GetLabel(), calo2->GetLabel(),
2257                                                                         GetReader(), ancPDG, ancStatus, momentum, prodVertex);
2258
2259             // printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Common ancestor label %d, pdg %d, name %s, status %d; \n",
2260             //                          ancLabel,ancPDG,TDatabasePDG::Instance()->GetParticle(ancPDG)->GetName(),ancStatus);
2261
2262             Int_t tag2 = GetMCAnalysisUtils()->CheckOrigin(calo2->GetLabels(),calo2->GetNLabels(),GetReader(), 0);
2263             if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCConversion)){
2264               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCConversion) && (ancPDG==22 || TMath::Abs(ancPDG)==11) && ancLabel > -1){
2265                 fhConvDeltaEtaMCConversion   ->Fill( pairM, dEta      );
2266                 fhConvDeltaPhiMCConversion   ->Fill( pairM, dPhi      );
2267                 fhConvAsymMCConversion       ->Fill( pairM, asymmetry );
2268                 fhConvDeltaEtaPhiMCConversion->Fill( dEta , dPhi      );
2269                 fhConvPtMCConversion         ->Fill( pairM, (mom+mom2).Pt());
2270                 fhConvDispersionMCConversion ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2271                 fhConvM02MCConversion        ->Fill( calo->GetM02(), calo2->GetM02());
2272                 fhConvDistMCConversion       ->Fill( convDist , prodVertex.Mag() );
2273                 fhConvDistMCConversion       ->Fill( convDist2, prodVertex.Mag() );
2274
2275                 if(dEta<0.05 && pairM<0.01 && asymmetry<0.1){
2276                   fhConvDistMCConversionCuts->Fill( convDist , prodVertex.Mag() );
2277                   fhConvDistMCConversionCuts->Fill( convDist2, prodVertex.Mag() );
2278                 }
2279
2280               }
2281             }
2282             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiNeutron)){
2283               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiNeutron) && ancPDG==-2112 && ancLabel > -1){
2284                 fhConvDeltaEtaMCAntiNeutron    ->Fill( pairM, dEta      );
2285                 fhConvDeltaPhiMCAntiNeutron    ->Fill( pairM, dPhi      );
2286                 fhConvAsymMCAntiNeutron        ->Fill( pairM, asymmetry );
2287                 fhConvDeltaEtaPhiMCAntiNeutron ->Fill( dEta , dPhi      );
2288                 fhConvPtMCAntiNeutron          ->Fill( pairM, (mom+mom2).Pt());
2289                 fhConvDispersionMCAntiNeutron  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2290                 fhConvM02MCAntiNeutron         ->Fill( calo->GetM02(), calo2->GetM02());
2291               }
2292             }
2293             else if(GetMCAnalysisUtils()->CheckTagBit(aodph.GetTag(),AliMCAnalysisUtils::kMCAntiProton)){
2294               if(GetMCAnalysisUtils()->CheckTagBit(tag2,AliMCAnalysisUtils::kMCAntiProton) && ancPDG==-2212 && ancLabel > -1){
2295                 fhConvDeltaEtaMCAntiProton    ->Fill( pairM, dEta      );
2296                 fhConvDeltaPhiMCAntiProton    ->Fill( pairM, dPhi      );
2297                 fhConvAsymMCAntiProton        ->Fill( pairM, asymmetry );
2298                 fhConvDeltaEtaPhiMCAntiProton ->Fill( dEta , dPhi      );
2299                 fhConvPtMCAntiProton          ->Fill( pairM, (mom+mom2).Pt());
2300                 fhConvDispersionMCAntiProton  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2301                 fhConvM02MCAntiProton         ->Fill( calo->GetM02(), calo2->GetM02());
2302               }
2303             }
2304
2305             //Pairs coming from fragmenting pairs.
2306             if(ancPDG < 22 && ancLabel > 7 && (ancStatus == 11 || ancStatus == 12) ){
2307               fhConvDeltaEtaMCString    ->Fill( pairM, dPhi);
2308               fhConvDeltaPhiMCString    ->Fill( pairM, dPhi);
2309               fhConvAsymMCString        ->Fill( pairM, TMath::Abs(mom.E()-mom2.E())/(mom.E()+mom2.E()) );
2310               fhConvDeltaEtaPhiMCString ->Fill( dPhi, dPhi );
2311               fhConvPtMCString          ->Fill( pairM, (mom+mom2).Pt());
2312               fhConvDispersionMCString  ->Fill( calo->GetDispersion(), calo2->GetDispersion());
2313               fhConvM02MCString         ->Fill( calo->GetM02(), calo2->GetM02());
2314             }
2315
2316           }// Data MC
2317
2318           break;
2319         }
2320
2321       }//Mass loop
2322
2323       //..........................................................................................................
2324       //Pair selected as converted, remove both clusters or recombine them into a photon and put them in the AOD
2325       if(bConverted){
2326         //Add to AOD
2327         if(fAddConvertedPairsToAOD){
2328           //Create AOD of pair analysis
2329           TLorentzVector mpair = mom+mom2;
2330           AliAODPWG4Particle aodpair = AliAODPWG4Particle(mpair);
2331           aodpair.SetLabel(aodph.GetLabel());
2332           //aodpair.SetInputFileIndex(input);
2333
2334           //printf("Index %d, Id %d\n",icalo, calo->GetID());
2335           //Set the indeces of the original caloclusters
2336           aodpair.SetCaloLabel(calo->GetID(),id2);
2337           aodpair.SetDetector(fCalorimeter);
2338           aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
2339           aodpair.SetTag(aodph.GetTag());
2340           aodpair.SetTagged(kTRUE);
2341           //Add AOD with pair object to aod branch
2342           AddAODParticle(aodpair);
2343           //printf("\t \t both added pair\n");
2344         }
2345
2346         //Do not add the current calocluster
2347         if(fRemoveConvertedPair) continue;
2348         else {
2349           //printf("TAGGED\n");
2350           //Tag this cluster as likely conversion
2351           aodph.SetTagged(kTRUE);
2352         }
2353       }//converted pair
2354     }//check conversion
2355     //printf("\t \t added single cluster %d\n",icalo);
2356
2357     //FIXME, this to MakeAnalysisFillHistograms ...
2358     fhNCellsE->Fill(aodph.E(),calo->GetNCells());
2359
2360     //Add AOD with photon object to aod branch
2361     AddAODParticle(aodph);
2362
2363   }//loop
2364
2365   delete [] indexConverted;
2366
2367   if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
2368
2369 }
2370
2371 //__________________________________________________________________
2372 void  AliAnaPhoton::MakeAnalysisFillHistograms()
2373 {
2374   //Fill histograms
2375
2376   //-------------------------------------------------------------------
2377   // Access MC information in stack if requested, check that it exists.
2378   AliStack         * stack       = 0x0;
2379   TParticle        * primary     = 0x0;
2380   TClonesArray     * mcparticles = 0x0;
2381   AliAODMCParticle * aodprimary  = 0x0;
2382
2383   if(IsDataMC()){
2384
2385     if(GetReader()->ReadStack()){
2386       stack =  GetMCStack() ;
2387       if(!stack) {
2388         printf("AliAnaPhoton::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
2389         abort();
2390       }
2391
2392     }
2393     else if(GetReader()->ReadAODMCParticles()){
2394
2395       //Get the list of MC particles
2396       mcparticles = GetReader()->GetAODMCParticles(0);
2397       if(!mcparticles && GetDebug() > 0)        {
2398         printf("AliAnaPhoton::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
2399       }
2400     }
2401   }// is data and MC
2402
2403
2404   // Get vertex
2405   Double_t v[3] = {0,0,0}; //vertex ;
2406   GetReader()->GetVertex(v);
2407   //fhVertex->Fill(v[0],v[1],v[2]);
2408   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
2409
2410   //----------------------------------
2411   //Loop on stored AOD photons
2412   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
2413   if(GetDebug() > 0) printf("AliAnaPhoton::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
2414
2415   for(Int_t iaod = 0; iaod < naod ; iaod++){
2416     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
2417     Int_t pdg = ph->GetIdentifiedParticleType();
2418
2419     if(GetDebug() > 3)
2420       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
2421
2422     //If PID used, fill histos with photons in Calorimeter fCalorimeter
2423     if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue;
2424     if(ph->GetDetector() != fCalorimeter) continue;
2425
2426     if(GetDebug() > 2)
2427       printf("AliAnaPhoton::MakeAnalysisFillHistograms() - ID Photon: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
2428
2429     //................................
2430     //Fill photon histograms
2431     Float_t ptcluster  = ph->Pt();
2432     Float_t phicluster = ph->Phi();
2433     Float_t etacluster = ph->Eta();
2434     Float_t ecluster   = ph->E();
2435
2436     fhEPhoton   ->Fill(ecluster);
2437     fhPtPhoton  ->Fill(ptcluster);
2438     fhPhiPhoton ->Fill(ptcluster,phicluster);
2439     fhEtaPhoton ->Fill(ptcluster,etacluster);
2440     if     (ecluster > 0.5)   fhEtaPhiPhoton  ->Fill(etacluster, phicluster);
2441     else if(GetMinPt() < 0.5) fhEtaPhi05Photon->Fill(etacluster, phicluster);
2442
2443     if(fCheckConversion &&ph->IsTagged()){
2444       fhPtPhotonConv->Fill(ptcluster);
2445       if(ecluster > 0.5)        fhEtaPhiPhotonConv  ->Fill(etacluster, phicluster);
2446       else if(GetMinPt() < 0.5) fhEtaPhi05PhotonConv->Fill(etacluster, phicluster);
2447     }
2448
2449     //.......................................
2450     //Play with the MC data if available
2451     if(IsDataMC()){
2452
2453       FillAcceptanceHistograms();
2454
2455       Int_t tag =ph->GetTag();
2456
2457       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
2458       {
2459         fhMCE  [mcPhoton] ->Fill(ecluster);
2460         fhPtMC [mcPhoton] ->Fill(ptcluster);
2461         fhPhiMC[mcPhoton] ->Fill(ecluster,phicluster);
2462         fhEtaMC[mcPhoton] ->Fill(ecluster,etacluster);
2463
2464         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion))
2465         {
2466           fhMCE  [mcConversion] ->Fill(ecluster);
2467           fhPtMC [mcConversion] ->Fill(ptcluster);
2468           fhPhiMC[mcConversion] ->Fill(ecluster,phicluster);
2469           fhEtaMC[mcConversion] ->Fill(ecluster,etacluster);
2470
2471           if(fCheckConversion){
2472             if(ph->IsTagged()) fhPtConversionTagged ->Fill(ptcluster);
2473             if(ptcluster > 0.5)fhEtaPhiConversion   ->Fill(etacluster,phicluster);
2474             else               fhEtaPhi05Conversion ->Fill(etacluster,phicluster);
2475           }
2476         }
2477
2478         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPrompt)){
2479           fhMCE  [mcPrompt] ->Fill(ecluster);
2480           fhPtMC [mcPrompt] ->Fill(ptcluster);
2481           fhPhiMC[mcPrompt] ->Fill(ecluster,phicluster);
2482           fhEtaMC[mcPrompt] ->Fill(ecluster,etacluster);
2483         }
2484         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCFragmentation))
2485         {
2486           fhMCE  [mcFragmentation] ->Fill(ecluster);
2487           fhPtMC [mcFragmentation] ->Fill(ptcluster);
2488           fhPhiMC[mcFragmentation] ->Fill(ecluster,phicluster);
2489           fhEtaMC[mcFragmentation] ->Fill(ecluster,etacluster);
2490
2491         }
2492         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCISR))
2493         {
2494           fhMCE  [mcISR] ->Fill(ecluster);
2495           fhPtMC [mcISR] ->Fill(ptcluster);
2496           fhPhiMC[mcISR] ->Fill(ecluster,phicluster);
2497           fhEtaMC[mcISR] ->Fill(ecluster,etacluster);
2498         }
2499         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
2500                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2501         {
2502           fhMCE  [mcPi0Decay] ->Fill(ecluster);
2503           fhPtMC [mcPi0Decay] ->Fill(ptcluster);
2504           fhPhiMC[mcPi0Decay] ->Fill(ecluster,phicluster);
2505           fhEtaMC[mcPi0Decay] ->Fill(ecluster,etacluster);
2506         }
2507         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
2508                 GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay))
2509         {
2510           fhMCE  [mcOtherDecay] ->Fill(ecluster);
2511           fhPtMC [mcOtherDecay] ->Fill(ptcluster);
2512           fhPhiMC[mcOtherDecay] ->Fill(ecluster,phicluster);
2513           fhEtaMC[mcOtherDecay] ->Fill(ecluster,etacluster);
2514         }
2515         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0))
2516         {
2517           fhMCE  [mcPi0] ->Fill(ecluster);
2518           fhPtMC [mcPi0] ->Fill(ptcluster);
2519           fhPhiMC[mcPi0] ->Fill(ecluster,phicluster);
2520           fhEtaMC[mcPi0] ->Fill(ecluster,etacluster);
2521         }
2522       }
2523       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron))
2524       {
2525         fhMCE  [mcAntiNeutron] ->Fill(ecluster);
2526         fhPtMC [mcAntiNeutron] ->Fill(ptcluster);
2527         fhPhiMC[mcAntiNeutron] ->Fill(ecluster,phicluster);
2528         fhEtaMC[mcAntiNeutron] ->Fill(ecluster,etacluster);
2529         if(ph->IsTagged() && fCheckConversion) fhPtAntiNeutronTagged ->Fill(ptcluster);
2530
2531       }
2532       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton))
2533       {
2534         fhMCE  [mcAntiProton] ->Fill(ecluster);
2535         fhPtMC [mcAntiProton] ->Fill(ptcluster);
2536         fhPhiMC[mcAntiProton] ->Fill(ecluster,phicluster);
2537         fhEtaMC[mcAntiProton] ->Fill(ecluster,etacluster);
2538         if(ph->IsTagged() && fCheckConversion) fhPtAntiProtonTagged ->Fill(ptcluster);
2539       }
2540       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron))
2541       {
2542         fhMCE  [mcElectron] ->Fill(ecluster);
2543         fhPtMC [mcElectron] ->Fill(ptcluster);
2544         fhPhiMC[mcElectron] ->Fill(ecluster,phicluster);
2545         fhEtaMC[mcElectron] ->Fill(ecluster,etacluster);
2546       }
2547       else{
2548         fhMCE  [mcOther] ->Fill(ecluster);
2549         fhPtMC [mcOther] ->Fill(ptcluster);
2550         fhPhiMC[mcOther] ->Fill(ecluster,phicluster);
2551         fhEtaMC[mcOther] ->Fill(ecluster,etacluster);
2552         if(ph->IsTagged() && fCheckConversion) fhPtUnknownTagged ->Fill(ptcluster);
2553
2554
2555         //               printf(" AliAnaPhoton::MakeAnalysisFillHistograms() - Label %d, pT %2.3f Unknown, bits set: ",
2556         //                                      ph->GetLabel(),ph->Pt());
2557         //                for(Int_t i = 0; i < 20; i++) {
2558         //                        if(GetMCAnalysisUtils()->CheckTagBit(tag,i)) printf(" %d, ",i);
2559         //                }
2560         //                printf("\n");
2561
2562       }
2563
2564       //....................................................................
2565       // Access MC information in stack if requested, check that it exists.
2566       Int_t label =ph->GetLabel();
2567       if(label < 0) {
2568         if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
2569         continue;
2570       }
2571
2572       Float_t eprim   = 0;
2573       Float_t ptprim  = 0;
2574       if(GetReader()->ReadStack()){
2575
2576         if(label >=  stack->GetNtrack()) {
2577           if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
2578           continue ;
2579         }
2580
2581         primary = stack->Particle(label);
2582         if(!primary){
2583           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
2584           continue;
2585         }
2586         eprim   = primary->Energy();
2587         ptprim  = primary->Pt();
2588
2589       }
2590       else if(GetReader()->ReadAODMCParticles()){
2591         //Check which is the input
2592         if(ph->GetInputFileIndex() == 0){
2593           if(!mcparticles) continue;
2594           if(label >=  mcparticles->GetEntriesFast()) {
2595             if(GetDebug() > 2)  printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
2596                                        label, mcparticles->GetEntriesFast());
2597             continue ;
2598           }
2599           //Get the particle
2600           aodprimary = (AliAODMCParticle*) mcparticles->At(label);
2601
2602         }
2603
2604         if(!aodprimary){
2605           printf("AliAnaPhoton::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
2606           continue;
2607         }
2608
2609         eprim   = aodprimary->E();
2610         ptprim  = aodprimary->Pt();
2611
2612       }
2613
2614       fh2E     ->Fill(ecluster, eprim);
2615       fh2Pt    ->Fill(ptcluster, ptprim);
2616       fhDeltaE ->Fill(eprim-ecluster);
2617       fhDeltaPt->Fill(ptprim-ptcluster);
2618       if(eprim > 0)  fhRatioE  ->Fill(ecluster/eprim);
2619       if(ptprim > 0) fhRatioPt ->Fill(ptcluster/ptprim);
2620
2621     }//Histograms with MC
2622
2623   }// aod loop
2624
2625 }
2626
2627
2628 //__________________________________________________________________
2629 void AliAnaPhoton::Print(const Option_t * opt) const
2630 {
2631   //Print some relevant parameters set for the analysis
2632
2633   if(! opt)
2634     return;
2635
2636   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
2637   AliAnaPartCorrBaseClass::Print(" ");
2638
2639   printf("Calorimeter            =     %s\n", fCalorimeter.Data()) ;
2640   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
2641   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
2642   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
2643   printf("Reject clusters with a track matched = %d\n",fRejectTrackMatch);
2644   printf("Check Pair Conversion                = %d\n",fCheckConversion);
2645   printf("Add conversion pair to AOD           = %d\n",fAddConvertedPairsToAOD);
2646   printf("Conversion pair mass cut             = %f\n",fMassCut);
2647   printf("Conversion selection cut : A < %1.2f; %1.3f < Dphi < %1.3f; Deta < %1.3f\n",
2648          fConvAsymCut,fConvDPhiMinCut, fConvDPhiMaxCut, fConvDEtaCut);
2649   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
2650   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
2651   printf("    \n") ;
2652
2653 }