]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaElectron.cxx
Merge branch 'master' into TPCdev
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaElectron.cxx
index 2ddf6fedc350f9e79d73b640517e893acc0d28f8..6338f1822e91437a50c727a41a992502f3d3dc41 100755 (executable)
@@ -50,16 +50,22 @@ ClassImp(AliAnaElectron)
   
 //________________________________
 AliAnaElectron::AliAnaElectron() : 
-    AliAnaCaloTrackCorrBaseClass(),            fCalorimeter(""), 
+    AliAnaCaloTrackCorrBaseClass(),       fCalorimeter(""), 
     fMinDist(0.),                         fMinDist2(0.),                         fMinDist3(0.), 
     fTimeCutMin(-1),                      fTimeCutMax(999999),         
-    fNCellsCut(0),                        fFillSSHistograms(kFALSE),             
+    fNCellsCut(0),                        fNLMCutMin(-1),                        fNLMCutMax(10),
+    fFillSSHistograms(kFALSE),             fFillOnlySimpleSSHisto(1),
     fFillWeightHistograms(kFALSE),        fNOriginHistograms(8), 
     fdEdxMin(0.),                         fdEdxMax (200.), 
     fEOverPMin(0),                        fEOverPMax (2),
+    fAODParticle(0),
     // Histograms
     fhdEdxvsE(0),                         fhdEdxvsP(0),                 
     fhEOverPvsE(0),                       fhEOverPvsP(0),
+    fhdEdxvsECutM02(0),                   fhdEdxvsPCutM02(0),
+    fhEOverPvsECutM02(0),                 fhEOverPvsPCutM02(0),
+    fhdEdxvsECutEOverP(0),                fhdEdxvsPCutEOverP(0),
+    fhEOverPvsECutM02CutdEdx(0),          fhEOverPvsPCutM02CutdEdx(0),
     // Weight studies
     fhECellClusterRatio(0),               fhECellClusterLogRatio(0),                 
     fhEMaxCellClusterRatio(0),            fhEMaxCellClusterLogRatio(0),    
@@ -74,10 +80,11 @@ AliAnaElectron::AliAnaElectron() :
     fhEmbedElectronELambda0MostlyBkg(0),  fhEmbedElectronELambda0FullBkg(0)        
 {
   //default ctor
-  for(Int_t index = 0; index < 2; index++){
-    
-    fhNCellsE [index] = 0;    
-    fhTimeE   [index] = 0;  
+  for(Int_t index = 0; index < 2; index++)
+  {
+    fhNCellsE [index] = 0;
+    fhNLME    [index] = 0;
+    fhTimeE   [index] = 0;
     fhMaxCellDiffClusterE[index] = 0;
     fhE       [index] = 0;    
     fhPt      [index] = 0;                        
@@ -100,23 +107,47 @@ AliAnaElectron::AliAnaElectron() :
     fhEtaLam0HighE   [index] = 0;             
     fhPhiLam0HighE   [index] = 0; 
     
-    for(Int_t i = 0; i < 10; i++){
+    fhDispEtaE       [index] = 0;                
+    fhDispPhiE       [index] = 0;
+    fhSumEtaE        [index] = 0;                
+    fhSumPhiE        [index] = 0;                
+    fhSumEtaPhiE     [index] = 0;
+    fhDispEtaPhiDiffE[index] = 0;         
+    fhSphericityE    [index] = 0;
+    
+    for(Int_t i = 0; i < 10; i++)
+    {
       fhMCPt     [index][i] = 0;
       fhMCE      [index][i] = 0;
       fhMCPhi    [index][i] = 0;
       fhMCEta    [index][i] = 0;
       fhMCDeltaE [index][i] = 0;                
-      fhMC2E     [index][i] = 0;              
+      fhMC2E     [index][i] = 0;
+      fhMCdEdxvsE       [i] = 0;
+      fhMCdEdxvsP       [i] = 0;
+      fhMCEOverPvsE     [i] = 0;
+      fhMCEOverPvsP     [i] = 0;
     }
     
-    for(Int_t i = 0; i < 6; i++){
-      fhMCELambda0[index][i]= 0;
+    for(Int_t i = 0; i < 6; i++)
+    {
+      fhMCELambda0       [index][i] = 0;
+      fhMCEDispEta       [index][i] = 0;
+      fhMCEDispPhi       [index][i] = 0;
+      fhMCESumEtaPhi     [index][i] = 0;
+      fhMCEDispEtaPhiDiff[index][i] = 0;
+      fhMCESphericity    [index][i] = 0;
     }
     
+    for(Int_t i = 0; i < 5; i++)
+    {
+      fhDispEtaDispPhiEBin[index][i] = 0 ;
+    }
   }
   
   //Weight studies
-  for(Int_t i =0; i < 14; i++){
+  for(Int_t i =0; i < 14; i++)
+  {
     fhLambda0ForW0[i] = 0;
     //fhLambda1ForW0[i] = 0;
   }
@@ -127,7 +158,7 @@ AliAnaElectron::AliAnaElectron() :
 }
 
 //____________________________________________________________________________
-Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom
+Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom, Int_t nMaxima)
 {
   //Select clusters if they pass different cuts
   if(GetDebug() > 2) 
@@ -166,6 +197,11 @@ Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
   }
   else if(GetDebug() > 2)  printf(" Track-matching cut passed \n");
   
+  //...........................................
+  // skip clusters with too many maxima
+  if(nMaxima < fNLMCutMin || nMaxima > fNLMCutMax) return kFALSE ;
+  if(GetDebug() > 2) printf(" \t Cluster %d pass NLM %d of out of range \n",calo->GetID(), nMaxima);
+  
   //.......................................
   //Check Distance to Bad channel, set bit.
   Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
@@ -186,8 +222,8 @@ Bool_t  AliAnaElectron::ClusterSelected(AliVCluster* calo, TLorentzVector mom)
     
 }
 
-//__________________________________________________________________________________________________________
-void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_t mcTag, const Int_t pidTag)
+//______________________________________________________________________________________________
+void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, Int_t mcTag, Int_t pidTag)
 {
   
   //Fill cluster Shower Shape histograms
@@ -205,8 +241,13 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
   Float_t lambda1 = cluster->GetM20();
   Float_t disp    = cluster->GetDispersion()*cluster->GetDispersion();
   
+  Float_t l0   = 0., l1   = 0.;
+  Float_t dispp= 0., dEta = 0., dPhi    = 0.; 
+  Float_t sEta = 0., sPhi = 0., sEtaPhi = 0.;    
+  
   TLorentzVector mom;
-  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
+  if(GetReader()->GetDataType() != AliCaloTrackReader::kMC)
+  {
     cluster->GetMomentum(mom,GetVertex(0)) ;}//Assume that come from vertex in straight line
   else{
     Double_t vertex[]={0,0,0};
@@ -221,25 +262,51 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
   fhLam1E[pidIndex] ->Fill(energy,lambda1);
   fhDispE[pidIndex] ->Fill(energy,disp);
   
-  if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5){
+  if(fCalorimeter == "EMCAL" && GetModuleNumber(cluster) > 5)
+  {
     fhLam0ETRD[pidIndex]->Fill(energy,lambda0);
     fhLam1ETRD[pidIndex]->Fill(energy,lambda1);
     fhDispETRD[pidIndex]->Fill(energy,disp);
   }
   
-  if(energy < 2){
-    fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
-    fhEtaLam0LowE[pidIndex]    ->Fill(eta,   lambda0);
-    fhPhiLam0LowE[pidIndex]    ->Fill(phi,   lambda0);
-  }
-  else {
-    fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
-    fhEtaLam0HighE[pidIndex]   ->Fill(eta,   lambda0);
-    fhPhiLam0HighE[pidIndex]   ->Fill(phi,   lambda0);
-  }
-    
-  if(IsDataMC()){
+  if(!fFillOnlySimpleSSHisto)
+  {
+    if(energy < 2)
+    {
+      fhNCellsLam0LowE[pidIndex] ->Fill(ncells,lambda0);
+      fhEtaLam0LowE[pidIndex]    ->Fill(eta,   lambda0);
+      fhPhiLam0LowE[pidIndex]    ->Fill(phi,   lambda0);
+    }
+    else 
+    {
+      fhNCellsLam0HighE[pidIndex]->Fill(ncells,lambda0);
+      fhEtaLam0HighE[pidIndex]   ->Fill(eta,   lambda0);
+      fhPhiLam0HighE[pidIndex]   ->Fill(phi,   lambda0);
+    }
     
+    if(fCalorimeter == "EMCAL")
+    {
+      GetCaloUtils()->GetEMCALRecoUtils()->RecalculateClusterShowerShapeParameters(GetEMCALGeometry(), GetReader()->GetInputEvent()->GetEMCALCells(), cluster,
+                                                                                   l0, l1, dispp, dEta, dPhi, sEta, sPhi, sEtaPhi);
+      fhDispEtaE        [pidIndex]-> Fill(energy,dEta);
+      fhDispPhiE        [pidIndex]-> Fill(energy,dPhi);
+      fhSumEtaE         [pidIndex]-> Fill(energy,sEta);
+      fhSumPhiE         [pidIndex]-> Fill(energy,sPhi);
+      fhSumEtaPhiE      [pidIndex]-> Fill(energy,sEtaPhi);
+      fhDispEtaPhiDiffE [pidIndex]-> Fill(energy,dPhi-dEta);
+      if(dEta+dPhi>0)fhSphericityE     [pidIndex]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
+      
+      if      (energy < 2 ) fhDispEtaDispPhiEBin[pidIndex][0]->Fill(dEta,dPhi);
+      else if (energy < 4 ) fhDispEtaDispPhiEBin[pidIndex][1]->Fill(dEta,dPhi);
+      else if (energy < 6 ) fhDispEtaDispPhiEBin[pidIndex][2]->Fill(dEta,dPhi);
+      else if (energy < 10) fhDispEtaDispPhiEBin[pidIndex][3]->Fill(dEta,dPhi);
+      else                  fhDispEtaDispPhiEBin[pidIndex][4]->Fill(dEta,dPhi);
+      
+    }
+  }
+  
+  if(IsDataMC())
+  {
     AliVCaloCells* cells = 0;
     if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
     else                        cells = GetPHOSCells();
@@ -265,26 +332,24 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
       
     }  // embedded fraction    
     
-    // Get the fraction of the cluster energy that carries the cell with highest energy
-    Int_t absID             =-1 ;
-    Float_t maxCellFraction = 0.;
-    
-    absID = GetCaloUtils()->GetMaxEnergyCell(cells, cluster,maxCellFraction);
-    
     // Check the origin and fill histograms
-    if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) && 
+    Int_t index = -1;
+
+    if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) &&
        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) &&
        !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0) &&
