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