Index: /trunk/Delphes.cpp
===================================================================
--- /trunk/Delphes.cpp	(revision 320)
+++ /trunk/Delphes.cpp	(revision 321)
@@ -571,7 +571,5 @@
 	  // sorts the vector and smears duplicates 
       list_of_towers_with_photon.mergeDuplicates(); 
-      
-      for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) 
-	{
+      for(unsigned int i=0; i<list_of_towers_with_photon.size(); i++) { 
 	  float eta = list_of_towers_with_photon[i].getEta(); 
 	  float phi = list_of_towers_with_photon[i].getPhi();
@@ -588,5 +586,5 @@
       // 2.2.1 ********************* sorting collections by decreasing pt
       DET->SortedVector(electron); 
-      float iPhiEl=0,iEtaEl=0;
+      float iPhiEl=0,iEtaEl=0,ptisoEl=0;
       for(unsigned int i=0; i < electron.size(); i++) 
 	{
@@ -596,5 +594,6 @@
 	  elementElec->PhiCalo = electron[i].PhiCalo();
 	  elementElec->Charge = sign(electron[i].PID());
-	  elementElec->IsolFlag = DET->Isolation(electron[i].Phi(),electron[i].Eta(),TrackCentral,DET->ISOL_PT,DET->ISOL_Cone);      
+	  elementElec->IsolFlag = DET->Isolation(electron[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoEl);
+	  elementElec->IsolPt = ptisoEl;
 	  DET->BinEtaPhi(elementElec->PhiCalo,elementElec->EtaCalo,iPhiEl,iEtaEl); 
  	  D_CaloTower calElec(list_of_active_towers.getElement(iEtaEl,iPhiEl));
@@ -603,5 +602,5 @@
       
       DET->SortedVector(muon);
-      float iPhiMu=0,iEtaMu=0;
+      float iPhiMu=0,iEtaMu=0,ptisoMu=0;
       for(unsigned int i=0; i < muon.size(); i++) 
 	{
@@ -611,8 +610,11 @@
 	  elementMu->EtaCalo = muon[i].EtaCalo(); 
 	  elementMu->PhiCalo = muon[i].PhiCalo();
-	  elementMu->IsolFlag = DET->Isolation(muon[i].Phi(),muon[i].Eta(),TrackCentral,DET->ISOL_PT,DET->ISOL_Cone); 
-	  DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu); 
+	  elementMu->IsolFlag = DET->Isolation(muon[i],TrackCentral,DET->ISOL_PT,DET->ISOL_Cone,ptisoMu); 
+	  elementMu->IsolPt = ptisoMu;
+	  DET->BinEtaPhi(elementMu->PhiCalo,elementMu->EtaCalo,iPhiMu,iEtaMu);
 	  D_CaloTower calMuon(list_of_active_towers.getElement(iEtaMu,iPhiMu));
-	  elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
+	  if( calMuon.getEem() !=0 ) elementMu->EHoverEE = calMuon.getEhad()/calMuon.getEem();
+	  else elementMu->EHoverEE = UNDEFINED;
+	  elementMu->EtRatio = DET->CaloIsolation(muon[i], list_of_active_towers);
 	}
 
Index: /trunk/interface/SmearUtil.h
===================================================================
--- /trunk/interface/SmearUtil.h	(revision 320)
+++ /trunk/interface/SmearUtil.h	(revision 321)
@@ -155,4 +155,7 @@
   float ISOL_PT;      //minimal pt of tracks for isolation criteria
   float ISOL_Cone;    //Cone  for isolation criteria
+  float ISOL_Calo_ET;  //minimal tower energy for isolation criteria
+  float ISOL_Calo_Cone;//Cone for calorimetric isolation
+  unsigned int ISOL_Calo_Grid; //Grid size (N x N) for calorimetric isolation
   
   //General jet variable
@@ -230,6 +233,9 @@
   bool Btaggedjet(const TLorentzVector &JET, const TSimpleArray<GenParticle> &subarray);
 
-  /// Lepton isolation
-  bool Isolation(const float phi, const float eta,const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone);
+  /// Lepton isolation based on tracking
+  bool Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso);
+
+  /// Lepton isolation based on calorimetry (optional. Default: off)
+  float CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers);
 
   //********************* returns a segmented value for eta and phi, for calo towers *****
Index: /trunk/src/LHCOConverter.cc
===================================================================
--- /trunk/src/LHCOConverter.cc	(revision 320)
+++ /trunk/src/LHCOConverter.cc	(revision 321)
@@ -205,11 +205,7 @@
 	  case lhcoMuonID: {
 	   cout << "ERROR: LHCOConverter: for muons use LHCOConverter::BranchReaderMuon!\n";
-		TRootMuon muon(*((TRootMuon*) branch->At(i))); 
-		ntrk = muon.Charge;  ratioE = muon.EHoverEE; 
 	  } break;
 	  case lhcoETmisID: {
 	   cout << "ERROR: LHCOConverter: for ETmis use LHCOConverter::BranchReaderETmis!\n";
-		TRootMuon muon(*((TRootMuon*) branch->At(i))); 
-		ntrk = muon.Charge;  ratioE = muon.EHoverEE; 
 	  } break;
 
@@ -284,4 +280,10 @@
         } // loop on jets
 
