]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGGA/CaloTrackCorrelations/AliAnaParticleIsolation.cxx
fix wrong initialization of int with a float
[u/mrichter/AliRoot.git] / PWGGA / CaloTrackCorrelations / AliAnaParticleIsolation.cxx
index 1d5e01cb8f94942c32247b5a9eb2e7e75cd079ec..ef5a02e7351c94ee8a0be8792c04286202969dc5 100755 (executable)
@@ -25,7 +25,6 @@
 //-Yaxian Mao (add the possibility for different IC method with different pt range, 01/10/2010)
 //////////////////////////////////////////////////////////////////////////////
 
-
 // --- ROOT system ---
 #include <TClonesArray.h>
 #include <TList.h>
@@ -36,8 +35,6 @@
 #include "TParticle.h"
 #include "TDatabasePDG.h"
 
-
-
 // --- Analysis system ---
 #include "AliAnaParticleIsolation.h"
 #include "AliCaloTrackReader.h"
@@ -68,8 +65,10 @@ fFillPileUpHistograms(0),
 fFillTMHisto(0),                  fFillSSHisto(1),
 fFillUEBandSubtractHistograms(1), fFillCellHistograms(0),
 fFillHighMultHistograms(0),       fFillTaggedDecayHistograms(0),
+fNDecayBits(0),                   fDecayBits(),
 fFillNLMHistograms(0),
 fLeadingOnly(0),                  fCheckLeadingWithNeutralClusters(0),
+fFillBackgroundBinHistograms(0),  fNBkgBin(0),
 // Several IC
 fNCones(0),                       fNPtThresFrac(0),
 fConeSizes(),                     fPtThresholds(),
@@ -81,8 +80,6 @@ fhPtNLocMaxIso(0),
 fhPhiIso(0),                      fhEtaIso(0),                              fhEtaPhiIso(0),
 fhEtaPhiNoIso(0),
 fhENoIso(0),                      fhPtNoIso(0),                             fhPtNLocMaxNoIso(0),
-fhPtDecayIso(0),                  fhPtDecayNoIso(0),
-fhEtaPhiDecayIso(0),              fhEtaPhiDecayNoIso(0),
 fhPtInCone(0),
 fhPtClusterInCone(0),             fhPtCellInCone(0),                        fhPtTrackInCone(0),
 fhPtTrackInConeOtherBC(0),        fhPtTrackInConeOtherBCPileUpSPD(0),
@@ -95,6 +92,8 @@ fhEtaPhiInConeTrack(0),           fhEtaPhiTrack(0),
 fhEtaBandCluster(0),              fhPhiBandCluster(0),
 fhEtaBandTrack(0),                fhPhiBandTrack(0),
 fhEtaBandCell(0),                 fhPhiBandCell(0),
+fhConePtLead(0),                  fhConePtLeadCluster(0),                   fhConePtLeadTrack(0),
+fhConePtLeadClustervsTrack(0),    fhConePtLeadClusterTrackFrac(0),
 fhConeSumPt(0),                   fhConeSumPtCellTrack(0),
 fhConeSumPtCell(0),               fhConeSumPtCluster(0),                    fhConeSumPtTrack(0),
 fhConeSumPtEtaBandUECluster(0),             fhConeSumPtPhiBandUECluster(0),
@@ -120,7 +119,7 @@ fhFractionClusterOutConeEta(0),             fhFractionClusterOutConeEtaTrigEtaPh
 fhFractionClusterOutConePhi(0),             fhFractionClusterOutConePhiTrigEtaPhi(0),
 fhFractionCellOutConeEta(0),                fhFractionCellOutConeEtaTrigEtaPhi(0),
 fhFractionCellOutConePhi(0),                fhFractionCellOutConePhiTrigEtaPhi(0),
-fhConeSumPtClustervsTrack(0),
+fhConeSumPtClustervsTrack(0),               fhConeSumPtClusterTrackFrac(0),
 fhConeSumPtEtaUESubClustervsTrack(0),       fhConeSumPtPhiUESubClustervsTrack(0),
 fhConeSumPtCellvsTrack(0),
 fhConeSumPtEtaUESubCellvsTrack(0),          fhConeSumPtPhiUESubCellvsTrack(0),
@@ -136,7 +135,8 @@ fhConeSumPtSubvsConeSumPtTotPhiCell(0),     fhConeSumPtSubNormvsConeSumPtTotPhiC
 fhConeSumPtSubvsConeSumPtTotEtaCell(0),     fhConeSumPtSubNormvsConeSumPtTotEtaCell(0),
 fhConeSumPtVSUETracksEtaBand(0),            fhConeSumPtVSUETracksPhiBand(0),
 fhConeSumPtVSUEClusterEtaBand(0),           fhConeSumPtVSUEClusterPhiBand(0),
-
+fhPtLeadConeBinLambda0(0),                  fhSumPtConeBinLambda0(0),
+fhPtLeadConeBinLambda0MC(0),                fhSumPtConeBinLambda0MC(0),
 // Number of local maxima in cluster
 fhNLocMax(),
 fhELambda0LocMax1(),              fhELambda1LocMax1(),
@@ -195,6 +195,21 @@ fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
     }
   }
   
+  for(Int_t ibit =0; ibit< 4; ibit++)
+  {
+    fhPtDecayIso       [ibit] = 0;
+    fhPtLambda0Decay[0][ibit] = 0;
+    fhPtLambda0Decay[1][ibit] = 0;
+    fhPtDecayNoIso     [ibit] = 0;
+    fhEtaPhiDecayIso   [ibit] = 0;
+    fhEtaPhiDecayNoIso [ibit] = 0;
+    for(Int_t imc = 0; imc < 9; imc++)
+    {
+      fhPtDecayIsoMC  [ibit][imc]    = 0;
+      fhPtDecayNoIsoMC[ibit][imc]    = 0;
+    }
+  }
+  
   for(Int_t i = 0; i < 5 ; i++)
   {
     fPtFractions    [i] = 0 ;
@@ -215,9 +230,6 @@ fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
     fhEtaIsoMC   [imc]    = 0;
     fhPtLambda0MC[imc][0] = 0;
     fhPtLambda0MC[imc][1] = 0;
-    
-    fhPtDecayIsoMC  [imc]    = 0;
-    fhPtDecayNoIsoMC[imc]    = 0;
   }
   
   for(Int_t i = 0; i < 2 ; i++)
@@ -253,7 +265,6 @@ fhTimePileUpMainVertexZDistance(0), fhTimePileUpMainVertexZDiamond(0)
     fhENoIsoPileUp  [i] = 0 ;
     fhPtNoIsoPileUp [i] = 0 ;
   }
-  
 }
 
 //_______________________________________________________________________________________________
@@ -842,11 +853,13 @@ void AliAnaParticleIsolation::CalculateNormalizeUEBandPerUnitArea(AliAODPWG4Part
 }
 
 
-//__________________________________________________________________________________________________
+//______________________________________________________________________________________________________________
 void AliAnaParticleIsolation::CalculateCaloSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