-       !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)){
-      fhMCELambda0[pidIndex][kmcssPhoton]    ->Fill(energy, lambda0);
-
+       !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta))
+    {
+      index = kmcssPhoton;
             
     }//photon   no conversion
     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron && 
-              !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion))){
-      fhMCELambda0[pidIndex][kmcssElectron]    ->Fill(energy, lambda0);
+              !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)))
+    {
+      index = kmcssElectron;
        
-      if(!GetReader()->IsEmbeddedClusterSelectionOn()){
+      if(!GetReader()->IsEmbeddedClusterSelectionOn())
+      {
         //Check particle overlaps in cluster
         
         //Compare the primary depositing more energy with the rest, if no photon/electron as comon ancestor (conversions), count as other particle
@@ -307,13 +372,13 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
           fhMCElectronELambda0NOverlap   ->Fill(energy, lambda0);
         }
         else {
-          printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d!!", noverlaps);
+          printf("AliAnaElectron::FillShowerShapeHistogram() - n overlaps = %d for ancestor %d!!", noverlaps, ancLabel);
         }
       }//No embedding
       
       //Fill histograms to check shape of embedded clusters
-      if(GetReader()->IsEmbeddedClusterSelectionOn()){
-        
+      if(GetReader()->IsEmbeddedClusterSelectionOn())
+      {
         if     (fraction > 0.9) 
         {
           fhEmbedElectronELambda0FullSignal   ->Fill(energy, lambda0);
@@ -333,20 +398,34 @@ void  AliAnaElectron::FillShowerShapeHistograms(AliVCluster* cluster, const Int_
       } // embedded      
     }//electron
     else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron) && 
-               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) ){
-      fhMCELambda0[pidIndex][kmcssConversion]    ->Fill(energy, lambda0);      
+               GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion) )
+    {
+      index = kmcssConversion;
     }//conversion photon
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  ){
-      fhMCELambda0[pidIndex][kmcssPi0]    ->Fill(energy, lambda0);      
+    else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)  )
+    {
+      index = kmcssPi0;
     }//pi0
-    else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  ){
-      fhMCELambda0[pidIndex][kmcssEta]    ->Fill(energy, lambda0);
-      
+    else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)  )
+    {
+      index = kmcssEta;      
     }//eta    
-    else {
-      fhMCELambda0[pidIndex][kmcssOther]    ->Fill(energy, lambda0);      
+    else 
+    {
+      index = kmcssOther;
     }//other particles 
     
+    fhMCELambda0[pidIndex][index]    ->Fill(energy, lambda0);
+    
+    if(fCalorimeter == "EMCAL" && !fFillOnlySimpleSSHisto)
+    {
+      fhMCEDispEta        [pidIndex][index]-> Fill(energy,dEta);
+      fhMCEDispPhi        [pidIndex][index]-> Fill(energy,dPhi);
+      fhMCESumEtaPhi      [pidIndex][index]-> Fill(energy,sEtaPhi);
+      fhMCEDispEtaPhiDiff [pidIndex][index]-> Fill(energy,dPhi-dEta);
+      if(dEta+dPhi>0)fhMCESphericity     [pidIndex][index]-> Fill(energy,(dPhi-dEta)/(dEta+dPhi));
+    }
+    
   }//MC data
   
 }
@@ -381,7 +460,7 @@ TObjString *  AliAnaElectron::GetAnalysisCuts()
   parList += GetCaloPID()->GetPIDParametersList() ;
   
   //Get parameters set in FiducialCut class (not available yet)
-  //parlist += GetFidCut()->GetFidCutParametersList() 
+  //parlist += GetFidCut()->GetFidCutParametersList()
   
   return new TObjString(parList) ;
 }
@@ -403,7 +482,18 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
   Int_t nPoverEbins = GetHistogramRanges()->GetHistoPOverEBins();       Float_t pOverEmax = GetHistogramRanges()->GetHistoPOverEMax();       Float_t pOverEmin = GetHistogramRanges()->GetHistoPOverEMin();
   Int_t tbins       = GetHistogramRanges()->GetHistoTimeBins() ;        Float_t tmax      = GetHistogramRanges()->GetHistoTimeMax();         Float_t tmin      = GetHistogramRanges()->GetHistoTimeMin();
 
