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