+	int ratio = 0; char ratio_string[3];
+	if(muon->EtRatio != UNDEFINED) ratio =  static_cast<int>(muon->EtRatio);
+        if(ratio<0) sprintf(ratio_string,"00");
+	else if(ratio<10) sprintf(ratio_string,"0%d",ratio);
+	else if(ratio<100)sprintf(ratio_string,"%d",ratio);
+	else { cout << "Error: muon EtRatio > 99" << endl; sprintf(ratio_string,"99"); }
 	// output
 	outfile << fixed << setprecision(3)
@@ -293,6 +295,8 @@
                 << setw(7) << jmass <<" "	// invariant mass
                 << setw(7) << (float)muon->Charge  <<" "	// -1 for mu-, +1 for mu+
-                << setw(7) << (float)closest_jet_ID << " "
-                << setw(7) << muon->EHoverEE <<" "	// E_had/E_em
+                << setw(7) << (float)closest_jet_ID << " ";
+	outfile << fixed << setprecision(0)
+                << setw(4) << floor(muon->IsolPt) << "." << ratio_string<<" ";	// E_had/E_em
+	outfile << fixed << setprecision(3)
                 << setw(7) << 0.0 << setw(7) << 0.0  // dummy/dummy
                 << endl;
Index: /trunk/src/SmearUtil.cc
===================================================================
--- /trunk/src/SmearUtil.cc	(revision 320)
+++ /trunk/src/SmearUtil.cc	(revision 321)
@@ -102,7 +102,10 @@
   PTCUT_taujet    = 10.0;
 
+  // Isolation
   ISOL_PT         = 2.0;      //minimal pt of tracks for isolation criteria
   ISOL_Cone       = 0.5;      //Cone  for isolation criteria
-
+  ISOL_Calo_ET    = 1E99;     //minimal tower energy for isolation criteria. Default off = 1E99
+  ISOL_Calo_Cone  = 0.5;      //Cone for calorimetric isolation
+  ISOL_Calo_Grid  = 3;        //Grid size (N x N) for calorimetric isolation
 
   // General jet variable
@@ -228,6 +231,10 @@
   PTCUT_taujet    = DET.PTCUT_taujet;
 
-  ISOL_PT         = DET.ISOL_PT;      //minimal pt of tracks for isolation criteria
-  ISOL_Cone       = DET.ISOL_Cone;      //Cone  for isolation criteria
+  // Isolation
+  ISOL_PT         = DET.ISOL_PT;     // tracking isolation
+  ISOL_Cone       = DET.ISOL_Cone;  
+  ISOL_Calo_ET    = DET.ISOL_Calo_ET;  // calorimeter isolation, defaut off
+  ISOL_Calo_Cone  = DET.ISOL_Calo_Cone;
+  ISOL_Calo_Grid  = DET.ISOL_Calo_Grid; 
 
 
@@ -347,7 +354,10 @@
   PTCUT_taujet    = DET.PTCUT_taujet;
 
-  ISOL_PT         = DET.ISOL_PT;      //minimal pt of tracks for isolation criteria
-  ISOL_Cone       = DET.ISOL_Cone;      //Cone  for isolation criteria
-
+  // Isolation
+  ISOL_PT         = DET.ISOL_PT;       // tracking isolation
+  ISOL_Cone       = DET.ISOL_Cone;    
+  ISOL_Calo_ET    = DET.ISOL_Calo_ET;  // calorimeter isolation, defaut off
+  ISOL_Calo_Cone  = DET.ISOL_Calo_Cone;
+  ISOL_Calo_Grid  = DET.ISOL_Calo_Grid;
 
   // General jet variable
@@ -497,7 +507,9 @@
     else if(strstr(temp_string.c_str(),"PTCUT_taujet"))     {curstring >> varname >> value; PTCUT_taujet      = value;}
 
-    else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT        = value;}
-    else if(strstr(temp_string.c_str(),"ISOL_Cone"))        {curstring >> varname >> value; ISOL_Cone        = value;}
-
+    else if(strstr(temp_string.c_str(),"ISOL_PT"))          {curstring >> varname >> value; ISOL_PT           = value;}
+    else if(strstr(temp_string.c_str(),"ISOL_Cone"))        {curstring >> varname >> value; ISOL_Cone         = value;}
+    else if(strstr(temp_string.c_str(),"ISOL_Calo_ET"))     {curstring >> varname >> value; ISOL_Calo_ET      = value;}
+    else if(strstr(temp_string.c_str(),"ISOL_Calo_Cone"))   {curstring >> varname >> value; ISOL_Calo_Cone    = value;}
+    else if(strstr(temp_string.c_str(),"ISOL_Calo_Grid"))   {curstring >> varname >> ivalue; ISOL_Calo_Grid    = ivalue;}
 
     else if(strstr(temp_string.c_str(),"JET_coneradius"))   {curstring >> varname >> value; JET_coneradius    = value;}
