PWGGA/CaloTrackCorrelations/AliAnaElectron.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
  16 //_________________________________________________________________________
  17 //
  18 // Class for the photon identification.
  19 // Clusters from calorimeters are identified as photons
  20 // and kept in the AOD. Few histograms produced.
  21 // Copy of AliAnaPhoton just add electron id.
  22 //
  23 // -- Author: Gustavo Conesa (LPSC-IN2P3-CRNS)
  24 //////////////////////////////////////////////////////////////////////////////
  25
  26
  27 // --- ROOT system ---
  28 #include <TH2F.h>
  29 #include <TH3D.h>
  30 #include <TClonesArray.h>
  31 #include <TObjString.h>
  32 //#include <Riostream.h>
  33 #include "TParticle.h"
  34 #include "TDatabasePDG.h"
  35 #include "AliVTrack.h"
  36
  37 // --- Analysis system ---
  38 #include "AliAnaElectron.h"
  39 #include "AliCaloTrackReader.h"
  40 #include "AliStack.h"
  41 #include "AliCaloPID.h"
  42 #include "AliMCAnalysisUtils.h"
  43 #include "AliFiducialCut.h"
  44 #include "AliVCluster.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliMixedEvent.h"
  47
  48
  49 ClassImp(AliAnaElectron)
  50
  51 //________________________________
  52 AliAnaElectron::AliAnaElectron() :
  53     AliAnaCaloTrackCorrBaseClass(),
  54     fMinDist(0.),                         fMinDist2(0.),                         fMinDist3(0.),
  55     fTimeCutMin(-1),                      fTimeCutMax(999999),
  56     fNCellsCut(0),                        fNLMCutMin(-1),                        fNLMCutMax(10),
  57     fFillSSHistograms(kFALSE),             fFillOnlySimpleSSHisto(1),
  58     fFillWeightHistograms(kFALSE),        fNOriginHistograms(8),
  59     fdEdxMin(0.),                         fdEdxMax (200.),
  60     fEOverPMin(0),                        fEOverPMax (2),
  61     fAODParticle(0),
  62     // Histograms
  63     fhdEdxvsE(0),                         fhdEdxvsP(0),
  64     fhEOverPvsE(0),                       fhEOverPvsP(0),
  65     fhdEdxvsECutM02(0),                   fhdEdxvsPCutM02(0),
  66     fhEOverPvsECutM02(0),                 fhEOverPvsPCutM02(0),
  67     fhdEdxvsECutEOverP(0),                fhdEdxvsPCutEOverP(0),
  68     fhEOverPvsECutM02CutdEdx(0),          fhEOverPvsPCutM02CutdEdx(0),
  69     // Weight studies
  70     fhECellClusterRatio(0),               fhECellClusterLogRatio(0),
  71     fhEMaxCellClusterRatio(0),            fhEMaxCellClusterLogRatio(0),
  72     // MC histograms
  73     // Electron SS MC histograms
  74     fhMCElectronELambda0NoOverlap(0),
  75     fhMCElectronELambda0TwoOverlap(0),    fhMCElectronELambda0NOverlap(0),
  76
  77     //Embedding
  78     fhEmbeddedSignalFractionEnergy(0),
  79     fhEmbedElectronELambda0FullSignal(0), fhEmbedElectronELambda0MostlySignal(0),
  80     fhEmbedElectronELambda0MostlyBkg(0),  fhEmbedElectronELambda0FullBkg(0)
  81 {
  82   //default ctor
  83   for(Int_t index = 0; index < 2; index++)
  84   {
  85     fhNCellsE [index] = 0;
  86     fhNLME    [index] = 0;
  87     fhTimeE   [index] = 0;
  88     fhMaxCellDiffClusterE[index] = 0;
  89     fhE       [index] = 0;
  90     fhPt      [index] = 0;
  91     fhPhi     [index] = 0;
  92     fhEta     [index] = 0;
  93     fhEtaPhi  [index] = 0;
  94     fhEtaPhi05[index] = 0;
  95
  96     // Shower shape histograms
  97     fhDispE   [index] = 0;
  98     fhLam0E   [index] = 0;
  99     fhLam1E   [index] = 0;
 100     fhDispETRD[index] = 0;
 101     fhLam0ETRD[index] = 0;
 102     fhLam1ETRD[index] = 0;
 103     fhNCellsLam0LowE [index] = 0;
 104     fhNCellsLam0HighE[index] = 0;
 105     fhEtaLam0LowE    [index] = 0;
 106     fhPhiLam0LowE    [index] = 0;
 107     fhEtaLam0HighE   [index] = 0;
 108     fhPhiLam0HighE   [index] = 0;
 109
 110     fhDispEtaE       [index] = 0;
 111     fhDispPhiE       [index] = 0;
 112     fhSumEtaE        [index] = 0;
 113     fhSumPhiE        [index] = 0;
 114     fhSumEtaPhiE     [index] = 0;
 115     fhDispEtaPhiDiffE[index] = 0;
 116     fhSphericityE    [index] = 0;
 117
 118     for(Int_t i = 0; i < 10; i++)
 119     {
 120       fhMCPt     [index][i] = 0;
 121       fhMCE      [index][i] = 0;
 122       fhMCPhi    [index][i] = 0;
 123       fhMCEta    [index][i] = 0;
 124       fhMCDeltaE [index][i] = 0;
 125       fhMC2E     [index][i] = 0;
 126       fhMCdEdxvsE       [i] = 0;
 127       fhMCdEdxvsP       [i] = 0;
 128       fhMCEOverPvsE     [i] = 0;
 129       fhMCEOverPvsP     [i] = 0;
 130     }
 131
 132     for(Int_t i = 0; i < 6; i++)
 133     {
 134       fhMCELambda0       [index][i] = 0;
 135       fhMCEDispEta       [index][i] = 0;
 136       fhMCEDispPhi       [index][i] = 0;
 137       fhMCESumEtaPhi     [index][i] = 0;
 138       fhMCEDispEtaPhiDiff[index][i] = 0;
 139       fhMCESphericity    [index][i] = 0;
 140     }
 141
 142     for(Int_t i = 0; i < 5; i++)
 143     {
 144       fhDispEtaDispPhiEBin[index][i] = 0 ;
 145     }
 146   }
 147
 148   //Weight studies
 149   for(Int_t i =0; i < 14; i++)
 150   {
 151     fhLambda0ForW0[i] = 0;
 152     //fhLambda1ForW0[i] = 0;
 153   }
 154
 155   //Initialize parameters
 156   InitParameters();
 157
 158 }
 159
 160 //____________________________________________________________________________
 161 Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
 162 {
 163   //Select clusters if they pass different cuts
 164   if(GetDebug() > 2)
 165     printf("AliAnaElectron::ClusterSelected() Current Event %d; Before selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 166            GetReader()->GetEventNumber(),
 167            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 168
 169   //.......................................
 170   //If too small or big energy, skip it
 171   if(mom.E() < GetMinEnergy() || mom.E() > GetMaxEnergy() ) return kFALSE ;
 172   if(GetDebug() > 2) printf("\t Cluster %d Pass E Cut \n",calo->GetID());
 173
 174   //.......................................
 175   // TOF cut, BE CAREFUL WITH THIS CUT
 176   Double_t tof = calo->GetTOF()*1e9;
 177   if(tof < fTimeCutMin || tof > fTimeCutMax) return kFALSE;
 178   if(GetDebug() > 2)  printf("\t Cluster %d Pass Time Cut \n",calo->GetID());
 179
 180   //.......................................
 181   if(calo->GetNCells() <= fNCellsCut && GetReader()->GetDataType() != AliCaloTrackReader::kMC) return kFALSE;
 182   if(GetDebug() > 2) printf("\t Cluster %d Pass NCell Cut \n",calo->GetID());
 183
 184   //.......................................
 185   //Check acceptance selection
 186   if(IsFiducialCutOn()){
 187     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,GetCalorimeter()) ;
 188     if(! in ) return kFALSE ;
 189   }
 190   if(GetDebug() > 2) printf("Fiducial cut passed \n");
 191
 192   //.......................................
 193   //Skip not matched clusters with tracks
 194   if(!IsTrackMatched(calo, GetReader()->GetInputEvent())) {
 195       if(GetDebug() > 2) printf("\t Reject non track-matched clusters\n");
 196       return kFALSE ;
 197   }
 198   else if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
 199
 200   //...........................................
 201   // skip clusters with too many maxima
 202   if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
 203   if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
 204
 205   //.......................................
 206   //Check Distance to Bad channel, set bit.
 207   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
 208   if(distBad < 0.) distBad=9999. ; //workout strange convension dist = -1. ;
 209   if(distBad < fMinDist) {//In bad channel (PHOS cristal size 2.2x2.2 cm), EMCAL ( cell units )
 210     return kFALSE ;
 211   }
 212   else if(GetDebug() > 2) printf("\t Bad channel cut passed %4.2f > %2.2f \n",distBad, fMinDist);
 213   //printf("Cluster %d Pass Bad Dist Cut \n",icalo);
 214
 215   if(GetDebug() > 0)
 216     printf("AliAnaElectron::ClusterSelected() Current Event %d; After  selection : E %2.2f, pT %2.2f, Ecl %2.2f, phi %2.2f, eta %2.2f\n",
 217            GetReader()->GetEventNumber(),
 218            mom.E(), mom.Pt(),calo->E(),mom.Phi()*TMath::RadToDeg(),mom.Eta());
 219
 220   //All checks passed, cluster selected
 221   return kTRUE;
 222
 223 }
 224
 225 //______________________________________________________________________________________________
 226 void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
 227 {
 228
 229   //Fill cluster Shower Shape histograms
 230
 231   if(!fFillSSHistograms || GetMixedEvent()) return;
 232
 233   Int_t pidIndex = 0;// Electron
 234   if     (pidTag == AliCaloPID::kElectron)      pidIndex = 0;
 235   else if(pidTag == AliCaloPID::kChargedHadron) pidIndex = 1;
 236   else return;
 237
 238   Float_t energy  = cluster->E();
 239   Int_t   ncells  = cluster->GetNCells();
 240   Float_t lambda0 = cluster->GetM02();
 241   Float_t lambda1 = cluster->GetM20();
 242   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
 243
 244   Float_t l0   = 0., l1   = 0.;
 245   Float_t dispp= 0., dEta = 0., dPhi    = 0.;
 246   Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;
 247
 248   TLorentzVector mom;
 249   if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
 250   {
 251     cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
 252   else{
 253     Double_t vertex[]={0,0,0};
 254     cluster->GetMomentum(mom,vertex) ;
 255   }
 256
 257   Float_t eta = mom.Eta();
 258   Float_t phi = mom.Phi();
 259   if(phi < 0) phi+=TMath::TwoPi();
 260
 261   fhLam0E[pidIndex] ->Fill(energy,lambda0);
 262   fhLam1E[pidIndex] ->Fill(energy,lambda1);
 263   fhDispE[pidIndex] ->Fill(energy,disp);
 264
 265   if(GetCalorimeter() == "EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
 266      GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD() )
 267   {
 268     fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
 269     fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
 270     fhDispETRD[pidIndex]->Fill(energy,disp);
 271   }
 272
 273   if(!fFillOnlySimpleSSHisto)
 274   {
 275     if(energy < 2)
 276     {
 277       fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
 278       fhEtaLam0LowE[pidIndex]    ->Fill(eta,   lambda0);
 279       fhPhiLam0LowE[pidIndex]    ->Fill(phi,   lambda0);
 280     }
 281     else
 282     {
 283       fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
 284       fhEtaLam0HighE[pidIndex]   ->Fill(eta,   lambda0);
 285       fhPhiLam0HighE[pidIndex]   ->Fill(phi,   lambda0);
 286     }
 287
 288     if(GetCalorimeter() == "EMCAL")
 289     {
 290       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
 291                                                                                    l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
 292       fhDispEtaE        [pidIndex]-> Fill(energy,dEta);
 293       fhDispPhiE        [pidIndex]-> Fill(energy,dPhi);
 294       fhSumEtaE         [pidIndex]-> Fill(energy,sEta);
 295       fhSumPhiE         [pidIndex]-> Fill(energy,sPhi);
 296       fhSumEtaPhiE      [pidIndex]-> Fill(energy,sEtaPhi);
 297       fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
 298       if(dEta+dPhi>0)fhSphericityE     [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
 299
 300       if      (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
 301       else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
 302       else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
 303       else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
 304       else                  fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
 305
 306     }
 307   }
 308
 309   if(IsDataMC())
 310   {
 311     AliVCaloCells* cells = 0;
 312     if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
 313     else                        cells = GetPHOSCells();
 314
 315     //Fill histograms to check shape of embedded clusters
 316     Float_t fraction = 0;
 317     if(GetReader()->IsEmbeddedClusterSelectionOn()){//Only working for EMCAL
 318
 319       Float_t clusterE = 0; // recalculate in case corrections applied.
 320       Float_t cellE    = 0;
 321       for(Int_t icell = 0; icell < cluster->GetNCells(); icell++){
 322         cellE    = cells->GetCellAmplitude(cluster->GetCellAbsId(icell));
 323         clusterE+=cellE;
 324         fraction+=cellE*cluster->GetCellAmplitudeFraction(icell);
 325       }
 326
 327       //Fraction of total energy due to the embedded signal
 328       fraction/=clusterE;
 329
 330       if(GetDebug() > 1 ) printf("AliAnaElectron::FillShowerShapeHistogram() - Energy fraction of embedded signal %2.3f, Energy %2.3f\n",fraction, clusterE);
 331
 332       fhEmbeddedSignalFractionEnergy->Fill(clusterE,fraction);
 333
 334     }  // embedded fraction
 335
 336     // Check the origin and fill histograms
 337     Int_t index = -1;
 338
 339     if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
 340        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
 341        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
 342        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
 343     {
 344       index = kmcssPhoton;
 345
 346     }//photon   no conversion
 347     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron &&
 348               !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
 349     {
 350       index = kmcssElectron;
 351
 352       if(!GetReader()->IsEmbeddedClusterSelectionOn())
 353       {
 354         //Check particle overlaps in cluster
 355
 356         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
 357         Int_t ancPDG = 0, ancStatus = -1;
 358         TLorentzVector momentum; TVector3 prodVertex;
 359         Int_t ancLabel = 0;
 360         Int_t noverlaps = 1;
 361         for (UInt_t ilab = 0; ilab < cluster->GetNLabels(); ilab++ ) {
 362           ancLabel = GetMCAnalysisUtils()->CheckCommonAncestor(cluster->GetLabels()[0],cluster->GetLabels()[ilab], GetReader(),ancPDG,ancStatus,momentum,prodVertex);
 363           if(ancPDG!=22 && TMath::Abs(ancPDG)!=11) noverlaps++;
 364         }
 365
 366         if(noverlaps == 1){
 367           fhMCElectronELambda0NoOverlap  ->Fill(energy, lambda0);
 368         }
 369         else if(noverlaps == 2){
 370           fhMCElectronELambda0TwoOverlap ->Fill(energy, lambda0);
 371         }
 372         else if(noverlaps > 2){
 373           fhMCElectronELambda0NOverlap   ->Fill(energy, lambda0);
 374         }
 375         else {
 376           printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d for ancestor %d!!", noverlaps, ancLabel);
 377         }
 378       }//No embedding
 379
 380       //Fill histograms to check shape of embedded clusters
 381       if(GetReader()->IsEmbeddedClusterSelectionOn())
 382       {
 383         if     (fraction > 0.9)
 384         {
 385           fhEmbedElectronELambda0FullSignal   ->Fill(energy, lambda0);
 386         }
 387         else if(fraction > 0.5)
 388         {
 389           fhEmbedElectronELambda0MostlySignal ->Fill(energy, lambda0);
 390         }
 391         else if(fraction > 0.1)
 392         {
 393           fhEmbedElectronELambda0MostlyBkg    ->Fill(energy, lambda0);
 394         }
 395         else
 396         {
 397           fhEmbedElectronELambda0FullBkg      ->Fill(energy, lambda0);
 398         }
 399       } // embedded
 400     }//electron
 401     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) &&
 402                GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
 403     {
 404       index = kmcssConversion;
 405     }//conversion photon
 406     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  )
 407     {
 408       index = kmcssPi0;
 409     }//pi0
 410     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  )
 411     {
 412       index = kmcssEta;
 413     }//eta
 414     else
 415     {
 416       index = kmcssOther;
 417     }//other particles
 418
 419     fhMCELambda0[pidIndex][index]    ->Fill(energy, lambda0);
 420
 421     if(GetCalorimeter() == "EMCAL" && !fFillOnlySimpleSSHisto)
 422     {
 423       fhMCEDispEta        [pidIndex][index]-> Fill(energy,dEta);
 424       fhMCEDispPhi        [pidIndex][index]-> Fill(energy,dPhi);
 425       fhMCESumEtaPhi      [pidIndex][index]-> Fill(energy,sEtaPhi);
 426       fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
 427       if(dEta+dPhi>0)fhMCESphericity     [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
 428     }
 429
 430   }//MC data
 431
 432 }
 433
 434 //_____________________________________________
 435 TObjString *  AliAnaElectron::GetAnalysisCuts()
 436 {
 437   //Save parameters used for analysis
 438   TString parList ; //this will be list of parameters used for this analysis.
 439   const Int_t buffersize = 255;
 440   char onePar[buffersize] ;
 441
 442   snprintf(onePar,buffersize,"--- AliAnaElectron ---\n") ;
 443   parList+=onePar ;
 444   snprintf(onePar,buffersize,"Calorimeter: %s\n",GetCalorimeter().Data()) ;
 445   parList+=onePar ;
 446   snprintf(onePar,buffersize," %2.2f < dEdx < %2.2f  \n",fdEdxMin,fdEdxMax) ;
 447   parList+=onePar ;
 448   snprintf(onePar,buffersize," %2.2f <  E/P < %2.2f  \n",fEOverPMin, fEOverPMax) ;
 449   parList+=onePar ;
 450   snprintf(onePar,buffersize,"fMinDist =%2.2f (Minimal distance to bad channel to accept cluster) \n",fMinDist) ;
 451   parList+=onePar ;
 452   snprintf(onePar,buffersize,"fMinDist2=%2.2f (Cuts on Minimal distance to study acceptance evaluation) \n",fMinDist2) ;
 453   parList+=onePar ;
 454   snprintf(onePar,buffersize,"fMinDist3=%2.2f (One more cut on distance used for acceptance-efficiency study) \n",fMinDist3) ;
 455   parList+=onePar ;
 456
 457   //Get parameters set in base class.
 458   parList += GetBaseParametersList() ;
 459
 460   //Get parameters set in PID class.
 461   parList += GetCaloPID()->GetPIDParametersList() ;
 462
 463   //Get parameters set in FiducialCut class (not available yet)
 464   //parlist += GetFidCut()->GetFidCutParametersList()
 465
 466   return new TObjString(parList) ;
 467 }
 468
 469 //_______________________________________________
 470 TList *  AliAnaElectron::GetCreateOutputObjects()
 471 {
 472   // Create histograms to be saved in output file and
 473   // store them in outputContainer
 474   TList * outputContainer = new TList() ;
 475   outputContainer->SetName("ElectronHistos") ;
 476
 477   Int_t nptbins     = GetHistogramRanges()->GetHistoPtBins();           Float_t ptmax     = GetHistogramRanges()->GetHistoPtMax();           Float_t ptmin     = GetHistogramRanges()->GetHistoPtMin();
 478   Int_t nphibins    = GetHistogramRanges()->GetHistoPhiBins();          Float_t phimax    = GetHistogramRanges()->GetHistoPhiMax();          Float_t phimin    = GetHistogramRanges()->GetHistoPhiMin();
 479   Int_t netabins    = GetHistogramRanges()->GetHistoEtaBins();          Float_t etamax    = GetHistogramRanges()->GetHistoEtaMax();          Float_t etamin    = GetHistogramRanges()->GetHistoEtaMin();
 480   Int_t ssbins      = GetHistogramRanges()->GetHistoShowerShapeBins();  Float_t ssmax     = GetHistogramRanges()->GetHistoShowerShapeMax();  Float_t ssmin     = GetHistogramRanges()->GetHistoShowerShapeMin();
 481   Int_t nbins       = GetHistogramRanges()->GetHistoNClusterCellBins(); Int_t   nmax      = GetHistogramRanges()->GetHistoNClusterCellMax(); Int_t   nmin      = GetHistogramRanges()->GetHistoNClusterCellMin();
 482   Int_t ndedxbins   = GetHistogramRanges()->GetHistodEdxBins();         Float_t dedxmax   = GetHistogramRanges()->GetHistodEdxMax();         Float_t dedxmin   = GetHistogramRanges()->GetHistodEdxMin();
 483   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
 484   Int_t tbins       = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax      = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin      = GetHistogramRanges()->GetHistoTimeMin();
 485
 486
 487   // MC labels, titles, for originator particles
 488   TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
 489   TString pnamess[] = { "Photon","Hadron" ,"Pi0"    ,"Eta" ,"Conversion"     ,"Electron"} ;
 490   TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
 491     "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P"                    } ;
 492
 493   TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
 494     "Conversion", "Hadron", "AntiNeutron","AntiProton"                        } ;
 495
 496
 497   fhdEdxvsE  = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 498   fhdEdxvsE->SetXTitle("E (GeV)");
 499   fhdEdxvsE->SetYTitle("<dE/dx>");
 500   outputContainer->Add(fhdEdxvsE);
 501
 502   fhdEdxvsP  = new TH2F ("hdEdxvsP","matched track <dE/dx> vs track P ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 503   fhdEdxvsP->SetXTitle("P (GeV/c)");
 504   fhdEdxvsP->SetYTitle("<dE/dx>");
 505   outputContainer->Add(fhdEdxvsP);
 506
 507   fhEOverPvsE  = new TH2F ("hEOverPvsE","matched track E/p vs cluster E ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 508   fhEOverPvsE->SetXTitle("E (GeV)");
 509   fhEOverPvsE->SetYTitle("E/p");
 510   outputContainer->Add(fhEOverPvsE);
 511
 512   fhEOverPvsP  = new TH2F ("hEOverPvsP","matched track E/p vs track P ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 513   fhEOverPvsP->SetXTitle("P (GeV/c)");
 514   fhEOverPvsP->SetYTitle("E/p");
 515   outputContainer->Add(fhEOverPvsP);
 516
 517
 518   fhdEdxvsECutM02  = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 519   fhdEdxvsECutM02->SetXTitle("E (GeV)");
 520   fhdEdxvsECutM02->SetYTitle("<dE/dx>");
 521   outputContainer->Add(fhdEdxvsECutM02);
 522
 523   fhdEdxvsPCutM02  = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 524   fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
 525   fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
 526   outputContainer->Add(fhdEdxvsPCutM02);
 527
 528   fhEOverPvsECutM02  = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 529   fhEOverPvsECutM02->SetXTitle("E (GeV)");
 530   fhEOverPvsECutM02->SetYTitle("E/p");
 531   outputContainer->Add(fhEOverPvsECutM02);
 532
 533   fhEOverPvsPCutM02  = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 534   fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
 535   fhEOverPvsPCutM02->SetYTitle("E/p");
 536   outputContainer->Add(fhEOverPvsPCutM02);
 537
 538
 539   fhdEdxvsECutEOverP  = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 540   fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
 541   fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
 542   outputContainer->Add(fhdEdxvsECutEOverP);
 543
 544   fhdEdxvsPCutEOverP  = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 545   fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
 546   fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
 547   outputContainer->Add(fhdEdxvsPCutEOverP);
 548
 549   fhEOverPvsECutM02CutdEdx  = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 550   fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
 551   fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
 552   outputContainer->Add(fhEOverPvsECutM02CutdEdx);
 553
 554   fhEOverPvsPCutM02CutdEdx  = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 555   fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
 556   fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
 557   outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
 558
 559   if(IsDataMC())
 560   {
 561     for(Int_t i = 0; i < fNOriginHistograms; i++)
 562     {
 563       fhMCdEdxvsE[i]  = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
 564                                      Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
 565                                      nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 566       fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
 567       fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
 568       outputContainer->Add(fhMCdEdxvsE[i]) ;
 569
 570       fhMCdEdxvsP[i]  = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
 571                                  Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
 572                                  nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
 573       fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
 574       fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
 575       outputContainer->Add(fhMCdEdxvsP[i]) ;
 576
 577
 578       fhMCEOverPvsE[i]  = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
 579                                  Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
 580                                  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 581       fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
 582       fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
 583       outputContainer->Add(fhMCEOverPvsE[i]) ;
 584
 585       fhMCEOverPvsP[i]  = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
 586                                  Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
 587                                  nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
 588       fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
 589       fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
 590       outputContainer->Add(fhMCEOverPvsP[i]) ;
 591
 592     }
 593   }
 594
 595   TString pidParticle[] = {"Electron","ChargedHadron"} ;
 596
 597   if(fFillWeightHistograms)
 598   {
 599
 600     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
 601                                      nptbins,ptmin,ptmax, 100,0,1.);
 602     fhECellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
 603     fhECellClusterRatio->SetYTitle("E_{cell i}/E_{cluster}");
 604     outputContainer->Add(fhECellClusterRatio);
 605
 606     fhECellClusterLogRatio  = new TH2F ("hECellClusterLogRatio"," Log(cell energy / cluster energy) vs cluster energy, for selected electrons",
 607                                         nptbins,ptmin,ptmax, 100,-10,0);
 608     fhECellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
 609     fhECellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
 610     outputContainer->Add(fhECellClusterLogRatio);
 611
 612     fhEMaxCellClusterRatio  = new TH2F ("hEMaxCellClusterRatio"," max cell energy / cluster energy vs cluster energy, for selected electrons",
 613                                         nptbins,ptmin,ptmax, 100,0,1.);
 614     fhEMaxCellClusterRatio->SetXTitle("E_{cluster} (GeV) ");
 615     fhEMaxCellClusterRatio->SetYTitle("E_{max cell}/E_{cluster}");
 616     outputContainer->Add(fhEMaxCellClusterRatio);
 617
 618     fhEMaxCellClusterLogRatio  = new TH2F ("hEMaxCellClusterLogRatio"," Log(max cell energy / cluster energy) vs cluster energy, for selected electrons",
 619                                            nptbins,ptmin,ptmax, 100,-10,0);
 620     fhEMaxCellClusterLogRatio->SetXTitle("E_{cluster} (GeV) ");
 621     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
 622     outputContainer->Add(fhEMaxCellClusterLogRatio);
 623
 624     for(Int_t iw = 0; iw < 14; iw++)
 625     {
 626       fhLambda0ForW0[iw]  = new TH2F (Form("hLambda0ForW0%d",iw),Form("shower shape, #lambda^{2}_{0} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
 627                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 628       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
 629       fhLambda0ForW0[iw]->SetYTitle("#lambda^{2}_{0}");
 630       outputContainer->Add(fhLambda0ForW0[iw]);
 631
 632       //        fhLambda1ForW0[iw]  = new TH2F (Form("hLambda1ForW0%d",iw),Form("shower shape, #lambda^{2}_{1} vs E, w0 = %1.1f, for selected electrons",1+0.5*iw),
 633       //                                        nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 634       //        fhLambda1ForW0[iw]->SetXTitle("E_{cluster}");
 635       //        fhLambda1ForW0[iw]->SetYTitle("#lambda^{2}_{1}");
 636       //        outputContainer->Add(fhLambda1ForW0[iw]);
 637
 638     }
 639   }
 640
 641   for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
 642   {
 643     //Shower shape
 644     if(fFillSSHistograms)
 645     {
 646       fhLam0E[pidIndex]  = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
 647                                      Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()),
 648                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 649       fhLam0E[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 650       fhLam0E[pidIndex]->SetXTitle("E (GeV)");
 651       outputContainer->Add(fhLam0E[pidIndex]);
 652
 653       fhLam1E[pidIndex]  = new TH2F (Form("h%sLam1E",pidParticle[pidIndex].Data()),
 654                                      Form("%s: #lambda_{1}^{2} vs E",pidParticle[pidIndex].Data()),
 655                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 656       fhLam1E[pidIndex]->SetYTitle("#lambda_{1}^{2}");
 657       fhLam1E[pidIndex]->SetXTitle("E (GeV)");
 658       outputContainer->Add(fhLam1E[pidIndex]);
 659
 660       fhDispE[pidIndex]  = new TH2F (Form("h%sDispE",pidParticle[pidIndex].Data()),
 661                                      Form("%s: dispersion^{2} vs E",pidParticle[pidIndex].Data()),
 662                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 663       fhDispE[pidIndex]->SetYTitle("D^{2}");
 664       fhDispE[pidIndex]->SetXTitle("E (GeV) ");
 665       outputContainer->Add(fhDispE[pidIndex]);
 666
 667       if(GetCalorimeter() == "EMCAL" &&  GetFirstSMCoveredByTRD() >=0 )
 668       {
 669         fhLam0ETRD[pidIndex]  = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
 670                                           Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
 671                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 672         fhLam0ETRD[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 673         fhLam0ETRD[pidIndex]->SetXTitle("E (GeV)");
 674         outputContainer->Add(fhLam0ETRD[pidIndex]);
 675
 676         fhLam1ETRD[pidIndex]  = new TH2F (Form("h%sLam1ETRD",pidParticle[pidIndex].Data()),
 677                                           Form("%s: #lambda_{1}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
 678                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 679         fhLam1ETRD[pidIndex]->SetYTitle("#lambda_{1}^{2}");
 680         fhLam1ETRD[pidIndex]->SetXTitle("E (GeV)");
 681         outputContainer->Add(fhLam1ETRD[pidIndex]);
 682
 683         fhDispETRD[pidIndex]  = new TH2F (Form("h%sDispETRD",pidParticle[pidIndex].Data()),
 684                                           Form("%s: dispersion^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()),
 685                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 686         fhDispETRD[pidIndex]->SetYTitle("Dispersion^{2}");
 687         fhDispETRD[pidIndex]->SetXTitle("E (GeV) ");
 688         outputContainer->Add(fhDispETRD[pidIndex]);
 689       }
 690
 691       if(!fFillOnlySimpleSSHisto)
 692       {
 693
 694         fhNCellsLam0LowE[pidIndex]  = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
 695                                                 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
 696                                                 nbins,nmin, nmax, ssbins,ssmin,ssmax);
 697         fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
 698         fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 699         outputContainer->Add(fhNCellsLam0LowE[pidIndex]);
 700
 701         fhNCellsLam0HighE[pidIndex]  = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
 702                                                  Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
 703                                                  nbins,nmin, nmax, ssbins,ssmin,ssmax);
 704         fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
 705         fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 706         outputContainer->Add(fhNCellsLam0HighE[pidIndex]);
 707
 708
 709         fhEtaLam0LowE[pidIndex]  = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
 710                                              Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
 711                                              netabins,etamin,etamax, ssbins,ssmin,ssmax);
 712         fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 713         fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
 714         outputContainer->Add(fhEtaLam0LowE[pidIndex]);
 715
 716         fhPhiLam0LowE[pidIndex]  = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
 717                                              Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
 718                                              nphibins,phimin,phimax, ssbins,ssmin,ssmax);
 719         fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 720         fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
 721         outputContainer->Add(fhPhiLam0LowE[pidIndex]);
 722
 723         fhEtaLam0HighE[pidIndex]  = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
 724                                               Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
 725                                               netabins,etamin,etamax, ssbins,ssmin,ssmax);
 726         fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 727         fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
 728         outputContainer->Add(fhEtaLam0HighE[pidIndex]);
 729
 730         fhPhiLam0HighE[pidIndex]  = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
 731                                               Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
 732                                               nphibins,phimin,phimax, ssbins,ssmin,ssmax);
 733         fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
 734         fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
 735         outputContainer->Add(fhPhiLam0HighE[pidIndex]);
 736
 737         if(GetCalorimeter() == "EMCAL")
 738         {
 739           fhDispEtaE[pidIndex]  = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
 740                                             Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
 741                                             nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 742           fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
 743           fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
 744           outputContainer->Add(fhDispEtaE[pidIndex]);
 745
 746           fhDispPhiE[pidIndex]  = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
 747                                             Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),
 748                                             nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 749           fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
 750           fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
 751           outputContainer->Add(fhDispPhiE[pidIndex]);
 752
 753           fhSumEtaE[pidIndex]  = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
 754                                            Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),
 755                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 756           fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
 757           fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
 758           outputContainer->Add(fhSumEtaE[pidIndex]);
 759
 760           fhSumPhiE[pidIndex]  = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
 761                                            Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),
 762                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 763           fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
 764           fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
 765           outputContainer->Add(fhSumPhiE[pidIndex]);
 766
 767           fhSumEtaPhiE[pidIndex]  = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
 768                                               Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),
 769                                               nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
 770           fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
 771           fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
 772           outputContainer->Add(fhSumEtaPhiE[pidIndex]);
 773
 774           fhDispEtaPhiDiffE[pidIndex]  = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
 775                                                    Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()),
 776                                                    nptbins,ptmin,ptmax,200, -10,10);
 777           fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
 778           fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
 779           outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);
 780
 781           fhSphericityE[pidIndex]  = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
 782                                                Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),
 783                                                nptbins,ptmin,ptmax, 200, -1,1);
 784           fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
 785           fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
 786           outputContainer->Add(fhSphericityE[pidIndex]);
 787
 788           Int_t bin[] = {0,2,4,6,10,1000};
 789           for(Int_t i = 0; i < 5; i++)
 790           {
 791             fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
 792                                                           Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]),
 793                                                           ssbins,ssmin,ssmax , ssbins,ssmin,ssmax);
 794             fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
 795             fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
 796             outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]);
 797           }
 798         }
 799       }
 800     } // Shower shape
 801
 802     if(IsDataMC())
 803     {
 804       if(fFillSSHistograms)
 805       {
 806         for(Int_t i = 0; i < 6; i++)
 807         {
 808           fhMCELambda0[pidIndex][i]  = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 809                                                 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
 810                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 811           fhMCELambda0[pidIndex][i]->SetYTitle("#lambda_{0}^{2}");
 812           fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
 813           outputContainer->Add(fhMCELambda0[pidIndex][i]) ;
 814
 815           if(GetCalorimeter()=="EMCAL" && !fFillOnlySimpleSSHisto)
 816           {
 817             fhMCEDispEta[pidIndex][i]  = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 818                                                    Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
 819                                                    nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 820             fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
 821             fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
 822             outputContainer->Add(fhMCEDispEta[pidIndex][i]);
 823
 824             fhMCEDispPhi[pidIndex][i]  = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 825                                                    Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
 826                                                    nptbins,ptmin,ptmax, ssbins,ssmin,ssmax);
 827             fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
 828             fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
 829             outputContainer->Add(fhMCEDispPhi[pidIndex][i]);
 830
 831             fhMCESumEtaPhi[pidIndex][i]  = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 832                                                      Form("cluster from %s : %s like, #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
 833                                                      nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax);
 834             fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
 835             fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
 836             outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
 837
 838             fhMCEDispEtaPhiDiff[pidIndex][i]  = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 839                                                           Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
 840                                                           nptbins,ptmin,ptmax,200,-10,10);
 841             fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
 842             fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
 843             outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);
 844
 845             fhMCESphericity[pidIndex][i]  = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
 846                                                       Form("cluster from %s : %s like, (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
 847                                                       nptbins,ptmin,ptmax, 200,-1,1);
 848             fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
 849             fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
 850             outputContainer->Add(fhMCESphericity[pidIndex][i]);
 851           }
 852
 853         }// loop
 854       }
 855     }
 856
 857     //if(IsCaloPIDOn() && pidIndex > 0) continue;
 858
 859     fhNCellsE[pidIndex]  = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
 860                                      Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
 861                                      nptbins,ptmin,ptmax, nbins,nmin,nmax);
 862     fhNCellsE[pidIndex]->SetXTitle("E (GeV)");
 863     fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
 864     outputContainer->Add(fhNCellsE[pidIndex]);
 865
 866     fhNLME[pidIndex]  = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
 867                                      Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
 868                                      nptbins,ptmin,ptmax, 10,0,10);
 869     fhNLME[pidIndex]->SetXTitle("E (GeV)");
 870     fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
 871     outputContainer->Add(fhNLME[pidIndex]);
 872
 873     fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
 874                                  Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
 875                                  ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
 876     fhTimeE[pidIndex]->SetXTitle("E (GeV)");
 877     fhTimeE[pidIndex]->SetYTitle(" t (ns)");
 878     outputContainer->Add(fhTimeE[pidIndex]);
 879
 880     fhMaxCellDiffClusterE[pidIndex]  = new TH2F (Form("h%sMaxCellDiffClusterE",pidParticle[pidIndex].Data()),
 881                                                  Form("%s: energy vs difference of cluster energy - max cell energy / cluster energy, good clusters",pidParticle[pidIndex].Data()),
 882                                                  nptbins,ptmin,ptmax, 500,0,1.);
 883     fhMaxCellDiffClusterE[pidIndex]->SetXTitle("E_{cluster} (GeV) ");
 884     fhMaxCellDiffClusterE[pidIndex]->SetYTitle("(E_{cluster} - E_{cell max})/ E_{cluster}");
 885     outputContainer->Add(fhMaxCellDiffClusterE[pidIndex]);
 886
 887     fhE[pidIndex]  = new TH1F(Form("h%sE",pidParticle[pidIndex].Data()),
 888                               Form("Number of %s over calorimeter vs energy",pidParticle[pidIndex].Data()),
 889                               nptbins,ptmin,ptmax);
 890     fhE[pidIndex]->SetYTitle("N");
 891     fhE[pidIndex]->SetXTitle("E_{#gamma}(GeV)");
 892     outputContainer->Add(fhE[pidIndex]) ;
 893
 894     fhPt[pidIndex]  = new TH1F(Form("h%sPtElectron",pidParticle[pidIndex].Data()),
 895                                Form("Number of %s over calorimeter vs p_{T}",pidParticle[pidIndex].Data()),
 896                                nptbins,ptmin,ptmax);
 897     fhPt[pidIndex]->SetYTitle("N");
 898     fhPt[pidIndex]->SetXTitle("p_{T #gamma}(GeV/c)");
 899     outputContainer->Add(fhPt[pidIndex]) ;
 900
 901     fhPhi[pidIndex]  = new TH2F(Form("h%sPhiElectron",pidParticle[pidIndex].Data()),
 902                                 Form("%s: #phi vs p_{T}",pidParticle[pidIndex].Data()),
 903                                 nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 904     fhPhi[pidIndex]->SetYTitle("#phi (rad)");
 905     fhPhi[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
 906     outputContainer->Add(fhPhi[pidIndex]) ;
 907
 908     fhEta[pidIndex]  = new TH2F(Form("h%sEta",pidParticle[pidIndex].Data()),
 909                                 Form("%s: #eta vs p_{T}",pidParticle[pidIndex].Data()),
 910                                 nptbins,ptmin,ptmax,netabins,etamin,etamax);
 911     fhEta[pidIndex]->SetYTitle("#eta");
 912     fhEta[pidIndex]->SetXTitle("p_{T #gamma} (GeV/c)");
 913     outputContainer->Add(fhEta[pidIndex]) ;
 914
 915     fhEtaPhi[pidIndex]  = new TH2F(Form("h%sEtaPhi",pidParticle[pidIndex].Data()),
 916                                    Form("%s: #eta vs #phi",pidParticle[pidIndex].Data()),
 917                                    netabins,etamin,etamax,nphibins,phimin,phimax);
 918     fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
 919     fhEtaPhi[pidIndex]->SetXTitle("#eta");
 920     outputContainer->Add(fhEtaPhi[pidIndex]) ;
 921     if(GetMinPt() < 0.5)
 922     {
 923       fhEtaPhi05[pidIndex]  = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
 924                                        Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
 925                                        netabins,etamin,etamax,nphibins,phimin,phimax);
 926       fhEtaPhi05[pidIndex]->SetYTitle("#phi (rad)");
 927       fhEtaPhi05[pidIndex]->SetXTitle("#eta");
 928       outputContainer->Add(fhEtaPhi05[pidIndex]) ;
 929     }
 930
 931
 932     if(IsDataMC())
 933     {
 934       for(Int_t i = 0; i < fNOriginHistograms; i++)
 935       {
 936         fhMCE[pidIndex][i]  = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 937                                        Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
 938                                        nptbins,ptmin,ptmax);
 939         fhMCE[pidIndex][i]->SetXTitle("E (GeV)");
 940         outputContainer->Add(fhMCE[pidIndex][i]) ;
 941
 942         fhMCPt[pidIndex][i]  = new TH1F(Form("h%sPt_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 943                                         Form("%s like cluster from %s : p_{T} ",pidParticle[pidIndex].Data(),ptype[i].Data()),
 944                                         nptbins,ptmin,ptmax);
 945         fhMCPt[pidIndex][i]->SetXTitle("p_{T} (GeV/c)");
 946         outputContainer->Add(fhMCPt[pidIndex][i]) ;
 947
 948         fhMCEta[pidIndex][i]  = new TH2F(Form("h%sEta_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 949                                          Form("%s like cluster from %s : #eta ",pidParticle[pidIndex].Data(),ptype[i].Data()),
 950                                          nptbins,ptmin,ptmax,netabins,etamin,etamax);
 951         fhMCEta[pidIndex][i]->SetYTitle("#eta");
 952         fhMCEta[pidIndex][i]->SetXTitle("E (GeV)");
 953         outputContainer->Add(fhMCEta[pidIndex][i]) ;
 954
 955         fhMCPhi[pidIndex][i]  = new TH2F(Form("h%sPhi_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 956                                          Form("%s like cluster from %s : #phi ",pidParticle[pidIndex].Data(),ptype[i].Data()),
 957                                          nptbins,ptmin,ptmax,nphibins,phimin,phimax);
 958         fhMCPhi[pidIndex][i]->SetYTitle("#phi (rad)");
 959         fhMCPhi[pidIndex][i]->SetXTitle("E (GeV)");
 960         outputContainer->Add(fhMCPhi[pidIndex][i]) ;
 961
 962
 963         fhMCDeltaE[pidIndex][i]  = new TH2F (Form("h%sDeltaE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 964                                              Form("%s like MC - Reco E from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
 965                                              nptbins,ptmin,ptmax, 200,-50,50);
 966         fhMCDeltaE[pidIndex][i]->SetXTitle("#Delta E (GeV)");
 967         outputContainer->Add(fhMCDeltaE[pidIndex][i]);
 968
 969         fhMC2E[pidIndex][i]  = new TH2F (Form("h%s2E_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
 970                                          Form("%s like E distribution, reconstructed vs generated from %s",pidParticle[pidIndex].Data(),pname[i].Data()),
 971                                          nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 972         fhMC2E[pidIndex][i]->SetXTitle("E_{rec} (GeV)");
 973         fhMC2E[pidIndex][i]->SetYTitle("E_{gen} (GeV)");
 974         outputContainer->Add(fhMC2E[pidIndex][i]);
 975
 976       }
 977     } // MC
 978
 979   }// pid Index
 980
 981
 982   if(fFillSSHistograms)
 983   {
 984     if(IsDataMC())
 985     {
 986       if(!GetReader()->IsEmbeddedClusterSelectionOn())
 987       {
 988         fhMCElectronELambda0NoOverlap  = new TH2F("hELambda0_MCElectron_NoOverlap",
 989                                                   "cluster from Electron : E vs #lambda_{0}^{2}",
 990                                                   nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 991         fhMCElectronELambda0NoOverlap->SetYTitle("#lambda_{0}^{2}");
 992         fhMCElectronELambda0NoOverlap->SetXTitle("E (GeV)");
 993         outputContainer->Add(fhMCElectronELambda0NoOverlap) ;
 994
 995         fhMCElectronELambda0TwoOverlap  = new TH2F("hELambda0_MCElectron_TwoOverlap",
 996                                                    "cluster from Electron : E vs #lambda_{0}^{2}",
 997                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
 998         fhMCElectronELambda0TwoOverlap->SetYTitle("#lambda_{0}^{2}");
 999         fhMCElectronELambda0TwoOverlap->SetXTitle("E (GeV)");
1000         outputContainer->Add(fhMCElectronELambda0TwoOverlap) ;
1001
1002         fhMCElectronELambda0NOverlap  = new TH2F("hELambda0_MCElectron_NOverlap",
1003                                                  "cluster from Electron : E vs #lambda_{0}^{2}",
1004                                                  nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1005         fhMCElectronELambda0NOverlap->SetYTitle("#lambda_{0}^{2}");
1006         fhMCElectronELambda0NOverlap->SetXTitle("E (GeV)");
1007         outputContainer->Add(fhMCElectronELambda0NOverlap) ;
1008
1009       } //No embedding
1010
1011       //Fill histograms to check shape of embedded clusters
1012       if(GetReader()->IsEmbeddedClusterSelectionOn())
1013       {
1014
1015         fhEmbeddedSignalFractionEnergy  = new TH2F("hEmbeddedSignal_FractionEnergy",
1016                                                    "Energy Fraction of embedded signal versus cluster energy",
1017                                                    nptbins,ptmin,ptmax,100,0.,1.);
1018         fhEmbeddedSignalFractionEnergy->SetYTitle("Fraction");
1019         fhEmbeddedSignalFractionEnergy->SetXTitle("E (GeV)");
1020         outputContainer->Add(fhEmbeddedSignalFractionEnergy) ;
1021
1022         fhEmbedElectronELambda0FullSignal  = new TH2F("hELambda0_EmbedElectron_FullSignal",
1023                                                       "cluster from Electron embedded with more than 90% energy in cluster : E vs #lambda_{0}^{2}",
1024                                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1025         fhEmbedElectronELambda0FullSignal->SetYTitle("#lambda_{0}^{2}");
1026         fhEmbedElectronELambda0FullSignal->SetXTitle("E (GeV)");
1027         outputContainer->Add(fhEmbedElectronELambda0FullSignal) ;
1028
1029         fhEmbedElectronELambda0MostlySignal  = new TH2F("hELambda0_EmbedElectron_MostlySignal",
1030                                                         "cluster from Electron embedded with 50% to 90% energy in cluster : E vs #lambda_{0}^{2}",
1031                                                         nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1032         fhEmbedElectronELambda0MostlySignal->SetYTitle("#lambda_{0}^{2}");
1033         fhEmbedElectronELambda0MostlySignal->SetXTitle("E (GeV)");
1034         outputContainer->Add(fhEmbedElectronELambda0MostlySignal) ;
1035
1036         fhEmbedElectronELambda0MostlyBkg  = new TH2F("hELambda0_EmbedElectron_MostlyBkg",
1037                                                      "cluster from Electron embedded with 10% to 50% energy in cluster : E vs #lambda_{0}^{2}",
1038                                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1039         fhEmbedElectronELambda0MostlyBkg->SetYTitle("#lambda_{0}^{2}");
1040         fhEmbedElectronELambda0MostlyBkg->SetXTitle("E (GeV)");
1041         outputContainer->Add(fhEmbedElectronELambda0MostlyBkg) ;
1042
1043         fhEmbedElectronELambda0FullBkg  = new TH2F("hELambda0_EmbedElectron_FullBkg",
1044                                                    "cluster from Electronm embedded with 0% to 10% energy in cluster : E vs #lambda_{0}^{2}",
1045                                                    nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
1046         fhEmbedElectronELambda0FullBkg->SetYTitle("#lambda_{0}^{2}");
1047         fhEmbedElectronELambda0FullBkg->SetXTitle("E (GeV)");
1048         outputContainer->Add(fhEmbedElectronELambda0FullBkg) ;
1049
1050
1051       }// embedded histograms
1052
1053     }//Histos with MC
1054
1055   }// Fill SS MC histograms
1056
1057   return outputContainer ;
1058
1059 }
1060
1061 //_________________________
1062 void AliAnaElectron::Init()
1063 {
1064
1065   //Init
1066   //Do some checks
1067   if(GetCalorimeter() == "PHOS" && !GetReader()->IsPHOSSwitchedOn() && NewOutputAOD()){
1068     printf("AliAnaElectron::Init() - !!STOP: You want to use PHOS in analysis but it is not read!! \n!!Check the configuration file!!\n");
1069     abort();
1070   }
1071   else  if(GetCalorimeter() == "EMCAL" && !GetReader()->IsEMCALSwitchedOn() && NewOutputAOD()){
1072     printf("AliAnaElectron::Init() - !!STOP: You want to use EMCAL in analysis but it is not read!! \n!!Check the configuration file!!\n");
1073     abort();
1074   }
1075
1076 }
1077
1078 //___________________________________
1079 void AliAnaElectron::InitParameters()
1080 {
1081
1082   //Initialize the parameters of the analysis.
1083   AddToHistogramsName("AnaElectron_");
1084
1085   fMinDist     = 2.;
1086   fMinDist2    = 4.;
1087   fMinDist3    = 5.;
1088
1089   fTimeCutMin  = -1;
1090   fTimeCutMax  = 9999999;
1091   fNCellsCut   = 0;
1092
1093   fdEdxMin     = 76.; // for LHC11a, but for LHC11c pass1 56.
1094   fdEdxMax     = 85.; // for LHC11a, but for LHC11c pass1 64.
1095
1096   fEOverPMin   = 0.8; // for LHC11a, but for LHC11c pass1 0.9
1097   fEOverPMax   = 1.2; // for LHC11a and LHC11c pass1
1098
1099 }
1100
1101 //_________________________________________
1102 void  AliAnaElectron::MakeAnalysisFillAOD()
1103 {
1104   //Do photon analysis and fill aods
1105
1106   //Get the vertex
1107   Double_t v[3] = {0,0,0}; //vertex ;
1108   GetReader()->GetVertex(v);
1109
1110   //Select the Calorimeter of the photon
1111   TObjArray * pl = 0x0;
1112   if(GetCalorimeter() == "PHOS")
1113     pl = GetPHOSClusters();
1114   else if (GetCalorimeter() == "EMCAL")
1115     pl = GetEMCALClusters();
1116
1117   if(!pl)
1118   {
1119     Info("MakeAnalysisFillAOD","TObjArray with %s clusters is NULL!\n",GetCalorimeter().Data());
1120     return;
1121   }
1122
1123   //Init arrays, variables, get number of clusters
1124   TLorentzVector mom, mom2 ;
1125   Int_t nCaloClusters = pl->GetEntriesFast();
1126   //List to be used in conversion analysis, to tag the cluster as candidate for conversion
1127
1128   if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillAOD() - input %s cluster entries %d\n", GetCalorimeter().Data(), nCaloClusters);
1129
1130   //----------------------------------------------------
1131   // Fill AOD with PHOS/EMCAL AliAODPWG4Particle objects
1132   //----------------------------------------------------
1133   // Loop on clusters
1134   for(Int_t icalo = 0; icalo < nCaloClusters; icalo++)
1135   {
1136           AliVCluster * calo =  (AliVCluster*) (pl->At(icalo));
1137     //printf("calo %d, %f\n",icalo,calo->E());
1138
1139     //Get the index where the cluster comes, to retrieve the corresponding vertex
1140     Int_t evtIndex = 0 ;
1141     if (GetMixedEvent())
1142     {
1143       evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
1144       //Get the vertex and check it is not too large in z
1145       if(TMath::Abs(GetVertex(evtIndex)[2])> GetZvertexCut()) continue;
1146     }
1147
1148     //Cluster selection, not charged, with photon id and in fiducial cut
1149     if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
1150     {
1151       calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
1152     else{
1153       Double_t vertex[]={0,0,0};
1154       calo->GetMomentum(mom,vertex) ;
1155     }
1156
1157     //--------------------------------------
1158     // Cluster selection
1159     //--------------------------------------
1160     AliVCaloCells* cells    = 0;
1161     if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
1162     else                        cells = GetPHOSCells();
1163
1164     Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
1165     if(!ClusterSelected(calo,mom,nMaxima)) continue;
1166
1167     //-------------------------------------
1168     // PID selection via dE/dx
1169     //-------------------------------------
1170
1171     AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
1172
1173     if(!track)
1174     {
1175       printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
1176       continue;
1177     }
1178
1179     //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
1180     Float_t dEdx = track->GetTPCsignal();
1181     Float_t eOverp = calo->E()/track->P();
1182
1183     fhdEdxvsE->Fill(calo->E(), dEdx);
1184     fhdEdxvsP->Fill(track->P(),dEdx);
1185
1186     if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1187     {
1188       fhdEdxvsECutEOverP  ->Fill(calo->E(), dEdx);
1189       fhdEdxvsPCutEOverP  ->Fill(track->P(),dEdx);
1190     }
1191
1192     //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
1193     Float_t m02 = calo->GetM02();
1194     if(m02 > 0.1 && m02 < 0.4)
1195     {
1196       fhdEdxvsECutM02  ->Fill(calo->E(), dEdx);
1197       fhdEdxvsPCutM02  ->Fill(track->P(),dEdx);
1198       fhEOverPvsECutM02->Fill(calo->E(),  eOverp);
1199       fhEOverPvsPCutM02->Fill(track->P(), eOverp);
1200     }
1201
1202     Int_t pid  = AliCaloPID::kChargedHadron;
1203
1204     if( dEdx < fdEdxMax && dEdx > fdEdxMin)
1205     {
1206       fhEOverPvsE->Fill(calo->E(),  eOverp);
1207       fhEOverPvsP->Fill(track->P(), eOverp);
1208
1209       if(m02 > 0.1 && m02 < 0.4)
1210       {
1211         fhEOverPvsECutM02CutdEdx->Fill(calo->E(),  eOverp);
1212         fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
1213       }
1214
1215       if( eOverp < fEOverPMax && eOverp > fEOverPMin)
1216       {
1217         pid  = AliCaloPID::kElectron;
1218       } // E/p
1219
1220     }// dE/dx
1221
1222     Int_t pidIndex = 0;// Electron
1223     if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
1224
1225     //--------------------------------------------------------------------------------------
1226     //Play with the MC stack if available
1227     //--------------------------------------------------------------------------------------
1228
1229     //Check origin of the candidates
1230     Int_t tag = -1 ;
1231     if(IsDataMC())
1232     {
1233       tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(),GetCalorimeter());
1234
1235       if(GetDebug() > 0)
1236         printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",tag);
1237
1238       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1239       {
1240         fhMCdEdxvsE  [kmcPhoton]->Fill(calo ->E(), dEdx);
1241         fhMCdEdxvsP  [kmcPhoton]->Fill(track->P(), dEdx);
1242         fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
1243         fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
1244
1245         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1246         {
1247           fhMCdEdxvsE  [kmcConversion]->Fill(calo ->E(), dEdx);
1248           fhMCdEdxvsP  [kmcConversion]->Fill(track->P(), dEdx);
1249           fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
1250           fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
1251         }
1252         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1253                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1254         {
1255           fhMCdEdxvsE  [kmcPi0Decay]->Fill(calo ->E(), dEdx);
1256           fhMCdEdxvsP  [kmcPi0Decay]->Fill(track->P(), dEdx);
1257           fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
1258           fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
1259         }
1260         else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1261                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1262         {
1263           fhMCdEdxvsE  [kmcOtherDecay]->Fill(calo ->E(), dEdx);
1264           fhMCdEdxvsP  [kmcOtherDecay]->Fill(track->P(), dEdx);
1265           fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
1266           fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
1267         }
1268         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [pidIndex][kmcPi0])
1269         {
1270           fhMCdEdxvsE  [kmcPi0]->Fill(calo ->E(), dEdx);
1271           fhMCdEdxvsP  [kmcPi0]->Fill(track->P(), dEdx);
1272           fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
1273           fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
1274         }
1275         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1276         {
1277           fhMCdEdxvsE  [kmcEta]->Fill(calo ->E(), dEdx);
1278           fhMCdEdxvsP  [kmcEta]->Fill(track->P(), dEdx);
1279           fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
1280           fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
1281         }
1282       }
1283       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1284       {
1285         fhMCdEdxvsE  [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
1286         fhMCdEdxvsP  [kmcAntiNeutron]->Fill(track->P(), dEdx);
1287         fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
1288         fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
1289       }
1290       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1291       {
1292         fhMCdEdxvsE  [kmcAntiProton]->Fill(calo ->E(), dEdx);
1293         fhMCdEdxvsP  [kmcAntiProton]->Fill(track->P(), dEdx);
1294         fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
1295         fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
1296       }
1297       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1298       {
1299         fhMCdEdxvsE  [kmcElectron]->Fill(calo ->E(), dEdx);
1300         fhMCdEdxvsP  [kmcElectron]->Fill(track->P(), dEdx);
1301         fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
1302         fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
1303       }
1304       else if( fhMCE[pidIndex][kmcOther])
1305       {
1306         fhMCdEdxvsE  [kmcOther]->Fill(calo ->E(), dEdx);
1307         fhMCdEdxvsP  [kmcOther]->Fill(track->P(), dEdx);
1308         fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
1309         fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
1310       }
1311     }// set MC tag and fill Histograms with MC
1312
1313     //---------------------------------
1314     //Fill some shower shape histograms
1315     //---------------------------------
1316
1317     FillShowerShapeHistograms(calo,tag,pid);
1318
1319     if(pid == AliCaloPID::kElectron)
1320       WeightHistograms(calo);
1321
1322     //-----------------------------------------
1323     //PID Shower Shape selection or bit setting
1324     //-----------------------------------------
1325
1326     // Data, PID check on
1327     if(IsCaloPIDOn())
1328     {
1329       // Get most probable PID, 2 options check bayesian PID weights or redo PID
1330       // By default, redo PID
1331
1332       if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
1333       {
1334         if(fAODParticle == AliCaloPID::kElectron)
1335           continue;
1336
1337         if(fAODParticle == 0 )
1338           pid = AliCaloPID::kChargedHadron ;
1339       }
1340       if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - PDG of identified particle %d\n",pid);
1341     }
1342
1343     if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
1344                               mom.Pt(), pid);
1345
1346     Float_t maxCellFraction = 0;
1347     Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
1348     if ( absID >= 0 )fhMaxCellDiffClusterE[pidIndex]->Fill(mom.E(),maxCellFraction);
1349
1350     fhNCellsE[pidIndex] ->Fill(mom.E(),calo->GetNCells());
1351     fhNLME   [pidIndex] ->Fill(mom.E(),nMaxima);
1352     fhTimeE  [pidIndex] ->Fill(mom.E(),calo->GetTOF()*1.e9);
1353
1354
1355     //----------------------------
1356     //Create AOD for analysis
1357     //----------------------------
1358
1359     //Add AOD with electron/hadron object to aod branch
1360     if ( pid == fAODParticle || fAODParticle == 0 )
1361     {
1362       AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
1363
1364       //...............................................
1365       //Set the indeces of the original caloclusters (MC, ID), and calorimeter
1366       Int_t label = calo->GetLabel();
1367       aodpart.SetLabel(label);
1368       aodpart.SetCaloLabel (calo ->GetID(),-1);
1369       aodpart.SetTrackLabel(track->GetID(),-1);
1370
1371       aodpart.SetDetector(GetCalorimeter());
1372       //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
1373
1374       //...............................................
1375       //Set bad channel distance bit
1376       Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
1377       if     (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
1378       else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ;
1379       else                         aodpart.SetDistToBad(0) ;
1380       //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
1381
1382       // MC tag
1383       aodpart.SetTag(tag);
1384
1385       // PID tag
1386       aodpart.SetIdentifiedParticleType(pid);
1387
1388       AddAODParticle(aodpart);
1389     }
1390
1391   }//loop
1392
1393   if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());
1394
1395 }
1396
1397 //________________________________________________
1398 void  AliAnaElectron::MakeAnalysisFillHistograms()
1399 {
1400   //Fill histograms
1401
1402   //-------------------------------------------------------------------
1403   // Access MC information in stack if requested, check that it exists.
1404   AliStack         * stack       = 0x0;
1405   TParticle        * primary     = 0x0;
1406   TClonesArray     * mcparticles = 0x0;
1407   AliAODMCParticle * aodprimary  = 0x0;
1408
1409   if(IsDataMC())
1410   {
1411     if(GetReader()->ReadStack())
1412     {
1413       stack =  GetMCStack() ;
1414       if(!stack)
1415       {
1416         printf("AliAnaElectron::MakeAnalysisFillHistograms() - Stack not available, is the MC handler called? STOP\n");
1417         abort();
1418       }
1419     }
1420     else if(GetReader()->ReadAODMCParticles())
1421     {
1422       //Get the list of MC particles
1423       mcparticles = GetReader()->GetAODMCParticles();
1424       if(!mcparticles)
1425       {
1426         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
1427         abort();
1428       }
1429     }
1430   }// is data and MC
1431
1432
1433   // Get vertex
1434   Double_t v[3] = {0,0,0}; //vertex ;
1435   GetReader()->GetVertex(v);
1436   //fhVertex->Fill(v[0],v[1],v[2]);
1437   if(TMath::Abs(v[2]) > GetZvertexCut()) return ; // done elsewhere for Single Event analysis, but there for mixed event
1438
1439   //----------------------------------
1440   //Loop on stored AOD photons
1441   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1442   if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1443
1444   for(Int_t iaod = 0; iaod < naod ; iaod++)
1445   {
1446     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1447     Int_t pdg = ph->GetIdentifiedParticleType();
1448
1449     Int_t pidIndex = 0;// Electron
1450     if     (pdg == AliCaloPID::kElectron)      pidIndex = 0;
1451     else if(pdg == AliCaloPID::kChargedHadron) pidIndex = 1;
1452     else                                       continue    ;
1453
1454     if(ph->GetDetector() != GetCalorimeter()) continue;
1455
1456     if(GetDebug() > 2)
1457       printf("AliAnaElectron::MakeAnalysisFillHistograms() - ID Electron: pt %f, phi %f, eta %f\n", ph->Pt(),ph->Phi(),ph->Eta()) ;
1458
1459     //................................
1460     //Fill photon histograms
1461     Float_t ptcluster  = ph->Pt();
1462     Float_t phicluster = ph->Phi();
1463     Float_t etacluster = ph->Eta();
1464     Float_t ecluster   = ph->E();
1465
1466     fhE[pidIndex]   ->Fill(ecluster);
1467     fhPt[pidIndex]  ->Fill(ptcluster);
1468     fhPhi[pidIndex] ->Fill(ptcluster,phicluster);
1469     fhEta[pidIndex] ->Fill(ptcluster,etacluster);
1470     if     (ecluster > 0.5)   fhEtaPhi[pidIndex]  ->Fill(etacluster, phicluster);
1471     else if(GetMinPt() < 0.5) fhEtaPhi05[pidIndex]->Fill(etacluster, phicluster);
1472
1473     //.......................................
1474     //Play with the MC data if available
1475     if(IsDataMC()){
1476
1477       //....................................................................
1478       // Access MC information in stack if requested, check that it exists.
1479       Int_t label =ph->GetLabel();
1480       if(label < 0) {
1481         if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillHistograms() *** bad label ***:  label %d \n", label);
1482         continue;
1483       }
1484
1485       Float_t eprim   = 0;
1486       //Float_t ptprim  = 0;
1487       if(GetReader()->ReadStack())
1488       {
1489         if(label >=  stack->GetNtrack())
1490         {
1491           if(GetDebug() > 2)  printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", label, stack->GetNtrack());
1492           continue ;
1493         }
1494
1495         primary = stack->Particle(label);
1496         if(!primary)
1497         {
1498           printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1499           continue;
1500         }
1501
1502         eprim   = primary->Energy();
1503         //ptprim  = primary->Pt();
1504
1505       }
1506       else if(GetReader()->ReadAODMCParticles())
1507       {
1508         //Check which is the input
1509         if(ph->GetInputFileIndex() == 0)
1510         {
1511           if(!mcparticles) continue;
1512
1513           if(label >=  mcparticles->GetEntriesFast())
1514           {
1515             if(GetDebug() > 2)  printf("AliAnaElectron::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n",
1516                                        label, mcparticles->GetEntriesFast());
1517             continue ;
1518           }
1519           //Get the particle
1520           aodprimary = (AliAODMCParticle*) mcparticles->At(label);
1521
1522         }
1523
1524         if(!aodprimary)
1525         {
1526           printf("AliAnaElectron::MakeAnalysisFillHistograms() *** no primary ***:  label %d \n", label);
1527           continue;
1528         }
1529
1530         eprim   = aodprimary->E();
1531         //ptprim  = aodprimary->Pt();
1532       }
1533
1534       Int_t tag =ph->GetTag();
1535
1536       if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
1537       {
1538         fhMCE  [pidIndex][kmcPhoton] ->Fill(ecluster);
1539         fhMCPt [pidIndex][kmcPhoton] ->Fill(ptcluster);
1540         fhMCPhi[pidIndex][kmcPhoton] ->Fill(ecluster,phicluster);
1541         fhMCEta[pidIndex][kmcPhoton] ->Fill(ecluster,etacluster);
1542
1543         fhMC2E[pidIndex][kmcPhoton]     ->Fill(ecluster, eprim);
1544         fhMCDeltaE[pidIndex][kmcPhoton] ->Fill(ecluster,eprim-ecluster);
1545
1546         if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
1547         {
1548           fhMCE  [pidIndex][kmcConversion] ->Fill(ecluster);
1549           fhMCPt [pidIndex][kmcConversion] ->Fill(ptcluster);
1550           fhMCPhi[pidIndex][kmcConversion] ->Fill(ecluster,phicluster);
1551           fhMCEta[pidIndex][kmcConversion] ->Fill(ecluster,etacluster);
1552
1553           fhMC2E[pidIndex][kmcConversion]     ->Fill(ecluster, eprim);
1554           fhMCDeltaE[pidIndex][kmcConversion] ->Fill(ecluster,eprim-ecluster);
1555
1556         }
1557         else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
1558                 !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
1559         {
1560           fhMCE  [pidIndex][kmcPi0Decay] ->Fill(ecluster);
1561           fhMCPt [pidIndex][kmcPi0Decay] ->Fill(ptcluster);
1562           fhMCPhi[pidIndex][kmcPi0Decay] ->Fill(ecluster,phicluster);
1563           fhMCEta[pidIndex][kmcPi0Decay] ->Fill(ecluster,etacluster);
1564
1565           fhMC2E[pidIndex][kmcPi0Decay]     ->Fill(ecluster, eprim);
1566           fhMCDeltaE[pidIndex][kmcPi0Decay] ->Fill(ecluster,eprim-ecluster);
1567         }
1568         else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
1569                   GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
1570         {
1571           fhMCE  [pidIndex][kmcOtherDecay] ->Fill(ecluster);
1572           fhMCPt [pidIndex][kmcOtherDecay] ->Fill(ptcluster);
1573           fhMCPhi[pidIndex][kmcOtherDecay] ->Fill(ecluster,phicluster);
1574           fhMCEta[pidIndex][kmcOtherDecay] ->Fill(ecluster,etacluster);
1575
1576           fhMC2E[pidIndex][kmcOtherDecay]     ->Fill(ecluster, eprim);
1577           fhMCDeltaE[pidIndex][kmcOtherDecay] ->Fill(ecluster,eprim-ecluster);
1578         }
1579         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [pidIndex][kmcPi0])
1580         {
1581           fhMCE  [pidIndex][kmcPi0] ->Fill(ecluster);
1582           fhMCPt [pidIndex][kmcPi0] ->Fill(ptcluster);
1583           fhMCPhi[pidIndex][kmcPi0] ->Fill(ecluster,phicluster);
1584           fhMCEta[pidIndex][kmcPi0] ->Fill(ecluster,etacluster);
1585
1586           fhMC2E[pidIndex][kmcPi0]     ->Fill(ecluster, eprim);
1587           fhMCDeltaE[pidIndex][kmcPi0] ->Fill(ecluster,eprim-ecluster);
1588
1589         }
1590         else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
1591         {
1592           fhMCE  [pidIndex][kmcEta] ->Fill(ecluster);
1593           fhMCPt [pidIndex][kmcEta] ->Fill(ptcluster);
1594           fhMCPhi[pidIndex][kmcEta] ->Fill(ecluster,phicluster);
1595           fhMCEta[pidIndex][kmcEta] ->Fill(ecluster,etacluster);
1596
1597           fhMC2E[pidIndex][kmcEta]     ->Fill(ecluster, eprim);
1598           fhMCDeltaE[pidIndex][kmcEta] ->Fill(ecluster,eprim-ecluster);
1599
1600         }
1601       }
1602       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
1603       {
1604         fhMCE  [pidIndex][kmcAntiNeutron] ->Fill(ecluster);
1605         fhMCPt [pidIndex][kmcAntiNeutron] ->Fill(ptcluster);
1606         fhMCPhi[pidIndex][kmcAntiNeutron] ->Fill(ecluster,phicluster);
1607         fhMCEta[pidIndex][kmcAntiNeutron] ->Fill(ecluster,etacluster);
1608
1609         fhMC2E[pidIndex][kmcAntiNeutron]     ->Fill(ecluster, eprim);
1610         fhMCDeltaE[pidIndex][kmcAntiNeutron] ->Fill(ecluster,eprim-ecluster);
1611
1612       }
1613       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
1614       {
1615         fhMCE  [pidIndex][kmcAntiProton] ->Fill(ecluster);
1616         fhMCPt [pidIndex][kmcAntiProton] ->Fill(ptcluster);
1617         fhMCPhi[pidIndex][kmcAntiProton] ->Fill(ecluster,phicluster);
1618         fhMCEta[pidIndex][kmcAntiProton] ->Fill(ecluster,etacluster);
1619
1620         fhMC2E[pidIndex][kmcAntiProton]     ->Fill(ecluster, eprim);
1621         fhMCDeltaE[pidIndex][kmcAntiProton] ->Fill(ecluster,eprim-ecluster);
1622
1623       }
1624       else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
1625       {
1626         fhMCE  [pidIndex][kmcElectron] ->Fill(ecluster);
1627         fhMCPt [pidIndex][kmcElectron] ->Fill(ptcluster);
1628         fhMCPhi[pidIndex][kmcElectron] ->Fill(ecluster,phicluster);
1629         fhMCEta[pidIndex][kmcElectron] ->Fill(ecluster,etacluster);
1630
1631         fhMC2E[pidIndex][kmcElectron]     ->Fill(ecluster, eprim);
1632         fhMCDeltaE[pidIndex][kmcElectron] ->Fill(ecluster,eprim-ecluster);
1633
1634       }
1635       else if( fhMCE[pidIndex][kmcOther])
1636       {
1637         fhMCE  [pidIndex][kmcOther] ->Fill(ecluster);
1638         fhMCPt [pidIndex][kmcOther] ->Fill(ptcluster);
1639         fhMCPhi[pidIndex][kmcOther] ->Fill(ecluster,phicluster);
1640         fhMCEta[pidIndex][kmcOther] ->Fill(ecluster,etacluster);
1641
1642         fhMC2E[pidIndex][kmcOther]     ->Fill(ecluster, eprim);
1643         fhMCDeltaE[pidIndex][kmcOther] ->Fill(ecluster,eprim-ecluster);
1644
1645       }
1646
1647     }//Histograms with MC
1648
1649   }// aod loop
1650
1651 }
1652
1653 //____________________________________________________
1654 void AliAnaElectron::Print(const Option_t * opt) const
1655 {
1656   //Print some relevant parameters set for the analysis
1657
1658   if(! opt)
1659     return;
1660
1661   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
1662   AliAnaCaloTrackCorrBaseClass::Print(" ");
1663
1664   printf("Calorimeter            =     %s\n", GetCalorimeter().Data()) ;
1665   printf(" %2.2f < dEdx < %2.2f  \n",fdEdxMin,fdEdxMax) ;
1666   printf(" %2.2f <  E/P < %2.2f  \n",fEOverPMin,fEOverPMax) ;
1667   printf("Min Distance to Bad Channel   = %2.1f\n",fMinDist);
1668   printf("Min Distance to Bad Channel 2 = %2.1f\n",fMinDist2);
1669   printf("Min Distance to Bad Channel 3 = %2.1f\n",fMinDist3);
1670   printf("Time Cut: %3.1f < TOF  < %3.1f\n", fTimeCutMin, fTimeCutMax);
1671   printf("Number of cells in cluster is        > %d \n", fNCellsCut);
1672   printf("    \n") ;
1673
1674 }
1675
1676 //______________________________________________________
1677 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
1678 {
1679   // Calculate weights and fill histograms
1680
1681   if(!fFillWeightHistograms || GetMixedEvent()) return;
1682
1683   AliVCaloCells* cells = 0;
1684   if(GetCalorimeter() == "EMCAL") cells = GetEMCALCells();
1685   else                        cells = GetPHOSCells();
1686
1687   // First recalculate energy in case non linearity was applied
1688   Float_t  energy = 0;
1689   Float_t  ampMax = 0;
1690   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1691
1692     Int_t id       = clus->GetCellsAbsId()[ipos];
1693
1694     //Recalibrate cell energy if needed
1695     Float_t amp = cells->GetCellAmplitude(id);
1696     GetCaloUtils()->RecalibrateCellAmplitude(amp,GetCalorimeter(), id);
1697
1698     energy    += amp;
1699
1700     if(amp> ampMax)
1701       ampMax = amp;
1702
1703   } // energy loop
1704
1705   if(energy <=0 ) {
1706     printf("AliAnaCalorimeterQA::WeightHistograms()- Wrong calculated energy %f\n",energy);
1707     return;
1708   }
1709
1710   //printf("AliAnaElectron::WeightHistograms() - energy %f, ampmax %f, rat %f, lograt %f\n",energy,ampMax,ampMax/energy,TMath::Log(ampMax/energy));
1711   fhEMaxCellClusterRatio   ->Fill(energy,ampMax/energy);
1712   fhEMaxCellClusterLogRatio->Fill(energy,TMath::Log(ampMax/energy));
1713
1714   //Get the ratio and log ratio to all cells in cluster
1715   for (Int_t ipos = 0; ipos < clus->GetNCells(); ipos++) {
1716     Int_t id       = clus->GetCellsAbsId()[ipos];
1717
1718     //Recalibrate cell energy if needed
1719     Float_t amp = cells->GetCellAmplitude(id);
1720     GetCaloUtils()->RecalibrateCellAmplitude(amp, GetCalorimeter(), id);
1721
1722     //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
1723     fhECellClusterRatio   ->Fill(energy,amp/energy);
1724     fhECellClusterLogRatio->Fill(energy,TMath::Log(amp/energy));
1725   }
1726
1727   //Recalculate shower shape for different W0
1728   if(GetCalorimeter()=="EMCAL"){
1729
1730     Float_t l0org = clus->GetM02();
1731     Float_t l1org = clus->GetM20();
1732     Float_t dorg  = clus->GetDispersion();
1733
1734     for(Int_t iw = 0; iw < 14; iw++){
1735
1736       GetCaloUtils()->GetEMCALRecoUtils()->SetW0(1+iw*0.5);
1737       GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), cells, clus);
1738
1739       fhLambda0ForW0[iw]->Fill(energy,clus->GetM02());
1740       //fhLambda1ForW0[iw]->Fill(energy,clus->GetM20());
1741
1742       //printf("\t w %1.1f, l0 %f, l1 %f,\n",3+iw*0.5,clus->GetM02(),clus->GetM20());
1743
1744     } // w0 loop
1745
1746     // Set the original values back
1747     clus->SetM02(l0org);
1748     clus->SetM20(l1org);
1749     clus->SetDispersion(dorg);
1750
1751   }// EMCAL
1752 }
1753
1754