-                                                        Float_t & coneptsumCluster)
+                                                        Float_t & coneptsumCluster, Float_t & coneptLeadCluster)
 {
   // Get the cluster pT or sum of pT in isolation cone
+  coneptLeadCluster = 0;
+  coneptsumCluster  = 0;
   
   if( GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kOnlyCharged ) return ;
   
@@ -884,9 +897,11 @@ void AliAnaParticleIsolation::CalculateCaloSignalInCone(AliAODPWG4ParticleCorrel
     if(fFillHighMultHistograms) fhPtInConeCent->Fill(GetEventCentrality(),mom.Pt());
     
     coneptsumCluster+=mom.Pt();
+    if(mom.Pt() > coneptLeadCluster) coneptLeadCluster = mom.Pt();
   }
   
-  fhConeSumPtCluster   ->Fill(ptTrig,     coneptsumCluster);
+  fhConeSumPtCluster ->Fill(ptTrig, coneptsumCluster );
+  fhConePtLeadCluster->Fill(ptTrig, coneptLeadCluster);
 }
 
 //______________________________________________________________________________________________________
@@ -971,9 +986,9 @@ void AliAnaParticleIsolation::CalculateCaloCellSignalInCone(AliAODPWG4ParticleCo
   
 }
 
-//___________________________________________________________________________________________________
+//___________________________________________________________________________________________________________
 void AliAnaParticleIsolation::CalculateTrackSignalInCone(AliAODPWG4ParticleCorrelation * aodParticle,
-                                                         Float_t & coneptsumTrack)
+                                                         Float_t & coneptsumTrack, Float_t & coneptLeadTrack)
 {
   // Get the track pT or sum of pT in isolation cone
   
@@ -1021,10 +1036,12 @@ void AliAnaParticleIsolation::CalculateTrackSignalInCone(AliAODPWG4ParticleCorre
     if(fFillHighMultHistograms) fhPtInConeCent->Fill(GetEventCentrality(),pTtrack);
     
     coneptsumTrack+=pTtrack;
+    if(pTtrack > coneptLeadTrack) coneptLeadTrack = pTtrack;
   }
-  
-  fhConeSumPtTrack->Fill(ptTrig, coneptsumTrack);
-  
+
+  fhConeSumPtTrack ->Fill(ptTrig, coneptsumTrack );
+  fhConePtLeadTrack->Fill(ptTrig, coneptLeadTrack);
+
 }
 
 //_________________________________________________________________