@@ -787,4 +799,16 @@
   f_out << left << setw(40) <<"* Cone for isolation criteria: "<<""
         << left << setw(20) <<ISOL_Cone             <<""<< right << setw(10)<<"*"<<"\n";
+
+  if(ISOL_Calo_ET > 1E98) f_out<<"# No Calorimetric isolation applied                                  *"<<"\n";
+  else {
+   f_out << left << setw(40) <<"* Minimum ET for towers [GeV]: "<<""
+         << left << setw(20) <<ISOL_Calo_ET               <<""<< right << setw(10)<<"*"<<"\n";
+   f_out << left << setw(40) <<"* Cone for calorimetric isolation: "<<""
+         << left << setw(20) <<ISOL_Calo_Cone             <<""<< right << setw(10)<<"*"<<"\n";
+   f_out << left << setw(40) <<"* Grid size (NxN) for calorimetric isolation: "<<""
+         << left << setw(20) <<ISOL_Calo_Grid             <<""<< right << setw(10)<<"*"<<"\n";
+  }
+
+
   f_out<<"*                                                                    *"<<"\n";
   f_out<<"#***************                                                     *"<<"\n";
@@ -1072,21 +1096,48 @@
 //***********************Isolation criteria***********************
 //****************************************************************
-bool RESOLution::Isolation(const float phi, const float eta,const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone)
+bool RESOLution::Isolation(const D_Particle& part, const vector<TRootTracks> &tracks, const float& pt_second_track, const float& isolCone, float& ptiso )
 {
    bool isolated = false;
-   float deltar=5000.; // Initial value; should be high; no further repercussion
-   // loop on all final charged particles, with p_t >2, close enough from the electron
-   for(unsigned int i=0; i < tracks.size(); i++) 
-      {
-       if(tracks[i].PT < pt_second_track)continue;
-       float genDeltaR = DeltaR(phi,eta,tracks[i].Phi,tracks[i].Eta); 
+   ptiso = 0; // sum of all track pt in isolation cone
+   float deltar=1E99; // Initial value; should be high; no further repercussion
+
+   // loop on all tracks, with p_t above threshold, close enough from the charged lepton
+   for(unsigned int i=0; i < tracks.size(); i++) {
+       if(tracks[i].PT < pt_second_track) continue; // ptcut on tracks
+       float genDeltaR = DeltaR(part.Phi(),part.Eta(),tracks[i].Phi,tracks[i].Eta); 
          if(
              (genDeltaR > deltar) ||
-             (genDeltaR==0)
+             (genDeltaR==0)  // rejets the track of the particle itself
            ) continue ;
-           deltar=genDeltaR;
+           deltar=genDeltaR;	// finds the closest track
+
+         // as long as (genDeltaR==0) is put above, the particle itself is not taken into account
+	 if( genDeltaR < ISOL_Cone) ptiso += tracks[i].PT; // dR cut on tracks
       }
    if(deltar > isolCone) isolated = true;
    return isolated;
+}
+
+// ******* Calorimetric isolation
+float RESOLution::CaloIsolation(const D_Particle& part,  const D_CaloTowerList & towers) {
+ // etrat, which is a percentage between 00 and 99. It is the ratio of the transverse energy 
+ // in a 3Ã3 grid surrounding the muon to the pT of the muon. For well-isolated muons, both ptiso and etrat will be small.
+   if(ISOL_Calo_ET>1E10) return UNDEFINED; // avoid doing anything unreasonable...
+   float etrat=0;
+   // available parameters: ISOL_Calo_ET , ISOL_Calo_Cone , 
+/*   for(unsigned int i=0; i < towers.size(); i++) {
+     if(towers[i].E > ISOL_Calo_ET) {
+        float genDeltaR = DeltaR(part.Phi(),part.Eta(),towers[i].getPhi(),towers[i].getEta());
+        if(genDeltaR < ISOL_Calo_Cone) {
+           ptiso += towers[i].getET();
+        }	    
+     }
+   } // loop on towers
+   ptiso -=
+*/
+  etrat = 100*etrat/part.Pt();
+  if(etrat<0) cout << "Error: negative etrat in CaloIsolation (" << etrat <<")\n";
+  else if(etrat>99) cout << "Error: etrat shoud be in [0;99] in CaloIsolation (" << etrat <<")\n";
+  return etrat; 
 }
 