-  fhdEdxvsE  = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax); 
+  
+  // MC labels, titles, for originator particles
+  TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ;
+  TString pnamess[] = { "Photon","Hadron" ,"Pi0"    ,"Eta" ,"Conversion"     ,"Electron"} ;
+  TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
+    "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P"                    } ;
+  
+  TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
+    "Conversion", "Hadron", "AntiNeutron","AntiProton"                        } ;
+
+  
+  fhdEdxvsE  = new TH2F ("hdEdxvsE","matched track <dE/dx> vs cluster E ", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
   fhdEdxvsE->SetXTitle("E (GeV)");
   fhdEdxvsE->SetYTitle("<dE/dx>");
   outputContainer->Add(fhdEdxvsE);  
@@ -424,9 +514,87 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
   outputContainer->Add(fhEOverPvsP);  
   
   
+  fhdEdxvsECutM02  = new TH2F ("hdEdxvsECutM02","matched track <dE/dx> vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+  fhdEdxvsECutM02->SetXTitle("E (GeV)");
+  fhdEdxvsECutM02->SetYTitle("<dE/dx>");
+  outputContainer->Add(fhdEdxvsECutM02);
+  
+  fhdEdxvsPCutM02  = new TH2F ("hdEdxvsPCutM02","matched track <dE/dx> vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+  fhdEdxvsPCutM02->SetXTitle("P (GeV/c)");
+  fhdEdxvsPCutM02->SetYTitle("<dE/dx>");
+  outputContainer->Add(fhdEdxvsPCutM02);
+  
+  fhEOverPvsECutM02  = new TH2F ("hEOverPvsECutM02","matched track E/p vs cluster E, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+  fhEOverPvsECutM02->SetXTitle("E (GeV)");
+  fhEOverPvsECutM02->SetYTitle("E/p");
+  outputContainer->Add(fhEOverPvsECutM02);
+  
+  fhEOverPvsPCutM02  = new TH2F ("hEOverPvsPCutM02","matched track E/p vs track P, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+  fhEOverPvsPCutM02->SetXTitle("P (GeV/c)");
+  fhEOverPvsPCutM02->SetYTitle("E/p");
+  outputContainer->Add(fhEOverPvsPCutM02);
+
+  
+  fhdEdxvsECutEOverP  = new TH2F ("hdEdxvsECutEOverP","matched track <dE/dx> vs cluster E, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+  fhdEdxvsECutEOverP->SetXTitle("E (GeV)");
+  fhdEdxvsECutEOverP->SetYTitle("<dE/dx>");
+  outputContainer->Add(fhdEdxvsECutEOverP);
+  
+  fhdEdxvsPCutEOverP  = new TH2F ("hdEdxvsPCutEOverP","matched track <dE/dx> vs track P, cut on E/p", nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+  fhdEdxvsPCutEOverP->SetXTitle("P (GeV/c)");
+  fhdEdxvsPCutEOverP->SetYTitle("<dE/dx>");
+  outputContainer->Add(fhdEdxvsPCutEOverP);
+  
+  fhEOverPvsECutM02CutdEdx  = new TH2F ("hEOverPvsECutM02CutdEdx","matched track E/p vs cluster E, dEdx cut, mild #lambda_{0}^{2} cut", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+  fhEOverPvsECutM02CutdEdx->SetXTitle("E (GeV)");
+  fhEOverPvsECutM02CutdEdx->SetYTitle("E/p");
+  outputContainer->Add(fhEOverPvsECutM02CutdEdx);
+  
+  fhEOverPvsPCutM02CutdEdx  = new TH2F ("hEOverPvsPCutM02CutdEdx","matched track E/p vs track P, dEdx cut, mild #lambda_{0}^{2} cut ", nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+  fhEOverPvsPCutM02CutdEdx->SetXTitle("P (GeV/c)");
+  fhEOverPvsPCutM02CutdEdx->SetYTitle("E/p");
+  outputContainer->Add(fhEOverPvsPCutM02CutdEdx);
+
+  if(IsDataMC())
+  {
+    for(Int_t i = 0; i < fNOriginHistograms; i++)
+    {
+      fhMCdEdxvsE[i]  = new TH2F(Form("hdEdxvsE_MC%s",pname[i].Data()),
+                                     Form("matched track <dE/dx> vs cluster E from %s : E ",ptype[i].Data()),
+                                     nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+      fhMCdEdxvsE[i]->SetXTitle("E (GeV)");
+      fhMCdEdxvsE[i]->SetYTitle("<dE/dx>");
+      outputContainer->Add(fhMCdEdxvsE[i]) ;
+      
+      fhMCdEdxvsP[i]  = new TH2F(Form("hdEdxvsP_MC%s",pname[i].Data()),
+                                 Form("matched track <dE/dx> vs track P from %s : E ",ptype[i].Data()),
+                                 nptbins,ptmin,ptmax,ndedxbins, dedxmin, dedxmax);
+      fhMCdEdxvsP[i]->SetXTitle("E (GeV)");
+      fhMCdEdxvsP[i]->SetYTitle("<dE/dx>");
+      outputContainer->Add(fhMCdEdxvsP[i]) ;
+
+      
+      fhMCEOverPvsE[i]  = new TH2F(Form("hEOverPvsE_MC%s",pname[i].Data()),
+                                 Form("matched track E/p vs cluster E from %s : E ",ptype[i].Data()),
+                                 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+      fhMCEOverPvsE[i]->SetXTitle("E (GeV)");
+      fhMCEOverPvsE[i]->SetYTitle("<dE/dx>");
+      outputContainer->Add(fhMCEOverPvsE[i]) ;
+      
+      fhMCEOverPvsP[i]  = new TH2F(Form("hEOverPvsP_MC%s",pname[i].Data()),
+                                 Form("matched track E/pvs track P from %s : E ",ptype[i].Data()),
+                                 nptbins,ptmin,ptmax,nPoverEbins,pOverEmin,pOverEmax);
+      fhMCEOverPvsP[i]->SetXTitle("E (GeV)");
+      fhMCEOverPvsP[i]->SetYTitle("<dE/dx>");
+      outputContainer->Add(fhMCEOverPvsP[i]) ;
+      
+    }
+  }
+  
   TString pidParticle[] = {"Electron","ChargedHadron"} ;
   
-  if(fFillWeightHistograms){
+  if(fFillWeightHistograms)
+  {
     
     fhECellClusterRatio  = new TH2F ("hECellClusterRatio"," cell energy / cluster energy vs cluster energy, for selected electrons",
                                      nptbins,ptmin,ptmax, 100,0,1.); 
@@ -452,7 +620,8 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     fhEMaxCellClusterLogRatio->SetYTitle("Log (E_{max cell}/E_{cluster})");
     outputContainer->Add(fhEMaxCellClusterLogRatio);
     
-    for(Int_t iw = 0; iw < 14; iw++){
+    for(Int_t iw = 0; iw < 14; iw++)
+    {
       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),
                                       nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
       fhLambda0ForW0[iw]->SetXTitle("E_{cluster}");
@@ -468,10 +637,11 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     }
   }  
   
-  for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++){
-    
+  for(Int_t pidIndex = 0; pidIndex < 2; pidIndex++)
+  {
     //Shower shape
-    if(fFillSSHistograms){
+    if(fFillSSHistograms)
+    {
       fhLam0E[pidIndex]  = new TH2F (Form("h%sLam0E",pidParticle[pidIndex].Data()),
                                      Form("%s: #lambda_{0}^{2} vs E",pidParticle[pidIndex].Data()), 
                                      nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
@@ -493,7 +663,8 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
       fhDispE[pidIndex]->SetXTitle("E (GeV) ");
       outputContainer->Add(fhDispE[pidIndex]);
       
-      if(fCalorimeter == "EMCAL"){
+      if(fCalorimeter == "EMCAL")
+      {
         fhLam0ETRD[pidIndex]  = new TH2F (Form("h%sLam0ETRD",pidParticle[pidIndex].Data()),
                                           Form("%s: #lambda_{0}^{2} vs E, EMCAL SM covered by TRD",pidParticle[pidIndex].Data()), 
                                           nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
@@ -516,61 +687,123 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
         outputContainer->Add(fhDispETRD[pidIndex]);   
       } 
       
-      fhNCellsLam0LowE[pidIndex]  = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
-                                              Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
-                                              nbins,nmin, nmax, ssbins,ssmin,ssmax); 
-      fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
-      fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      outputContainer->Add(fhNCellsLam0LowE[pidIndex]);  
-      
-      fhNCellsLam0HighE[pidIndex]  = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
-                                               Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()), 
-                                               nbins,nmin, nmax, ssbins,ssmin,ssmax); 
-      fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
-      fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      outputContainer->Add(fhNCellsLam0HighE[pidIndex]);  
-      
-      
-      fhEtaLam0LowE[pidIndex]  = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
-                                           Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()), 
-                                           netabins,etamin,etamax, ssbins,ssmin,ssmax); 
-      fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
-      outputContainer->Add(fhEtaLam0LowE[pidIndex]);  
-      
-      fhPhiLam0LowE[pidIndex]  = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
-                                           Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()), 
-                                           nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
-      fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
-      outputContainer->Add(fhPhiLam0LowE[pidIndex]);  
-      
-      fhEtaLam0HighE[pidIndex]  = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
-                                            Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
-                                            netabins,etamin,etamax, ssbins,ssmin,ssmax); 
-      fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
-      outputContainer->Add(fhEtaLam0HighE[pidIndex]);  
-      
-      fhPhiLam0HighE[pidIndex]  = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
-                                            Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()), 
-                                            nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
-      fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
-      fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
-      outputContainer->Add(fhPhiLam0HighE[pidIndex]);  
-      
-    } // Shower shape
+      if(!fFillOnlySimpleSSHisto)
+      {
         
-    if(IsDataMC()){
-      
-      if(fFillSSHistograms){
+        fhNCellsLam0LowE[pidIndex]  = new TH2F (Form("h%sNCellsLam0LowE",pidParticle[pidIndex].Data()),
+                                                Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()),
+                                                nbins,nmin, nmax, ssbins,ssmin,ssmax); 
+        fhNCellsLam0LowE[pidIndex]->SetXTitle("N_{cells}");
+        fhNCellsLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        outputContainer->Add(fhNCellsLam0LowE[pidIndex]);  
+        
+        fhNCellsLam0HighE[pidIndex]  = new TH2F (Form("h%sNCellsLam0HighE",pidParticle[pidIndex].Data()),
+                                                 Form("%s: N_{cells} in cluster vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()), 
+                                                 nbins,nmin, nmax, ssbins,ssmin,ssmax); 
+        fhNCellsLam0HighE[pidIndex]->SetXTitle("N_{cells}");
+        fhNCellsLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        outputContainer->Add(fhNCellsLam0HighE[pidIndex]);  
+        
+        
+        fhEtaLam0LowE[pidIndex]  = new TH2F (Form("h%sEtaLam0LowE",pidParticle[pidIndex].Data()),
+                                             Form("%s: #eta vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()), 
+                                             netabins,etamin,etamax, ssbins,ssmin,ssmax); 
+        fhEtaLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        fhEtaLam0LowE[pidIndex]->SetXTitle("#eta");
+        outputContainer->Add(fhEtaLam0LowE[pidIndex]);  
+        
+        fhPhiLam0LowE[pidIndex]  = new TH2F (Form("h%sPhiLam0LowE",pidParticle[pidIndex].Data()),
+                                             Form("%s: #phi vs #lambda_{0}^{2}, E < 2 GeV",pidParticle[pidIndex].Data()), 
+                                             nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
+        fhPhiLam0LowE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        fhPhiLam0LowE[pidIndex]->SetXTitle("#phi");
+        outputContainer->Add(fhPhiLam0LowE[pidIndex]);  
         
-        TString ptypess[] = { "#gamma","hadron?","#pi^{0}","#eta","#gamma->e^{#pm}","e^{#pm}"} ; 
+        fhEtaLam0HighE[pidIndex]  = new TH2F (Form("h%sEtaLam0HighE",pidParticle[pidIndex].Data()),
+                                              Form("%s: #eta vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()),
+                                              netabins,etamin,etamax, ssbins,ssmin,ssmax); 
+        fhEtaLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        fhEtaLam0HighE[pidIndex]->SetXTitle("#eta");
+        outputContainer->Add(fhEtaLam0HighE[pidIndex]);  
         
-        TString pnamess[] = { "Photon","Hadron","Pi0","Eta","Conversion","Electron"} ;
+        fhPhiLam0HighE[pidIndex]  = new TH2F (Form("h%sPhiLam0HighE",pidParticle[pidIndex].Data()),
+                                              Form("%s: #phi vs #lambda_{0}^{2}, E > 2 GeV",pidParticle[pidIndex].Data()), 
+                                              nphibins,phimin,phimax, ssbins,ssmin,ssmax); 
+        fhPhiLam0HighE[pidIndex]->SetYTitle("#lambda_{0}^{2}");
+        fhPhiLam0HighE[pidIndex]->SetXTitle("#phi");
+        outputContainer->Add(fhPhiLam0HighE[pidIndex]);  
         
-        for(Int_t i = 0; i < 6; i++){ 
+        if(fCalorimeter == "EMCAL")
+        {
+          fhDispEtaE[pidIndex]  = new TH2F (Form("h%sDispEtaE",pidParticle[pidIndex].Data()),
+                                            Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),  
+                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+          fhDispEtaE[pidIndex]->SetXTitle("E (GeV)");
+          fhDispEtaE[pidIndex]->SetYTitle("#sigma^{2}_{#eta #eta}");
+          outputContainer->Add(fhDispEtaE[pidIndex]);     
+          
+          fhDispPhiE[pidIndex]  = new TH2F (Form("h%sDispPhiE",pidParticle[pidIndex].Data()),
+                                            Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",pidParticle[pidIndex].Data()),  
+                                            nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+          fhDispPhiE[pidIndex]->SetXTitle("E (GeV)");
+          fhDispPhiE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}");
+          outputContainer->Add(fhDispPhiE[pidIndex]);  
+          
+          fhSumEtaE[pidIndex]  = new TH2F (Form("h%sSumEtaE",pidParticle[pidIndex].Data()),
+                                           Form("%s: #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i})^{2} / #Sigma w_{i} - <#eta>^{2} vs E",pidParticle[pidIndex].Data()),  
+                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+          fhSumEtaE[pidIndex]->SetXTitle("E (GeV)");
+          fhSumEtaE[pidIndex]->SetYTitle("#delta^{2}_{#eta #eta}");
+          outputContainer->Add(fhSumEtaE[pidIndex]);     
+          
+          fhSumPhiE[pidIndex]  = new TH2F (Form("h%sSumPhiE",pidParticle[pidIndex].Data()),
+                                           Form("%s: #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i})^{2}/ #Sigma w_{i} - <#phi>^{2} vs E",pidParticle[pidIndex].Data()),  
+                                           nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+          fhSumPhiE[pidIndex]->SetXTitle("E (GeV)");
+          fhSumPhiE[pidIndex]->SetYTitle("#delta^{2}_{#phi #phi}");
+          outputContainer->Add(fhSumPhiE[pidIndex]);  
+          
+          fhSumEtaPhiE[pidIndex]  = new TH2F (Form("h%sSumEtaPhiE",pidParticle[pidIndex].Data()),
+                                              Form("%s: #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",pidParticle[pidIndex].Data()),  
+                                              nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
+          fhSumEtaPhiE[pidIndex]->SetXTitle("E (GeV)");
+          fhSumEtaPhiE[pidIndex]->SetYTitle("#delta^{2}_{#eta #phi}");
+          outputContainer->Add(fhSumEtaPhiE[pidIndex]);
           
+          fhDispEtaPhiDiffE[pidIndex]  = new TH2F (Form("h%sDispEtaPhiDiffE",pidParticle[pidIndex].Data()),
+                                                   Form("%s: #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",pidParticle[pidIndex].Data()), 
+                                                   nptbins,ptmin,ptmax,200, -10,10); 
+          fhDispEtaPhiDiffE[pidIndex]->SetXTitle("E (GeV)");
+          fhDispEtaPhiDiffE[pidIndex]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+          outputContainer->Add(fhDispEtaPhiDiffE[pidIndex]);    
+          
+          fhSphericityE[pidIndex]  = new TH2F (Form("h%sSphericityE",pidParticle[pidIndex].Data()),
+                                               Form("%s: (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",pidParticle[pidIndex].Data()),  
+                                               nptbins,ptmin,ptmax, 200, -1,1); 
+          fhSphericityE[pidIndex]->SetXTitle("E (GeV)");
+          fhSphericityE[pidIndex]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+          outputContainer->Add(fhSphericityE[pidIndex]);
+          
+          Int_t bin[] = {0,2,4,6,10,1000};
+          for(Int_t i = 0; i < 5; i++)
+          {
+            fhDispEtaDispPhiEBin[pidIndex][i] = new TH2F (Form("h%sDispEtaDispPhi_EBin%d",pidParticle[pidIndex].Data(),i),
+                                                          Form("%s: #sigma^{2}_{#phi #phi} vs #sigma^{2}_{#eta #eta} for %d < E < %d GeV",pidParticle[pidIndex].Data(),bin[i],bin[i+1]), 
+                                                          ssbins,ssmin,ssmax , ssbins,ssmin,ssmax); 
+            fhDispEtaDispPhiEBin[pidIndex][i]->SetXTitle("#sigma^{2}_{#eta #eta}");
+            fhDispEtaDispPhiEBin[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+            outputContainer->Add(fhDispEtaDispPhiEBin[pidIndex][i]); 
+          }
+        }
+      }
+    } // Shower shape
+        
+    if(IsDataMC())
+    {
+      if(fFillSSHistograms)
+      {
+        for(Int_t i = 0; i < 6; i++)
+        { 
           fhMCELambda0[pidIndex][i]  = new TH2F(Form("h%sELambda0_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
                                                 Form("%s like cluster from %s : E vs #lambda_{0}^{2}",pidParticle[pidIndex].Data(),ptypess[i].Data()),
                                                 nptbins,ptmin,ptmax,ssbins,ssmin,ssmax); 
@@ -578,11 +811,49 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
           fhMCELambda0[pidIndex][i]->SetXTitle("E (GeV)");
           outputContainer->Add(fhMCELambda0[pidIndex][i]) ; 
           
+          if(fCalorimeter=="EMCAL" && !fFillOnlySimpleSSHisto)
+          {
+            fhMCEDispEta[pidIndex][i]  = new TH2F (Form("h%sEDispEtaE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
+                                                   Form("cluster from %s : %s like, #sigma^{2}_{#eta #eta} = #Sigma w_{i}(#eta_{i} - <#eta>)^{2}/ #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
+                                                   nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+            fhMCEDispEta[pidIndex][i]->SetXTitle("E (GeV)");
+            fhMCEDispEta[pidIndex][i]->SetYTitle("#sigma^{2}_{#eta #eta}");
+            outputContainer->Add(fhMCEDispEta[pidIndex][i]);     
+            
+            fhMCEDispPhi[pidIndex][i]  = new TH2F (Form("h%sEDispPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
+                                                   Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} = #Sigma w_{i}(#phi_{i} - <#phi>)^{2} / #Sigma w_{i} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),
+                                                   nptbins,ptmin,ptmax, ssbins,ssmin,ssmax); 
+            fhMCEDispPhi[pidIndex][i]->SetXTitle("E (GeV)");
+            fhMCEDispPhi[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}");
+            outputContainer->Add(fhMCEDispPhi[pidIndex][i]);  
+            
+            fhMCESumEtaPhi[pidIndex][i]  = new TH2F (Form("h%sESumEtaPhiE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
+                                                     Form("cluster from %s : %s like, #delta^{2}_{#eta #phi} = #Sigma w_{i}(#phi_{i} #eta_{i} ) / #Sigma w_{i} - <#phi><#eta> vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),  
+                                                     nptbins,ptmin,ptmax, 2*ssbins,-ssmax,ssmax); 
+            fhMCESumEtaPhi[pidIndex][i]->SetXTitle("E (GeV)");
+            fhMCESumEtaPhi[pidIndex][i]->SetYTitle("#delta^{2}_{#eta #phi}");
+            outputContainer->Add(fhMCESumEtaPhi[pidIndex][i]);
+            
+            fhMCEDispEtaPhiDiff[pidIndex][i]  = new TH2F (Form("h%sEDispEtaPhiDiffE_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
+                                                          Form("cluster from %s : %s like, #sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta} vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),  
+                                                          nptbins,ptmin,ptmax,200,-10,10); 
+            fhMCEDispEtaPhiDiff[pidIndex][i]->SetXTitle("E (GeV)");
+            fhMCEDispEtaPhiDiff[pidIndex][i]->SetYTitle("#sigma^{2}_{#phi #phi}-#sigma^{2}_{#eta #eta}");
+            outputContainer->Add(fhMCEDispEtaPhiDiff[pidIndex][i]);    
+            
+            fhMCESphericity[pidIndex][i]  = new TH2F (Form("h%sESphericity_MC%s",pidParticle[pidIndex].Data(),pnamess[i].Data()),
+                                                      Form("cluster from %s : %s like, (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi}) vs E",ptypess[i].Data(),pidParticle[pidIndex].Data()),  
+                                                      nptbins,ptmin,ptmax, 200,-1,1); 
+            fhMCESphericity[pidIndex][i]->SetXTitle("E (GeV)");
+            fhMCESphericity[pidIndex][i]->SetYTitle("s = (#sigma^{2}_{#phi #phi} - #sigma^{2}_{#eta #eta}) / (#sigma^{2}_{#eta #eta} + #sigma^{2}_{#phi #phi})");
+            outputContainer->Add(fhMCESphericity[pidIndex][i]);
+          }
+          
         }// loop    
       }   
     }
     
-    if(IsCaloPIDOn() && pidIndex > 0) continue;
+    //if(IsCaloPIDOn() && pidIndex > 0) continue;
     
     fhNCellsE[pidIndex]  = new TH2F (Form("h%sNCellsE",pidParticle[pidIndex].Data()),
                                      Form("N cells in %s cluster vs E ",pidParticle[pidIndex].Data()),
@@ -591,6 +862,13 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     fhNCellsE[pidIndex]->SetYTitle("# of cells in cluster");
     outputContainer->Add(fhNCellsE[pidIndex]);  
     
+    fhNLME[pidIndex]  = new TH2F (Form("h%sNLME",pidParticle[pidIndex].Data()),
+                                     Form("NLM in %s cluster vs E ",pidParticle[pidIndex].Data()),
+                                     nptbins,ptmin,ptmax, 10,0,10);
+    fhNLME[pidIndex]->SetXTitle("E (GeV)");
+    fhNLME[pidIndex]->SetYTitle("# of cells in cluster");
+    outputContainer->Add(fhNLME[pidIndex]);
+    
     fhTimeE[pidIndex] = new TH2F(Form("h%sTimeE",pidParticle[pidIndex].Data()),
                                  Form("Time in %s cluster vs E ",pidParticle[pidIndex].Data())
                                  ,nptbins,ptmin,ptmax, tbins,tmin,tmax);
@@ -639,7 +917,8 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     fhEtaPhi[pidIndex]->SetYTitle("#phi (rad)");
     fhEtaPhi[pidIndex]->SetXTitle("#eta");
     outputContainer->Add(fhEtaPhi[pidIndex]) ;
-    if(GetMinPt() < 0.5){
+    if(GetMinPt() < 0.5)
+    {
       fhEtaPhi05[pidIndex]  = new TH2F(Form("h%sEtaPhi05",pidParticle[pidIndex].Data()),
                                        Form("%s: #eta vs #phi, E > 0.5",pidParticle[pidIndex].Data()),
                                        netabins,etamin,etamax,nphibins,phimin,phimax); 
@@ -649,17 +928,10 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     }
     
     
-    if(IsDataMC()){
-      
-      
-      TString ptype[] = { "#gamma", "#gamma_{#pi decay}","#gamma_{other decay}", "#pi^{0}","#eta",
-        "e^{#pm}","#gamma->e^{#pm}","hadron?","Anti-N","Anti-P"                    } ; 
-      
-      TString pname[] = { "Photon","PhotonPi0Decay","PhotonOtherDecay","Pi0","Eta","Electron",
-        "Conversion", "Hadron", "AntiNeutron","AntiProton"                        } ;
-      
-      for(Int_t i = 0; i < fNOriginHistograms; i++){ 
-                
+    if(IsDataMC())
+    {      
+      for(Int_t i = 0; i < fNOriginHistograms; i++)
+      { 
         fhMCE[pidIndex][i]  = new TH1F(Form("h%sE_MC%s",pidParticle[pidIndex].Data(),pname[i].Data()),
                                        Form("%s like cluster from %s : E ",pidParticle[pidIndex].Data(),ptype[i].Data()),
                                        nptbins,ptmin,ptmax); 
@@ -706,10 +978,10 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
   }// pid Index
   
   
-  if(fFillSSHistograms){
-    
-    if(IsDataMC()){
-      
+  if(fFillSSHistograms)
+  {
+    if(IsDataMC())
+    {
       if(!GetReader()->IsEmbeddedClusterSelectionOn())
       {
         fhMCElectronELambda0NoOverlap  = new TH2F("hELambda0_MCElectron_NoOverlap",
@@ -781,15 +1053,6 @@ TList *  AliAnaElectron::GetCreateOutputObjects()
     
   }// Fill SS MC histograms
   
-  
-  //Store calo PID histograms
-  if(IsCaloPIDOn()){
-    TList * caloPIDHistos = GetCaloPID()->GetCreateOutputObjects() ;
-    for(Int_t i = 0; i < caloPIDHistos->GetEntries(); i++) {
-      outputContainer->Add(caloPIDHistos->At(i)) ;
-    }
-    delete caloPIDHistos;
-  }
   return outputContainer ;
   
 }
@@ -891,92 +1154,188 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     //--------------------------------------
     // Cluster selection
     //--------------------------------------
-    if(!ClusterSelected(calo,mom)) continue;
+    AliVCaloCells* cells    = 0;
+    if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
+    else                        cells = GetPHOSCells();
+    
+    Int_t nMaxima = GetCaloUtils()->GetNumberOfLocalMaxima(calo, cells); // NLM
+    if(!ClusterSelected(calo,mom,nMaxima)) continue;
     
     //----------------------------
     //Create AOD for analysis
     //----------------------------
-    AliAODPWG4Particle aodph = AliAODPWG4Particle(mom);
+    AliAODPWG4Particle aodpart = AliAODPWG4Particle(mom);
     
     //...............................................
     //Set the indeces of the original caloclusters (MC, ID), and calorimeter  
     Int_t label = calo->GetLabel();
-    aodph.SetLabel(label);
-    aodph.SetCaloLabel(calo->GetID(),-1);
-    aodph.SetDetector(fCalorimeter);
+    aodpart.SetLabel(label);
+    aodpart.SetCaloLabel(calo->GetID(),-1);
+    aodpart.SetDetector(fCalorimeter);
     //printf("Index %d, Id %d, iaod %d\n",icalo, calo->GetID(),GetOutputAODBranch()->GetEntriesFast());
     
     //...............................................
     //Set bad channel distance bit
     Double_t distBad=calo->GetDistanceToBadChannel() ; //Distance to bad channel
-    if     (distBad > fMinDist3) aodph.SetDistToBad(2) ;
-    else if(distBad > fMinDist2) aodph.SetDistToBad(1) ; 
-    else                         aodph.SetDistToBad(0) ;
-    //printf("DistBad %f Bit %d\n",distBad, aodph.DistToBad());
-    
-    //--------------------------------------------------------------------------------------
-    //Play with the MC stack if available
-    //--------------------------------------------------------------------------------------
-    
-    //Check origin of the candidates
-    if(IsDataMC()){
-      aodph.SetTag(GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader(), aodph.GetInputFileIndex()));
-      
-      if(GetDebug() > 0)
-        printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodph.GetTag());
-    }//Work with stack also   
-    
+    if     (distBad > fMinDist3) aodpart.SetDistToBad(2) ;
+    else if(distBad > fMinDist2) aodpart.SetDistToBad(1) ; 
+    else                         aodpart.SetDistToBad(0) ;
+    //printf("DistBad %f Bit %d\n",distBad, aodpart.DistToBad());
     
     //-------------------------------------
     //PID selection via dEdx
     //-------------------------------------
     
-    AliVTrack *track = 0;
-    if(!strcmp("AliESDCaloCluster",Form("%s",calo->ClassName()))){
-      Int_t iESDtrack = calo->GetTrackMatchedIndex();
-      if(iESDtrack<0) printf("AliAnaElectron::MakeAnalysisFillAOD() - Wrong track index\n");
-      AliVEvent * event = GetReader()->GetInputEvent();
-      track = dynamic_cast<AliVTrack*> (event->GetTrack(iESDtrack));
-    }
-    else {
-      track = dynamic_cast<AliVTrack*>(calo->GetTrackMatched(0));
-    }
-    
-    if(!track) {
+    AliVTrack *track = GetCaloUtils()->GetMatchedTrack(calo, GetReader()->GetInputEvent());
+
+    if(!track)
+    {
       printf("AliAnaElectron::MakeAnalysisFillAOD() - Null track");
       continue;
     }
     
+    //printf("track dedx %f, p %f, cluster E %f\n",track->GetTPCsignal(),track->P(),calo->E());
     Float_t dEdx = track->GetTPCsignal();
+    Float_t eOverp = calo->E()/track->P();
+    
     fhdEdxvsE->Fill(calo->E(), dEdx);
     fhdEdxvsP->Fill(track->P(),dEdx);
     
+    if( eOverp < fEOverPMax && eOverp > fEOverPMin)
+    {
+      fhdEdxvsECutEOverP  ->Fill(calo->E(), dEdx);
+      fhdEdxvsPCutEOverP  ->Fill(track->P(),dEdx);
+    }
+    
+    //Apply a mild cut on the cluster SS and check the value of dEdX and EOverP
+    Float_t m02 = calo->GetM02();
+    if(m02 > 0.1 && m02 < 0.4)
+    {
+      fhdEdxvsECutM02  ->Fill(calo->E(), dEdx);
+      fhdEdxvsPCutM02  ->Fill(track->P(),dEdx);
+      fhEOverPvsECutM02->Fill(calo->E(),  eOverp);
+      fhEOverPvsPCutM02->Fill(track->P(), eOverp);
+    }
+    
     Int_t pid  = AliCaloPID::kChargedHadron;
     
-    if( dEdx < fdEdxMax && dEdx > fdEdxMin) {
-      
-      Float_t eOverp = calo->E()/track->P();
+    if( dEdx < fdEdxMax && dEdx > fdEdxMin)
+    {
       fhEOverPvsE->Fill(calo->E(),  eOverp);
       fhEOverPvsP->Fill(track->P(), eOverp);
       
-      if( eOverp < fEOverPMax && eOverp > fEOverPMin) {
-        
+      if(m02 > 0.1 && m02 < 0.4)
+      {
+        fhEOverPvsECutM02CutdEdx->Fill(calo->E(),  eOverp);
+        fhEOverPvsPCutM02CutdEdx->Fill(track->P(), eOverp);
+      }
+      
+      if( eOverp < fEOverPMax && eOverp > fEOverPMin)
+      {
         pid  = AliCaloPID::kElectron;
       } //E/p
       
     }// dEdx
     
-    aodph.SetIdentifiedParticleType(pid);    
+    aodpart.SetIdentifiedParticleType(pid);    
     
     Int_t pidIndex = 0;// Electron
     if(pid == AliCaloPID::kChargedHadron) pidIndex = 1;
+  
+    //--------------------------------------------------------------------------------------
+    //Play with the MC stack if available
+    //--------------------------------------------------------------------------------------
+    
+    //Check origin of the candidates
+    if(IsDataMC())
+    {
+      Int_t tag = GetMCAnalysisUtils()->CheckOrigin(calo->GetLabels(),calo->GetNLabels(),GetReader());
+      aodpart.SetTag(tag);
+      
+      if(GetDebug() > 0)
+        printf("AliAnaElectron::MakeAnalysisFillAOD() - Origin of candidate, bit map %d\n",aodpart.GetTag());
+         
+      if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton) && fhMCE[pidIndex][kmcPhoton])
+      {
+        fhMCdEdxvsE  [kmcPhoton]->Fill(calo ->E(), dEdx);
+        fhMCdEdxvsP  [kmcPhoton]->Fill(track->P(), dEdx);
+        fhMCEOverPvsE[kmcPhoton]->Fill(calo ->E(), eOverp);
+        fhMCEOverPvsP[kmcPhoton]->Fill(track->P(), eOverp);
+        
+        if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCConversion) && fhMCE[pidIndex][kmcConversion])
+        {
+          fhMCdEdxvsE  [kmcConversion]->Fill(calo ->E(), dEdx);
+          fhMCdEdxvsP  [kmcConversion]->Fill(track->P(), dEdx);
+          fhMCEOverPvsE[kmcConversion]->Fill(calo ->E(), eOverp);
+          fhMCEOverPvsP[kmcConversion]->Fill(track->P(), eOverp);
+        }
+        else if( GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0Decay) &&
+                !GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE[pidIndex][kmcPi0Decay])
+        {
+          fhMCdEdxvsE  [kmcPi0Decay]->Fill(calo ->E(), dEdx);
+          fhMCdEdxvsP  [kmcPi0Decay]->Fill(track->P(), dEdx);
+          fhMCEOverPvsE[kmcPi0Decay]->Fill(calo ->E(), eOverp);
+          fhMCEOverPvsP[kmcPi0Decay]->Fill(track->P(), eOverp);
+        }
+        else if( (GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEtaDecay) ||
+                  GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCOtherDecay) ) && fhMCE[pidIndex][kmcOtherDecay])
+        {
+          fhMCdEdxvsE  [kmcOtherDecay]->Fill(calo ->E(), dEdx);
+          fhMCdEdxvsP  [kmcOtherDecay]->Fill(track->P(), dEdx);
+          fhMCEOverPvsE[kmcOtherDecay]->Fill(calo ->E(), eOverp);
+          fhMCEOverPvsP[kmcOtherDecay]->Fill(track->P(), eOverp);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPi0) && fhMCE  [pidIndex][kmcPi0])
+        {
+          fhMCdEdxvsE  [kmcPi0]->Fill(calo ->E(), dEdx);
+          fhMCdEdxvsP  [kmcPi0]->Fill(track->P(), dEdx);
+          fhMCEOverPvsE[kmcPi0]->Fill(calo ->E(), eOverp);
+          fhMCEOverPvsP[kmcPi0]->Fill(track->P(), eOverp);
+        }
+        else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCEta) && fhMCE[pidIndex][kmcEta])
+        {
+          fhMCdEdxvsE  [kmcEta]->Fill(calo ->E(), dEdx);
+          fhMCdEdxvsP  [kmcEta]->Fill(track->P(), dEdx);
+          fhMCEOverPvsE[kmcEta]->Fill(calo ->E(), eOverp);
+          fhMCEOverPvsP[kmcEta]->Fill(track->P(), eOverp);
+        }
+      }
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiNeutron) && fhMCE[pidIndex][kmcAntiNeutron])
+      {
+        fhMCdEdxvsE  [kmcAntiNeutron]->Fill(calo ->E(), dEdx);
+        fhMCdEdxvsP  [kmcAntiNeutron]->Fill(track->P(), dEdx);
+        fhMCEOverPvsE[kmcAntiNeutron]->Fill(calo ->E(), eOverp);
+        fhMCEOverPvsP[kmcAntiNeutron]->Fill(track->P(), eOverp);
+      }
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCAntiProton) && fhMCE[pidIndex][kmcAntiProton])
+      {
+        fhMCdEdxvsE  [kmcAntiProton]->Fill(calo ->E(), dEdx);
+        fhMCdEdxvsP  [kmcAntiProton]->Fill(track->P(), dEdx);
+        fhMCEOverPvsE[kmcAntiProton]->Fill(calo ->E(), eOverp);
+        fhMCEOverPvsP[kmcAntiProton]->Fill(track->P(), eOverp);
+      }
+      else if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCElectron) && fhMCE[pidIndex][kmcElectron])
+      {
+        fhMCdEdxvsE  [kmcElectron]->Fill(calo ->E(), dEdx);
+        fhMCdEdxvsP  [kmcElectron]->Fill(track->P(), dEdx);
+        fhMCEOverPvsE[kmcElectron]->Fill(calo ->E(), eOverp);
+        fhMCEOverPvsP[kmcElectron]->Fill(track->P(), eOverp);
+      }
+      else if( fhMCE[pidIndex][kmcOther])
+      {
+        fhMCdEdxvsE  [kmcOther]->Fill(calo ->E(), dEdx);
+        fhMCdEdxvsP  [kmcOther]->Fill(track->P(), dEdx);
+        fhMCEOverPvsE[kmcOther]->Fill(calo ->E(), eOverp);
+        fhMCEOverPvsP[kmcOther]->Fill(track->P(), eOverp);
+      }
+    }// set MC tag and fill Histograms with MC
     
     //---------------------------------
     //Fill some shower shape histograms
     //---------------------------------
 
-    FillShowerShapeHistograms(calo,aodph.GetTag(),pid);
-    
+    FillShowerShapeHistograms(calo,aodpart.GetTag(),pid);
+  
     if(pid == AliCaloPID::kElectron)
       WeightHistograms(calo);
     
@@ -985,35 +1344,37 @@ void  AliAnaElectron::MakeAnalysisFillAOD()
     //-----------------------------------------
     
     // Data, PID check on
-    if(IsCaloPIDOn()){
+    if(IsCaloPIDOn())
+    {
       // Get most probable PID, 2 options check bayesian PID weights or redo PID
       // By default, redo PID
     
-      if(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo)!=AliCaloPID::kPhoton) continue;
-      
-      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
+      if(GetCaloPID()->GetIdentifiedParticleType(calo)!=AliCaloPID::kPhoton)
+      {
+        continue;
+      }
+      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodpart.GetIdentifiedParticleType());
       
     }
         
-    if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
+    if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",
+                              aodpart.Pt(), aodpart.GetIdentifiedParticleType());
     
     //FIXME, this to MakeAnalysisFillHistograms ...
-    Int_t absID             = 0; 
     Float_t maxCellFraction = 0;
-    AliVCaloCells* cells    = 0;
-    
-    if(fCalorimeter == "EMCAL") cells = GetEMCALCells();
-    else                        cells = GetPHOSCells();
-    
-    absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
     
-    fhMaxCellDiffClusterE[pidIndex]->Fill(aodph.E(),maxCellFraction);
-    fhNCellsE[pidIndex]            ->Fill(aodph.E(),calo->GetNCells());
-    fhTimeE[pidIndex]              ->Fill(aodph.E(),calo->GetTOF()*1.e9);
+    Int_t absID = GetCaloUtils()->GetMaxEnergyCell(cells, calo,maxCellFraction);
+    if(absID>=0)fhMaxCellDiffClusterE[pidIndex]->Fill(aodpart.E(),maxCellFraction);
+    fhNCellsE[pidIndex]            ->Fill(aodpart.E(),calo->GetNCells());
+    fhNLME[pidIndex]               ->Fill(aodpart.E(),nMaxima);
+    fhTimeE[pidIndex]              ->Fill(aodpart.E(),calo->GetTOF()*1.e9);
+
+    //Add AOD with electron/hadron object to aod branch
+    if ( pid == fAODParticle || fAODParticle == 0 ) 
+    {
+      AddAODParticle(aodpart);
+    }
 
-    //Add AOD with photon object to aod branch
-    AddAODParticle(aodph);
-    
   }//loop
   
   if(GetDebug() > 1) printf("AliAnaElectron::MakeAnalysisFillAOD()  End fill AODs, with %d entries \n",GetOutputAODBranch()->GetEntriesFast());  
@@ -1045,7 +1406,7 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
     else if(GetReader()->ReadAODMCParticles()){
       
       //Get the list of MC particles
-      mcparticles = GetReader()->GetAODMCParticles(0);
+      mcparticles = GetReader()->GetAODMCParticles();
       if(!mcparticles && GetDebug() > 0)       {
         printf("AliAnaElectron::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
       }        
@@ -1064,7 +1425,8 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
   if(GetDebug() > 0) printf("AliAnaElectron::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
   
-  for(Int_t iaod = 0; iaod < naod ; iaod++){
+  for(Int_t iaod = 0; iaod < naod ; iaod++)
+  {
     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
     Int_t pdg = ph->GetIdentifiedParticleType();
 
@@ -1105,7 +1467,7 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
       }
       
       Float_t eprim   = 0;
-      Float_t ptprim  = 0;
+      //Float_t ptprim  = 0;
       if(GetReader()->ReadStack()){
         
         if(label >=  stack->GetNtrack()) {
@@ -1120,7 +1482,7 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
         }
         
         eprim   = primary->Energy();
-        ptprim  = primary->Pt();               
+        //ptprim  = primary->Pt();
         
       }
       else if(GetReader()->ReadAODMCParticles()){
@@ -1143,7 +1505,7 @@ void  AliAnaElectron::MakeAnalysisFillHistograms()
         }
         
         eprim   = aodprimary->E();
-        ptprim  = aodprimary->Pt();
+        //ptprim  = aodprimary->Pt();
         
       }
       
@@ -1288,24 +1650,6 @@ void AliAnaElectron::Print(const Option_t * opt) const
        
 } 
 
-//__________________________________________________________________________
-void AliAnaElectron::RecalibrateCellAmplitude(Float_t & amp, const Int_t id)
-{
-  //Recaculate cell energy if recalibration factor
-  
-  Int_t icol     = -1; Int_t irow     = -1; Int_t iRCU     = -1;
-  Int_t nModule  = GetModuleNumberCellIndexes(id,fCalorimeter, icol, irow, iRCU);
-  
-  if (GetCaloUtils()->IsRecalibrationOn()) {
-    if(fCalorimeter == "PHOS") {
-      amp *= GetCaloUtils()->GetPHOSChannelRecalibrationFactor(nModule,icol,irow);
-    }
-    else                                  {
-      amp *= GetCaloUtils()->GetEMCALChannelRecalibrationFactor(nModule,icol,irow);
-    }
-  }
-}
-
 //______________________________________________________
 void AliAnaElectron::WeightHistograms(AliVCluster *clus)
 {
@@ -1326,7 +1670,7 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
     
     //Recalibrate cell energy if needed
     Float_t amp = cells->GetCellAmplitude(id);
-    RecalibrateCellAmplitude(amp,id);
+    GetCaloUtils()->RecalibrateCellAmplitude(amp,fCalorimeter, id);
     
     energy    += amp;
     
@@ -1350,7 +1694,7 @@ void AliAnaElectron::WeightHistograms(AliVCluster *clus)
     
     //Recalibrate cell energy if needed
     Float_t amp = cells->GetCellAmplitude(id);
-    RecalibrateCellAmplitude(amp,id);
+    GetCaloUtils()->RecalibrateCellAmplitude(amp, fCalorimeter, id);
 
     //printf("energy %f, amp %f, rat %f, lograt %f\n",energy,amp,amp/energy,TMath::Log(amp/energy));
     fhECellClusterRatio   ->Fill(energy,amp/energy);