@@ -1122,10 +1139,11 @@ void AliAnaParticleIsolation::FillPileUpHistograms(Int_t clusterID)
 
 //_____________________________________________________________________________________________________________________
 void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(AliAODPWG4ParticleCorrelation  *pCandidate,
+                                                                            Float_t coneptsum, Float_t coneleadpt,
                                                                             Int_t mcIndex)
 {
   // Fill Track matching and Shower Shape control histograms
-  if(!fFillTMHisto &&  !fFillSSHisto) return;
+  if(!fFillTMHisto && !fFillSSHisto && !fFillBackgroundBinHistograms) return;
   
   Int_t  clusterID = pCandidate->GetCaloLabel(0) ;
   Int_t  nMaxima   = pCandidate->GetFiducialArea(); // bad name, just place holder for the moment
@@ -1143,110 +1161,163 @@ void AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms(AliA
   if     (fCalorimeter == "EMCAL") clusters = GetEMCALClusters();
   else if(fCalorimeter == "PHOS" ) clusters = GetPHOSClusters();
   
+  if(!clusters) return;
+  
+  AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
+
+  Float_t m02    = cluster->GetM02() ;
   Float_t energy = pCandidate->E();
   Float_t pt     = pCandidate->Pt();
-  
-  if(clusters)
+
+  // Get the max pt leading in cone or the sum of pt in cone
+  // assign a bin to the candidate, depending on both quantities
+  // see the shower shape in those bins.
+  if(fFillBackgroundBinHistograms)
   {
-    AliVCluster *cluster = FindCluster(clusters,clusterID,iclus);
+    // Get the background bin for this cone and trigger
+    Int_t ptsumBin  = -1;
+    Int_t leadptBin = -1;
+    
+    if( GetDebug() > 1 )
+      printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeControlHistograms() - pT cand: %2.2f, In cone pT: Sum %2.2f, Lead %2.2f, n bins %d\n",
+             pt,coneptsum,coneleadpt,fNBkgBin);
     
-    if(fFillSSHisto)
+    for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
     {
-      fhELambda0 [isolated]->Fill(energy, cluster->GetM02() );
-      fhPtLambda0[isolated]->Fill(pt,     cluster->GetM02() );
-      fhELambda1 [isolated]->Fill(energy, cluster->GetM20() );
-      
-      if(IsDataMC())
+      if( coneptsum  >= fBkgBinLimit[ibin] && coneptsum  < fBkgBinLimit[ibin+1]) ptsumBin  = ibin;
+      if( coneleadpt >= fBkgBinLimit[ibin] && coneleadpt < fBkgBinLimit[ibin+1]) leadptBin = ibin;
+    }
+    
+    // Fill the histograms per bin of pt lead or pt sum
+    if( GetDebug() > 1 && ptsumBin  >=0 ) printf("\t Sum bin %d [%2.2f,%2.2f]\n" , ptsumBin ,fBkgBinLimit[ptsumBin] ,fBkgBinLimit[ptsumBin +1]);
+    if( GetDebug() > 1 && leadptBin >=0 ) printf("\t Lead bin %d [%2.2f,%2.2f]\n", leadptBin,fBkgBinLimit[leadptBin],fBkgBinLimit[leadptBin+1]);
+    
+    if( leadptBin >=0 ) fhPtLeadConeBinLambda0[leadptBin]->Fill(pt,m02);
+    if( ptsumBin  >=0 ) fhSumPtConeBinLambda0 [ ptsumBin]->Fill(pt,m02);
+    
+    if( GetDebug() > 1 && leadptBin == 0 )
+      printf("No track/clusters in isolation cone: cand pt %2.2f GeV/c, track multiplicity %d, N clusters %d\n",
+             pt, GetTrackMultiplicity(),GetEMCALClusters()->GetEntriesFast());
+    
+    if(IsDataMC())
+    {
+      Int_t leadptBinMC = leadptBin+mcIndex*fNBkgBin;
+      Int_t  ptsumBinMC =  ptsumBin+mcIndex*fNBkgBin;
+      if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
+      if( ptsumBin  >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
+      if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
       {
-        if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
-          fhPtLambda0MC[kmcPhoton][isolated]->Fill(pt, cluster->GetM02());
-        
-        fhPtLambda0MC[mcIndex][isolated]->Fill(pt, cluster->GetM02());
+        leadptBinMC = leadptBin+kmcPhoton*fNBkgBin;
+        ptsumBinMC  =  ptsumBin+kmcPhoton*fNBkgBin;
+        if( leadptBin >=0 ) fhPtLeadConeBinLambda0MC[leadptBinMC]->Fill(pt,m02);
+        if( ptsumBin  >=0 ) fhSumPtConeBinLambda0MC [ ptsumBinMC]->Fill(pt,m02);
       }
-      
-      if(fCalorimeter == "EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
-         GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD()  )
+    }
+  }
+  
+  if(fFillSSHisto)
+  {
+    fhELambda0 [isolated]->Fill(energy, m02);
+    fhPtLambda0[isolated]->Fill(pt,     m02);
+    fhELambda1 [isolated]->Fill(energy, m02);
+    if(fFillTaggedDecayHistograms)
+    {
+      Int_t decayTag = pCandidate->GetBtag(); // temporary
+      if(decayTag < 0) decayTag = 0;    // temporary
+      for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
       {
-        fhELambda0TRD [isolated]->Fill(energy, cluster->GetM02() );
-        fhPtLambda0TRD[isolated]->Fill(pt    , cluster->GetM02() );
-        fhELambda1TRD [isolated]->Fill(energy, cluster->GetM20() );
+        if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+          fhPtLambda0Decay[isolated][ibit]->Fill(pt,m02);
       }
+    }
+    
+    if(IsDataMC())
+    {
+      if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
+        fhPtLambda0MC[kmcPhoton][isolated]->Fill(pt,m02);
       
-      if(fFillNLMHistograms)
-      {
-        fhNLocMax[isolated]->Fill(energy,nMaxima);
-        if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax1[isolated]->Fill(energy,cluster->GetM20()); }
-        else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMax2[isolated]->Fill(energy,cluster->GetM20()); }
-        else                { fhELambda0LocMaxN[isolated]->Fill(energy,cluster->GetM02()); fhELambda1LocMaxN[isolated]->Fill(energy,cluster->GetM20()); }
-      }
-    } // SS histo fill
+      fhPtLambda0MC[mcIndex][isolated]->Fill(pt,m02);
+    }
+    
+    if(fCalorimeter == "EMCAL" &&  GetFirstSMCoveredByTRD() >= 0 &&
+       GetModuleNumber(cluster) >= GetFirstSMCoveredByTRD()  )
+    {
+      fhELambda0TRD [isolated]->Fill(energy, m02 );
+      fhPtLambda0TRD[isolated]->Fill(pt    , m02 );
+      fhELambda1TRD [isolated]->Fill(energy, m02 );
+    }
+    
+    if(fFillNLMHistograms)
+    {
+      fhNLocMax[isolated]->Fill(energy,nMaxima);
+      if     (nMaxima==1) { fhELambda0LocMax1[isolated]->Fill(energy,m02); fhELambda1LocMax1[isolated]->Fill(energy,m02); }
+      else if(nMaxima==2) { fhELambda0LocMax2[isolated]->Fill(energy,m02); fhELambda1LocMax2[isolated]->Fill(energy,m02); }
+      else                { fhELambda0LocMaxN[isolated]->Fill(energy,m02); fhELambda1LocMaxN[isolated]->Fill(energy,m02); }
+    }
+  } // SS histo fill
+  
+  if(fFillTMHisto)
+  {
+    Float_t dZ  = cluster->GetTrackDz();
+    Float_t dR  = cluster->GetTrackDx();
+    
+    if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
+    {
+      dR = 2000., dZ = 2000.;
+      GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
+    }
     
+    //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
+    if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
+    {
+      fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
+      fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
+      if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
+    }
     
-    if(fFillTMHisto)
+    // Check dEdx and E/p of matched clusters
+    
+    if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
     {
-      Float_t dZ  = cluster->GetTrackDz();
-      Float_t dR  = cluster->GetTrackDx();
       
-      if(cluster->IsEMCAL() && GetCaloUtils()->IsRecalculationOfClusterTrackMatchingOn())
-      {
-        dR = 2000., dZ = 2000.;
-        GetCaloUtils()->GetEMCALRecoUtils()->GetMatchedResiduals(cluster->GetID(),dZ,dR);
-      }
+      AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
       
-      //printf("ParticleIsolation: dPhi %f, dEta %f\n",dR,dZ);
-      if(fhTrackMatchedDEta[isolated] && TMath::Abs(dR) < 999)
+      if(track)
       {
-        fhTrackMatchedDEta[isolated]->Fill(energy,dZ);
-        fhTrackMatchedDPhi[isolated]->Fill(energy,dR);
-        if(energy > 0.5) fhTrackMatchedDEtaDPhi[isolated]->Fill(dZ,dR);
+        Float_t dEdx = track->GetTPCsignal();
+        fhdEdx[isolated]->Fill(cluster->E(), dEdx);
+        
+        Float_t eOverp = cluster->E()/track->P();
+        fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
       }
+      //else
+      //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
       
-      // Check dEdx and E/p of matched clusters
       
-      if(TMath::Abs(dZ) < 0.05 && TMath::Abs(dR) < 0.05)
+      if(IsDataMC())
       {
-        
-        AliVTrack *track = GetCaloUtils()->GetMatchedTrack(cluster, GetReader()->GetInputEvent());
-        
-        if(track)
+        if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
         {
-          Float_t dEdx = track->GetTPCsignal();
-          fhdEdx[isolated]->Fill(cluster->E(), dEdx);
+          if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
+                     GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
+          else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
+          else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
+          else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
           
-          Float_t eOverp = cluster->E()/track->P();
-          fhEOverP[isolated]->Fill(cluster->E(),  eOverp);
         }
-        //else
-        //  printf("AliAnaParticleIsolation::FillTrackMatchingShowerShapeHistograms() - Residual OK but (dR, dZ)= (%2.4f,%2.4f) no track associated WHAT? \n", dR,dZ);
-        
-        
-        if(IsDataMC())
+        else
         {
-          if ( !GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCConversion)  )
-          {
-            if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
-                      GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 2.5 );
-            else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 0.5 );
-            else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 1.5 );
-            else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 3.5 );
-            
-          }
-          else
-          {
-            if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
-                      GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
-            else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
-            else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
-            else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
-          }
-          
-        }  // MC
-        
-      } // match window
+          if       ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPi0)      ||
+                     GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCEta)       ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 6.5 );
+          else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton)    ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 4.5 );
+          else if  ( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCElectron)  ) fhTrackMatchedMCParticle[isolated]->Fill(energy, 5.5 );
+          else                                                                                   fhTrackMatchedMCParticle[isolated]->Fill(energy, 7.5 );
+        }
+      }  // MC
       
-    }// TM histos fill
+    } // match window
     
-  } // clusters array available
+  }// TM histos fill
   
 }
 
