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