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