@@ -1473,61 +1544,69 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
   // Histograms for tagged candidates as decay
   if(fFillTaggedDecayHistograms)
   {
-    fhPtDecayNoIso  = new TH1F("hPtDecayNoIso",
-                               Form("Number of not isolated leading pi0 decay particles vs #it{p}_{T}, %s",parTitle.Data()),
-                               nptbins,ptmin,ptmax);
-    fhPtDecayNoIso->SetYTitle("#it{counts}");
-    fhPtDecayNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
-    outputContainer->Add(fhPtDecayNoIso) ;
-    
-    fhEtaPhiDecayNoIso  = new TH2F("hEtaPhiDecayNoIso",
-                                   Form("Number of not isolated leading Pi0 decay particles #eta vs #phi, %s",parTitle.Data()),
-                                   netabins,etamin,etamax,nphibins,phimin,phimax);
-    fhEtaPhiDecayNoIso->SetXTitle("#eta");
-    fhEtaPhiDecayNoIso->SetYTitle("#phi");
-    outputContainer->Add(fhEtaPhiDecayNoIso) ;
-    
-    if(!fMakeSeveralIC)
+    for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
     {
-      fhPtDecayIso  = new TH1F("hPtDecayIso",
-                               Form("Number of isolated #pi^{0} decay particles vs #it{p}_{T}, %s",parTitle.Data()),
-                               nptbins,ptmin,ptmax);
-      fhPtDecayIso->SetYTitle("#it{counts}");
-      fhPtDecayIso->SetXTitle("#it{p}_{T}(GeV/#it{c})");
-      outputContainer->Add(fhPtDecayIso) ;
-    
-      fhEtaPhiDecayIso  = new TH2F("hEtaPhiDecayIso",
-                                   Form("Number of isolated Pi0 decay particles #eta vs #phi, %s",parTitle.Data()),
-                                   netabins,etamin,etamax,nphibins,phimin,phimax);
-      fhEtaPhiDecayIso->SetXTitle("#eta");
-      fhEtaPhiDecayIso->SetYTitle("#phi");
-      outputContainer->Add(fhEtaPhiDecayIso) ;
-    }
-    
-    
-    if(IsDataMC())
-    {
-      for(Int_t imc = 0; imc < 9; imc++)
+      fhPtDecayNoIso[ibit]  =
+      new TH1F(Form("hPtDecayNoIso_bit%d",fDecayBits[ibit]),
+               Form("Number of not isolated leading pi0 decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
+               nptbins,ptmin,ptmax);
+      fhPtDecayNoIso[ibit]->SetYTitle("#it{counts}");
+      fhPtDecayNoIso[ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhPtDecayNoIso[ibit]) ;
+      
+      fhEtaPhiDecayNoIso[ibit]  =
+      new TH2F(Form("hEtaPhiDecayNoIso_bit%d",fDecayBits[ibit]),
+               Form("Number of not isolated leading Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
+               netabins,etamin,etamax,nphibins,phimin,phimax);
+      fhEtaPhiDecayNoIso[ibit]->SetXTitle("#eta");
+      fhEtaPhiDecayNoIso[ibit]->SetYTitle("#phi");
+      outputContainer->Add(fhEtaPhiDecayNoIso[ibit]) ;
+      
+      if(!fMakeSeveralIC)
       {
-        
-        fhPtDecayNoIsoMC[imc]  = new TH1F(Form("hPtDecayNoIsoMC%s",mcPartName[imc].Data()),
-                                          Form("#it{p}_{T} of NOT isolated, decay tag,  %s, %s",mcPartType[imc].Data(),parTitle.Data()),
-                                          nptbins,ptmin,ptmax);
-        fhPtDecayNoIsoMC[imc]->SetYTitle("#it{counts}");
-        fhPtDecayNoIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
-        outputContainer->Add(fhPtDecayNoIsoMC[imc]) ;
-        
-        if(!fMakeSeveralIC)
+        fhPtDecayIso[ibit]  =
+        new TH1F(Form("hPtDecayIso_bit%d",fDecayBits[ibit]),
+                 Form("Number of isolated #pi^{0} decay particles vs #it{p}_{T}, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
+                 nptbins,ptmin,ptmax);
+        fhPtDecayIso[ibit]->SetYTitle("#it{counts}");
+        fhPtDecayIso[ibit]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
+        outputContainer->Add(fhPtDecayIso[ibit]) ;
+        
+        fhEtaPhiDecayIso[ibit]  =
+        new TH2F(Form("hEtaPhiDecayIso_bit%d",fDecayBits[ibit]),
+                 Form("Number of isolated Pi0 decay particles #eta vs #phi, bit %d, %s",fDecayBits[ibit],parTitle.Data()),
+                 netabins,etamin,etamax,nphibins,phimin,phimax);
+        fhEtaPhiDecayIso[ibit]->SetXTitle("#eta");
+        fhEtaPhiDecayIso[ibit]->SetYTitle("#phi");
+        outputContainer->Add(fhEtaPhiDecayIso[ibit]) ;
+      }
+      
+      if(IsDataMC())
+      {
+        for(Int_t imc = 0; imc < 9; imc++)
         {
-          fhPtDecayIsoMC[imc]  = new TH1F(Form("hPtDecayMC%s",mcPartName[imc].Data()),
-                                          Form("#it{p}_{T} of isolated %s, decay tag, %s",mcPartType[imc].Data(),parTitle.Data()),
-                                          nptbins,ptmin,ptmax);
-          fhPtDecayIsoMC[imc]->SetYTitle("#it{counts}");
-          fhPtDecayIsoMC[imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
-          outputContainer->Add(fhPtDecayIsoMC[imc]) ;
-        }
-      }// loop
-    }// MC
+          
+          fhPtDecayNoIsoMC[ibit][imc]  =
+          new TH1F(Form("hPtDecayNoIso_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
+                   Form("#it{p}_{T} of NOT isolated, decay bit %d,  %s, %s",fDecayBits[ibit],mcPartType[imc].Data(),parTitle.Data()),
+                   nptbins,ptmin,ptmax);
+          fhPtDecayNoIsoMC[ibit][imc]->SetYTitle("#it{counts}");
+          fhPtDecayNoIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
+          outputContainer->Add(fhPtDecayNoIsoMC[ibit][imc]) ;
+          
+          if(!fMakeSeveralIC)
+          {
+            fhPtDecayIsoMC[ibit][imc]  =
+            new TH1F(Form("hPtDecay_bit%d_MC%s",fDecayBits[ibit],mcPartName[imc].Data()),
+                     Form("#it{p}_{T} of isolated %s, decay bit %d, %s",mcPartType[imc].Data(),fDecayBits[ibit],parTitle.Data()),
+                     nptbins,ptmin,ptmax);
+            fhPtDecayIsoMC[ibit][imc]->SetYTitle("#it{counts}");
+            fhPtDecayIsoMC[ibit][imc]->SetXTitle("#it{p}_{T}(GeV/#it{c})");
+            outputContainer->Add(fhPtDecayIsoMC[ibit][imc]) ;
+          }
+        }// MC particle loop
+      }// MC
+    } // bit loop
   }// decay
   
   if(!fMakeSeveralIC)
@@ -1602,6 +1681,13 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
       fhPtNLocMaxNoIso->SetXTitle("#it{p}_{T} (GeV/#it{c})");
       outputContainer->Add(fhPtNLocMaxNoIso) ;
     }
+
+    fhConePtLead  = new TH2F("hConePtLead",
+                            Form("Track or Cluster  leading #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
+                            nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+    fhConePtLead->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
+    fhConePtLead->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
+    outputContainer->Add(fhConePtLead) ;
     
     fhConeSumPt  = new TH2F("hConePtSum",
                             Form("Track and Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
@@ -1625,6 +1711,61 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
     fhPtInCone->SetXTitle("#it{p}_{T} (GeV/#it{c})");
     outputContainer->Add(fhPtInCone) ;
     
+    if(fFillBackgroundBinHistograms)
+    {
+      fhPtLeadConeBinLambda0 = new TH2F*[fNBkgBin];
+      fhSumPtConeBinLambda0  = new TH2F*[fNBkgBin];
+      
+      if(IsDataMC())
+      {
+        fhPtLeadConeBinLambda0MC = new TH2F*[fNBkgBin*9];
+        fhSumPtConeBinLambda0MC  = new TH2F*[fNBkgBin*9];
+      }
+      
+      for(Int_t ibin = 0; ibin < fNBkgBin; ibin++)
+      {
+        fhPtLeadConeBinLambda0[ibin]  = new TH2F
+        (Form("hPtLeadConeLambda0_Bin%d",ibin),
+         Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), %s",
+              fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+        fhPtLeadConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
+        fhPtLeadConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhPtLeadConeBinLambda0[ibin]) ;
+        
+        fhSumPtConeBinLambda0[ibin]  = new TH2F
+        (Form("hSumPtConeLambda0_Bin%d",ibin),
+         Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), %s",
+              fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+        fhSumPtConeBinLambda0[ibin]->SetYTitle("#lambda_{0}^{2}");
+        fhSumPtConeBinLambda0[ibin]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+        outputContainer->Add(fhSumPtConeBinLambda0[ibin]) ;
+        
+        if(IsDataMC())
+        {
+          for(Int_t imc = 0; imc < 9; imc++)
+          {
+            Int_t binmc = ibin+imc*fNBkgBin;
+            fhPtLeadConeBinLambda0MC[binmc]  = new TH2F
+            (Form("hPtLeadConeLambda0_Bin%d_MC%s",ibin, mcPartName[imc].Data()),
+             Form("#lambda_{0}, in cone %2.2f<#it{p}_{T}^{leading}<%2.2f (GeV/#it{c}), MC %s, %s",
+                  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+            fhPtLeadConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
+            fhPtLeadConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+            outputContainer->Add(fhPtLeadConeBinLambda0MC[binmc]) ;
+            
+            fhSumPtConeBinLambda0MC[binmc]  = new TH2F
+            (Form("hSumPtConeLambda0_Bin%d_MC%s",ibin,mcPartName[imc].Data()),
+             Form("#lambda_{0}, in cone %2.2f <#Sigma #it{p}_{T}< %2.2f (GeV/#it{c}), MC %s, %s",
+                  fBkgBinLimit[ibin],fBkgBinLimit[ibin+1], mcPartType[imc].Data(), parTitle.Data()),nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+            fhSumPtConeBinLambda0MC[binmc]->SetYTitle("#lambda_{0}^{2}");
+            fhSumPtConeBinLambda0MC[binmc]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+            outputContainer->Add(fhSumPtConeBinLambda0MC[binmc]) ;
+          } // MC particle loop
+        }
+        
+      }// pt bin loop
+    } // bkg cone pt bin histograms
+    
     if(fFillHighMultHistograms)
     {
       fhPtInConeCent  = new TH2F("hPtInConeCent",
@@ -1645,6 +1786,14 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
       fhConeSumPtCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
       outputContainer->Add(fhConeSumPtCluster) ;
       
+      fhConePtLeadCluster  = new TH2F("hConeLeadPtCluster",
+                                    Form("Cluster leading in isolation cone for #it{R} =  %2.2f",r),
+                                    nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhConePtLeadCluster->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
+      fhConePtLeadCluster->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
+      outputContainer->Add(fhConePtLeadCluster) ;
+
+      
       if(fFillCellHistograms)
       {
         fhConeSumPtCell  = new TH2F("hConePtSumCell",
@@ -1991,6 +2140,13 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
       fhConeSumPtTrack->SetYTitle("#Sigma #it{p}_{T}");
       fhConeSumPtTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
       outputContainer->Add(fhConeSumPtTrack) ;
+
+      fhConePtLeadTrack  = new TH2F("hConeLeadPtTrack",
+                                   Form("Track leading in isolation cone for #it{R} =  %2.2f",r),
+                                   nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhConePtLeadTrack->SetYTitle("#it{p}_{T, leading} (GeV/#it{c})");
+      fhConePtLeadTrack->SetXTitle("#it{p}_{T, trigger} (GeV/#it{c})");
+      outputContainer->Add(fhConePtLeadTrack) ;
       
       fhPtTrackInCone  = new TH2F("hPtTrackInCone",
                                   Form("#it{p}_{T} of tracks in isolation cone for #it{R} =  %2.2f",r),
@@ -2171,9 +2327,32 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
       fhConeSumPtClustervsTrack   = new TH2F("hConePtSumClustervsTrack",
                                              Form("Track vs Cluster #Sigma #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
                                              nptsumbins,ptsummin,ptsummax,nptsumbins,ptsummin,ptsummax);
-      fhConeSumPtClustervsTrack->SetXTitle("#Sigma #it{p}_{T} cluster");
-      fhConeSumPtClustervsTrack->SetYTitle("#Sigma #it{p}_{T} track");
+      fhConeSumPtClustervsTrack->SetXTitle("#Sigma #it{p}_{T}^{cluster} (GeV/#it{c})");
+      fhConeSumPtClustervsTrack->SetYTitle("#Sigma #it{p}_{T}^{track} (GeV/#it{c})");
       outputContainer->Add(fhConeSumPtClustervsTrack) ;
+
+      fhConeSumPtClusterTrackFrac   = new TH2F("hConePtSumClusterTrackFraction",
+                                             Form("#Sigma #it{p}_{T}^{cluster}/#Sigma #it{p}_{T}^{track} in isolation cone for #it{R} =  %2.2f",r),
+                                             nptbins,ptmin,ptmax,200,0,5);
+      fhConeSumPtClusterTrackFrac->SetYTitle("#Sigma #it{p}^{cluster}_{T} /#Sigma #it{p}_{T}^{track}");
+      fhConeSumPtClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhConeSumPtClusterTrackFrac) ;
+
+      
+      fhConePtLeadClustervsTrack   = new TH2F("hConePtLeadClustervsTrack",
+                                             Form("Track vs Cluster lead #it{p}_{T} in isolation cone for #it{R} =  %2.2f",r),
+                                             nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
+      fhConePtLeadClustervsTrack->SetXTitle("#it{p}^{leading cluster}_{T} (GeV/#it{c})");
+      fhConePtLeadClustervsTrack->SetYTitle("#it{p}^{leading track}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhConePtLeadClustervsTrack) ;
+      
+      fhConePtLeadClusterTrackFrac   = new TH2F("hConePtLeadClusterTrackFraction",
+                                               Form(" #it{p}^{leading cluster}_{T}/#it{p}^{leading track}_{T} in isolation cone for #it{R} =  %2.2f",r),
+                                               nptbins,ptmin,ptmax,200,0,5);
+      fhConePtLeadClusterTrackFrac->SetYTitle("#it{p}^{leading cluster}_{T}/ #it{p}^{leading track}_{T}");
+      fhConePtLeadClusterTrackFrac->SetXTitle("#it{p}^{trigger}_{T} (GeV/#it{c})");
+      outputContainer->Add(fhConePtLeadClusterTrackFrac) ;
+
       
       if(fFillCellHistograms)
       {
@@ -2457,6 +2636,20 @@ TList *  AliAnaParticleIsolation::GetCreateOutputObjects()
         fhPtLambda0[iso]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
         outputContainer->Add(fhPtLambda0[iso]) ;
         
+        if(fFillTaggedDecayHistograms)
+        {
+          for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
+          {
+            fhPtLambda0Decay[iso][ibit]  = new TH2F
+            (Form("hPtLambda0Decay%s_bit%d",isoName[iso].Data(),fDecayBits[ibit]),
+             Form("%s cluster : #it{p}_{T} vs #lambda_{0}, decay bit %d, %s",isoTitle[iso].Data(), fDecayBits[ibit], parTitle.Data()),
+             nptbins,ptmin,ptmax,ssbins,ssmin,ssmax);
+            fhPtLambda0Decay[iso][ibit]->SetYTitle("#lambda_{0}^{2}");
+            fhPtLambda0Decay[iso][ibit]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
+            outputContainer->Add(fhPtLambda0Decay[iso][ibit]) ;
+          }
+        }
+        
         if(IsDataMC())
         {
           for(Int_t imc = 0; imc < 9; imc++)
@@ -2958,7 +3151,8 @@ Int_t AliAnaParticleIsolation::GetMCIndex(Int_t mcTag)
   {
     return kmcPrompt;
   }
-  else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation))
+  else if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCFragmentation) ||
+          GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCISR))
   {
     return kmcFragment;
   }
@@ -3028,6 +3222,18 @@ void AliAnaParticleIsolation::InitParameters()
   fLeadingOnly = kTRUE;
   fCheckLeadingWithNeutralClusters = kTRUE;
   
+  fNDecayBits = 1;
+  fDecayBits[0] = AliNeutralMesonSelection::kPi0;
+  fDecayBits[1] = AliNeutralMesonSelection::kEta;
+  fDecayBits[2] = AliNeutralMesonSelection::kPi0Side;
+  fDecayBits[3] = AliNeutralMesonSelection::kEtaSide;
+  
+  fNBkgBin = 11;
+  fBkgBinLimit[ 0] = 00.0; fBkgBinLimit[ 1] = 00.2; fBkgBinLimit[ 2] = 00.3; fBkgBinLimit[ 3] = 00.4; fBkgBinLimit[ 4] = 00.5;
+  fBkgBinLimit[ 5] = 01.0; fBkgBinLimit[ 6] = 01.5; fBkgBinLimit[ 7] = 02.0; fBkgBinLimit[ 8] = 03.0; fBkgBinLimit[ 9] = 05.0;
+  fBkgBinLimit[10] = 10.0; fBkgBinLimit[11] = 100.;
+  for(Int_t ibin = 12; ibin < 20; ibin++) fBkgBinLimit[ibin] = 00.0;
+  
   //----------- Several IC-----------------
   fNCones             = 5 ;
   fNPtThresFrac       = 5 ;
@@ -3108,8 +3314,17 @@ Bool_t AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle(Int_t &
     
     //skip this event if near side associated particle pt larger than trigger
     
-    if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2())  return kFALSE;
-    
+    Float_t deltaPhi = phiTrig-phi;
+    //
+    // Calculate deltaPhi shift so that for the particles on the opposite side
+    // it is defined between 90 and 270 degrees
+    // Shift [-360,-90]  to [0, 270]
+    // and [270,360] to [-90,0]
+    if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+    if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+
+    if(pt > ptTrig && deltaPhi < TMath::PiOver2())  return kFALSE;
+  
   }// track loop
   
   // Compare if it is leading of all calorimeter clusters
@@ -3143,7 +3358,11 @@ Bool_t AliAnaParticleIsolation::IsTriggerTheNearSideEventLeadingParticle(Int_t &
       // skip this event if near side associated particle pt larger than trigger
       // not really needed for calorimeter, unless DCal is included
      
-      if(pt > ptTrig && TMath::Abs(phi-phiTrig) < TMath::PiOver2()) return kFALSE ;
+      Float_t deltaPhi = phiTrig-phi;
+      if(deltaPhi <= -TMath::PiOver2()) deltaPhi+=TMath::TwoPi();
+      if(deltaPhi > 3*TMath::PiOver2()) deltaPhi-=TMath::TwoPi();
+
+      if(pt > ptTrig && deltaPhi < TMath::PiOver2()) return kFALSE ;
 
     }// cluster loop
   } // check neutral clusters
@@ -3174,7 +3393,7 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
   
   Int_t n = 0, nfrac = 0;
   Bool_t  isolated  = kFALSE ;
-  Float_t coneptsum = 0 ;
+  Float_t coneptsum = 0, coneptlead = 0;
   TObjArray * pl    = 0x0; ;
   
   //Select the calorimeter for candidate isolation with neutral particles
@@ -3209,20 +3428,20 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
   for(Int_t iaod = iaod0; iaod < naod; iaod++)
   {
     AliAODPWG4ParticleCorrelation * aodinput =  (AliAODPWG4ParticleCorrelation*) (GetInputAODBranch()->At(iaod));
-    
-    if(IsLeadingOnlyOn()) aodinput->SetLeadingParticle(kTRUE);
-  
+
     // Check isolation only of clusters in fiducial region
+    
     if(IsFiducialCutOn())
     {
       Bool_t in = GetFiducialCut()->IsInFiducialCut(*aodinput->Momentum(),aodinput->GetDetector()) ;
-      if(! in ) return ;
+      if(! in ) continue ;
     }
     
     //If too small or too large pt, skip
     Float_t pt = aodinput->Pt();
     if(pt < GetMinPt() || pt > GetMaxPt() ) continue ;
 
+    
     //check if it is low pt trigger particle
     if( ( pt < GetIsolationCut()->GetPtThreshold() ||  pt < GetIsolationCut()->GetSumPtThreshold() ) &&
        !fMakeSeveralIC)
@@ -3231,20 +3450,21 @@ void  AliAnaParticleIsolation::MakeAnalysisFillAOD()
     }
     
     //After cuts, study isolation
-    n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0;
+    n=0; nfrac = 0; isolated = kFALSE; coneptsum = 0; coneptlead = 0;
     GetIsolationCut()->MakeIsolationCut(GetCTSTracks(),pl,
                                         GetReader(), GetCaloPID(),
                                         kTRUE, aodinput, GetAODObjArrayName(),
-                                        n,nfrac,coneptsum, isolated);
+                                        n,nfrac,coneptsum,coneptlead,isolated);
     
     if(!fMakeSeveralIC) aodinput->SetIsolated(isolated);
-  } // particle isolationloop
-  
-  if(GetDebug() > 1)
-  {
-    if(isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",idLeading);
-    printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
-  }
+
+    if(GetDebug() > 1)
+    {
+      if(isolated) printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() : Particle %d IS ISOLATED \n",iaod);
+      printf("AliAnaParticleIsolation::MakeAnalysisFillAOD() - End fill AODs \n");
+    }
+
+  } // particle isolation loop
   
 }
 
@@ -3299,7 +3519,7 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
       //In case a more strict IC is needed in the produced AOD
       Bool_t  isolated = kFALSE;
       Int_t   n = 0, nfrac = 0;
-      Float_t coneptsum = 0 ;
+      Float_t coneptsum = 0, coneptlead = 0;
       
       //Recover reference arrays with clusters and tracks
       TObjArray * refclusters = aod->GetObjArray(GetAODObjArrayName()+"Clusters");
@@ -3308,42 +3528,50 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
       GetIsolationCut()->MakeIsolationCut(reftracks,   refclusters,
                                           GetReader(), GetCaloPID(),
                                           kFALSE, aod, "",
-                                          n,nfrac,coneptsum, isolated);
+                                          n,nfrac,coneptsum,coneptlead,isolated);
     }
     
     Bool_t  isolated   = aod->IsIsolated();
-    Bool_t  decay      = aod->IsTagged();
     Float_t energy     = aod->E();
     Float_t phi        = aod->Phi();
     Float_t eta        = aod->Eta();
+
+    Int_t   decayTag = 0;
+    if(fFillTaggedDecayHistograms)
+    {
+      decayTag = aod->GetBtag(); // temporary
+      if(decayTag < 0) decayTag = 0; // temporary
+    }
     
-    if(GetDebug() > 0) printf(" AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f, Isolated %d\n",
-                              pt, eta, phi, isolated);
-    
-    //---------------------------------------------------------------
-    // Fill Shower shape and track matching histograms
-    //---------------------------------------------------------------
-    
-    FillTrackMatchingShowerShapeControlHistograms(aod,mcIndex);
+    if(GetDebug() > 0)
+      printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - pt %1.1f, eta %1.1f, phi %1.1f, Isolated %d\n",
+              pt, eta, phi, isolated);
     
     //---------------------------------------------------------------
     // Fill pt/sum pT distribution of particles in cone or in UE band
     //---------------------------------------------------------------
     
+    Float_t coneptLeadCluster= 0;
+    Float_t coneptLeadTrack  = 0;
     Float_t coneptsumCluster = 0;
     Float_t coneptsumTrack   = 0;
     Float_t coneptsumCell    = 0;
     Float_t etaBandptsumClusterNorm = 0;
     Float_t etaBandptsumTrackNorm   = 0;
     
-    CalculateTrackSignalInCone   (aod,coneptsumTrack  );
-    CalculateCaloSignalInCone    (aod,coneptsumCluster);
+    CalculateTrackSignalInCone   (aod,coneptsumTrack  , coneptLeadTrack  );
+    CalculateCaloSignalInCone    (aod,coneptsumCluster, coneptLeadCluster);
     if(fFillCellHistograms)
       CalculateCaloCellSignalInCone(aod,coneptsumCell   );
     
     if(GetIsolationCut()->GetParticleTypeInCone()==AliIsolationCut::kNeutralAndCharged)
     {
-      fhConeSumPtClustervsTrack     ->Fill(coneptsumCluster,coneptsumTrack);
+      fhConeSumPtClustervsTrack ->Fill(coneptsumCluster, coneptsumTrack );
+      fhConePtLeadClustervsTrack->Fill(coneptLeadCluster,coneptLeadTrack);
+
+      if(coneptsumTrack  > 0) fhConeSumPtClusterTrackFrac ->Fill(pt, coneptsumCluster /coneptsumTrack );
+      if(coneptLeadTrack > 0) fhConePtLeadClusterTrackFrac->Fill(pt, coneptLeadCluster/coneptLeadTrack);
+
       if(fFillCellHistograms)
       {
         fhConeSumPtCellvsTrack        ->Fill(coneptsumCell,   coneptsumTrack);
@@ -3355,8 +3583,13 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
     fhConeSumPt              ->Fill(pt,     coneptsumTrack+coneptsumCluster);
     fhConeSumPtTrigEtaPhi    ->Fill(eta,phi,coneptsumTrack+coneptsumCluster);
     
+    Float_t coneptLead = coneptLeadTrack;
+    if(coneptLeadCluster > coneptLeadTrack) coneptLead = coneptLeadCluster;
+    fhConePtLead->Fill(pt, coneptLead );
+    
     if(GetDebug() > 1)
-      printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f\n", iaod, coneptsumTrack+coneptsumCluster);
+      printf("AliAnaParticleIsolation::MakeAnalysisFillHistograms() - Particle %d Energy Sum in Isolation Cone %2.2f, Leading pT in cone %2.2f\n",
+             iaod, coneptsumTrack+coneptsumCluster, coneptLead);
     
     //normalize phi/eta band per area unit
     if(fFillUEBandSubtractHistograms)
@@ -3364,6 +3597,11 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
     
     //  printf("Histograms analysis : cluster pt = %f, etaBandTrack = %f, etaBandCluster = %f, isolation = %d\n",aod->Pt(),etaBandptsumTrackNorm,etaBandptsumClusterNorm,aod->IsIsolated());
     
+    //---------------------------------------------------------------
+    // Fill Shower shape and track matching histograms
+    //---------------------------------------------------------------
+    
+    FillTrackMatchingShowerShapeControlHistograms(aod, coneptsumTrack+coneptsumCluster, coneptLead, mcIndex);
     
     //---------------------------------------------------------------
     // Isolated/ Non isolated histograms
@@ -3396,20 +3634,26 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
       }//Histograms with MC
       
       // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
-      if(decay && fFillTaggedDecayHistograms)
+      if(fFillTaggedDecayHistograms && decayTag > 0)
       {
-        fhPtDecayIso    ->Fill(pt);
-        fhEtaPhiDecayIso->Fill(eta,phi);
-        
-        if(IsDataMC())
+        for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
         {
-          if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
-            fhPtDecayIsoMC[kmcPhoton]->Fill(pt);
-          
-          fhPtDecayIsoMC[mcIndex]->Fill(pt);
-        }
-      }
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+          {
+            fhPtDecayIso       [ibit]->Fill(pt);
+            fhEtaPhiDecayIso   [ibit]->Fill(eta,phi);
 
+            if(IsDataMC())
+            {
+              if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
+                fhPtDecayIsoMC[ibit][kmcPhoton]->Fill(pt);
+              
+              fhPtDecayIsoMC[ibit][mcIndex]->Fill(pt);
+            }
+          } // bit ok
+        } // bit loop
+      } // decay histograms
+      
       if(fFillNLMHistograms)
         fhPtNLocMaxIso ->Fill(pt,aod->GetFiducialArea()) ; // remember to change method name
       
@@ -3453,19 +3697,25 @@ void  AliAnaParticleIsolation::MakeAnalysisFillHistograms()
       }
       
       // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
-      if(decay && fFillTaggedDecayHistograms)
+      if(fFillTaggedDecayHistograms && decayTag > 0)
       {
-        fhPtDecayNoIso    ->Fill(pt);
-        fhEtaPhiDecayNoIso->Fill(eta,phi);
-        
-        if(IsDataMC())
+        for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
         {
-          if( GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton) )
-            fhPtDecayNoIsoMC[kmcPhoton]->Fill(pt);
-          
-          fhPtDecayNoIsoMC[mcIndex]->Fill(pt);
-        }
-      }// decay
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+          {
+            fhPtDecayNoIso[ibit]    ->Fill(pt);
+            fhEtaPhiDecayNoIso[ibit]->Fill(eta,phi);
+            
+            if(IsDataMC())
+            {
+              if(GetMCAnalysisUtils()->CheckTagBit(mcTag,AliMCAnalysisUtils::kMCPhoton))
+                fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(pt);
+              
+              fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(pt);
+            }
+          } // bit ok
+        } // bit loop
+      } // decay histograms
       
       if(fFillNLMHistograms)
         fhPtNLocMaxNoIso ->Fill(pt,aod->GetFiducialArea()); // remember to change method name
@@ -3752,10 +4002,16 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
   Float_t etaC  = ph->Eta();
   Float_t phiC  = ph->Phi();
   Int_t   tag   = ph->GetTag();
-  Bool_t  decay = ph->IsTagged();
-  
+
+  Int_t   decayTag = 0;
+  if(fFillTaggedDecayHistograms)
+  {
+    decayTag = ph->GetBtag(); // temporary
+    if(decayTag < 0) decayTag = 0; // temporary
+  }
+
   if(GetDebug() > 0)
-    printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f, decay? %d\n",ptC, decay);
+    printf("AliAnaParticleIsolation::MakeSeveralICAnalysis() - Isolate pT %2.2f, decay tag %d\n",ptC, decayTag);
   
   //Keep original setting used when filling AODs, reset at end of analysis
   Float_t ptthresorg = GetIsolationCut()->GetPtThreshold();
@@ -3763,7 +4019,7 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
   Float_t ptsumcorg  = GetIsolationCut()->GetSumPtThreshold();
   Float_t rorg       = GetIsolationCut()->GetConeSize();
   
-  Float_t coneptsum = 0 ;
+  Float_t coneptsum = 0, coneptlead = 0;
   Int_t   n    [10][10];//[fNCones][fNPtThresFrac];
   Int_t   nfrac[10][10];//[fNCones][fNPtThresFrac];
   Bool_t  isolated  = kFALSE;
@@ -3782,19 +4038,25 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
   }
   
   // Candidates tagged as decay in another analysis (AliAnaPi0EbE)
-  if(decay && fFillTaggedDecayHistograms)
+  if(fFillTaggedDecayHistograms && decayTag > 0)
   {
-    fhPtDecayNoIso    ->Fill(ptC);
-    fhEtaPhiDecayNoIso->Fill(etaC,phiC);
-    
-    if(IsDataMC())
+    for(Int_t ibit = 0; ibit < fNDecayBits; ibit++)
     {
-      if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
-        fhPtDecayNoIsoMC[kmcPhoton]->Fill(ptC);
-      
-      fhPtDecayNoIsoMC[mcIndex]->Fill(ptC);
-    }
-  }
+      if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[ibit]))
+      {
+        fhPtDecayNoIso[ibit]    ->Fill(ptC);
+        fhEtaPhiDecayNoIso[ibit]->Fill(etaC,phiC);
+        
+        if(IsDataMC())
+        {
+          if(GetMCAnalysisUtils()->CheckTagBit(tag,AliMCAnalysisUtils::kMCPhoton))
+            fhPtDecayNoIsoMC[ibit][kmcPhoton]->Fill(ptC);
+          
+          fhPtDecayNoIsoMC[ibit][mcIndex]->Fill(ptC);
+        }
+      } // bit ok
+    } // bit loop
+  } // decay histograms
   
   //Get vertex for photon momentum calculation
   Double_t vertex[] = {0,0,0} ; //vertex ;
@@ -3813,7 +4075,7 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
     
     //In case a more strict IC is needed in the produced AOD
     
-    isolated = kFALSE; coneptsum = 0;
+    isolated = kFALSE; coneptsum = 0; coneptlead = 0;
     
     GetIsolationCut()->SetSumPtThreshold(100);
     GetIsolationCut()->SetPtThreshold(100);
@@ -3916,7 +4178,8 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
       GetIsolationCut()->MakeIsolationCut(reftracks, refclusters,
                                           GetReader(), GetCaloPID(),
                                           kFALSE, ph, "",
-                                          n[icone][ipt],nfrac[icone][ipt],coneptsum, isolated);
+                                          n[icone][ipt],nfrac[icone][ipt],
+                                          coneptsum, coneptlead, isolated);
       
       // Normal pT threshold cut
       
@@ -3938,10 +4201,13 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
         fhPtThresIsolated [icone][ipt]->Fill(ptC);
         fhEtaPhiPtThresIso[icone][ipt]->Fill(etaC,phiC);
         
-        if( decay &&  fFillTaggedDecayHistograms )
+        if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
         {
-          fhPtPtThresDecayIso    [icone][ipt]->Fill(ptC);
-          fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
+          {
+            fhPtPtThresDecayIso    [icone][ipt]->Fill(ptC);
+            fhEtaPhiPtThresDecayIso[icone][ipt]->Fill(etaC,phiC);
+          }
         }
         
         if(IsDataMC())
@@ -3963,10 +4229,13 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
         fhPtFracIsolated [icone][ipt]->Fill(ptC);
         fhEtaPhiPtFracIso[icone][ipt]->Fill(etaC,phiC);
         
-        if( decay &&  fFillTaggedDecayHistograms )
+        if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
         {
-          fhPtPtFracDecayIso    [icone][ipt]->Fill(ptC);
-          fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
+          {
+            fhPtPtFracDecayIso    [icone][ipt]->Fill(ptC);
+            fhEtaPhiPtFracDecayIso[icone][ipt]->Fill(etaC,phiC);
+          }
         }
         
         if(IsDataMC())
@@ -3990,10 +4259,13 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
         fhSumPtIsolated [icone][ipt]->Fill(ptC) ;
         fhEtaPhiPtSumIso[icone][ipt]->Fill(etaC, phiC) ;
         
-        if( decay && fFillTaggedDecayHistograms )
+        if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
         {
-          fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
-          fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
+          {
+            fhPtPtSumDecayIso[icone][ipt]->Fill(ptC);
+            fhEtaPhiPtSumDecayIso[icone][ipt]->Fill(etaC, phiC) ;
+          }
         }
         
         if(IsDataMC())
@@ -4016,10 +4288,13 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
         fhPtFracPtSumIso    [icone][ipt]->Fill(ptC) ;
         fhEtaPhiFracPtSumIso[icone][ipt]->Fill(etaC,phiC) ;
         
-        if( decay &&  fFillTaggedDecayHistograms )
+        if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
         {
-          fhPtFracPtSumDecayIso    [icone][ipt]->Fill(ptC);
-          fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
+          {
+            fhPtFracPtSumDecayIso    [icone][ipt]->Fill(ptC);
+            fhEtaPhiFracPtSumDecayIso[icone][ipt]->Fill(etaC,phiC);
+          }
         }
       }
       
@@ -4033,12 +4308,14 @@ void  AliAnaParticleIsolation::MakeSeveralICAnalysis(AliAODPWG4ParticleCorrelati
         fhPtSumDensityIso    [icone][ipt]->Fill(ptC) ;
         fhEtaPhiSumDensityIso[icone][ipt]->Fill(etaC,phiC) ;
         
-        if( decay &&  fFillTaggedDecayHistograms )
+        if( fFillTaggedDecayHistograms && decayTag > 0 && fNDecayBits > 0)
         {
-          fhPtSumDensityDecayIso    [icone][ipt]->Fill(ptC);
-          fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
+          if(GetNeutralMesonSelection()->CheckDecayBit(decayTag,fDecayBits[0]))
+          {
+            fhPtSumDensityDecayIso    [icone][ipt]->Fill(ptC);
+            fhEtaPhiSumDensityDecayIso[icone][ipt]->Fill(etaC, phiC);
+          }
         }
-        
       }
     }//pt thresh loop