Index: modules/Calorimeter.cc
===================================================================
--- modules/Calorimeter.cc	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ modules/Calorimeter.cc	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
@@ -41,10 +41,11 @@
 
 Calorimeter::Calorimeter() :
-  fResolutionFormula(0),
+  fECalResolutionFormula(0), fHCalResolutionFormula(0),
   fItParticleInputArray(0), fItTrackInputArray(0),
   fTowerTrackArray(0), fItTowerTrackArray(0)
 {
-  fResolutionFormula = new DelphesFormula;
- 
+  fECalResolutionFormula = new DelphesFormula;
+  fHCalResolutionFormula = new DelphesFormula;
+
   fTowerTrackArray = new TObjArray;
   fItTowerTrackArray = fTowerTrackArray->MakeIterator();
@@ -55,6 +56,7 @@
 Calorimeter::~Calorimeter()
 {
-  if(fResolutionFormula) delete fResolutionFormula;
- 
+  if(fECalResolutionFormula) delete fECalResolutionFormula;
+  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
+
   if(fTowerTrackArray) delete fTowerTrackArray;
   if(fItTowerTrackArray) delete fItTowerTrackArray;
@@ -66,6 +68,6 @@
 {
   ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
-  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
-  Double_t fraction;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins, sizeFractions;
+  Double_t ecalFraction, hcalFraction;
   TBinMap::iterator itEtaBin;
   set< Double_t >::iterator itPhiBin;
@@ -114,11 +116,15 @@
   // set default energy fractions values
   fFractionMap.clear();
-  fFractionMap[0] = 1.0;
+  fFractionMap[0] = make_pair(0.0, 1.0);
 
   for(i = 0; i < size/2; ++i)
   {
     paramFractions = param[i*2 + 1];
-    fraction = paramFractions[0].GetDouble();
-    fFractionMap[param[i*2].GetInt()] = fraction;
+    sizeFractions = paramFractions.GetSize();
+
+    ecalFraction = paramFractions[0].GetDouble();
+    hcalFraction = paramFractions[1].GetDouble();
+
+    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
   }
 /*
@@ -130,6 +136,7 @@
 */
   // read resolution formulas
-  fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
- 
+  fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
+  fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "0"));
+
   // import array with output from other modules
   fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
@@ -141,6 +148,11 @@
   // create output arrays
   fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
-  fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
- 
+  fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
+  
+  fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
+  fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
+  fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
+
+
 }
 
@@ -167,6 +179,6 @@
   Int_t number;
   Long64_t towerHit, towerEtaPhi, hitEtaPhi;
-  Double_t fraction;
-  Double_t energy;
+  Double_t ecalFraction, hcalFraction;
+  Double_t ecalEnergy, hcalEnergy;
   Int_t pdgCode;
 
@@ -181,7 +193,9 @@
   DelphesFactory *factory = GetFactory();
   fTowerHits.clear();
-  fTowerFractions.clear();
-  fTrackFractions.clear();
- 
+  fTowerECalFractions.clear();
+  fTowerHCalFractions.clear();
+  fTrackECalFractions.clear();
+  fTrackHCalFractions.clear();
+
   // loop over all particles
   fItParticleInputArray->Reset();
@@ -200,8 +214,11 @@
     }
 
-    fraction = itFractionMap->second;
-    fTowerFractions.push_back(fraction);
-   
-    if(fraction < 1.0E-9) continue;
+    ecalFraction = itFractionMap->second.first;
+    hcalFraction = itFractionMap->second.second;
+
+    fTowerECalFractions.push_back(ecalFraction);
+    fTowerHCalFractions.push_back(hcalFraction);
+
+    if(ecalFraction < 1.0E-9 && hcalFraction < 1.0E-9) continue;
 
     // find eta bin [1, fEtaBins.size - 1]
@@ -243,8 +260,10 @@
     }
 
-    fraction = itFractionMap->second;
-  
-    fTrackFractions.push_back(fraction);
-  
+    ecalFraction = itFractionMap->second.first;
+    hcalFraction = itFractionMap->second.second;
+
+    fTrackECalFractions.push_back(ecalFraction);
+    fTrackHCalFractions.push_back(hcalFraction);
+
     // find eta bin [1, fEtaBins.size - 1]
     itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
@@ -308,12 +327,19 @@
       fTowerEdges[3] = (*phiBins)[phiBin];
 
-      fTowerEnergy = 0.0;
-      fTrackEnergy = 0.0;
-      
-      fTowerTime = 0.0;
-      fTrackTime = 0.0;
+      fTowerECalEnergy = 0.0;
+      fTowerHCalEnergy = 0.0;
+
+      fTrackECalEnergy = 0.0;
+      fTrackHCalEnergy = 0.0;
+
+      fTowerECalTime = 0.0;
+      fTowerHCalTime = 0.0;
+
+      fTrackECalTime = 0.0;
+      fTrackHCalTime = 0.0;
+
+      fTowerECalWeightTime = 0.0; 
+      fTowerHCalWeightTime = 0.0;
      
-      fTowerWeightTime = 0.0; 
-      
       fTowerTrackHits = 0;
       fTowerPhotonHits = 0;
@@ -330,12 +356,18 @@
       momentum = track->Momentum;
       position = track->Position;
- 
-      energy = momentum.E() * fTrackFractions[number];
+
       
-      fTrackEnergy += energy;
+      ecalEnergy = momentum.E() * fTrackECalFractions[number];
+      hcalEnergy = momentum.E() * fTrackHCalFractions[number];
+
+      fTrackECalEnergy += ecalEnergy;
+      fTrackHCalEnergy += hcalEnergy;
       
-      fTrackTime += TMath::Sqrt(energy)*position.T();
-      fTrackWeightTime += TMath::Sqrt(energy);
-    
+      fTrackECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+      fTrackHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+       
+      fTrackECalWeightTime += TMath::Sqrt(ecalEnergy);
+      fTrackHCalWeightTime += TMath::Sqrt(hcalEnergy);
+
       fTowerTrackArray->Add(track);
 
@@ -351,11 +383,17 @@
 
     // fill current tower
-    energy = momentum.E() * fTowerFractions[number];
-   
-    fTowerEnergy += energy;
+    ecalEnergy = momentum.E() * fTowerECalFractions[number];
+    hcalEnergy = momentum.E() * fTowerHCalFractions[number];
+
+    fTowerECalEnergy += ecalEnergy;
+    fTowerHCalEnergy += hcalEnergy;
+
+    fTowerECalTime += TMath::Sqrt(ecalEnergy)*position.T();
+    fTowerHCalTime += TMath::Sqrt(hcalEnergy)*position.T();
+
+    fTowerECalWeightTime += TMath::Sqrt(ecalEnergy);
+    fTowerHCalWeightTime += TMath::Sqrt(hcalEnergy);
     
-    fTowerTime += TMath::Sqrt(energy)*position.T();
-    fTowerWeightTime += TMath::Sqrt(energy);
-    
+
     fTower->AddCandidate(particle);
   }
@@ -369,18 +407,33 @@
 void Calorimeter::FinalizeTower()
 {
-  Candidate *tower;
+  Candidate *track, *tower;
   Double_t energy, pt, eta, phi;
-  Double_t sigma;
-  Double_t time;
+  Double_t ecalEnergy, hcalEnergy;
+  Double_t ecalSigma, hcalSigma;
+  Double_t ecalTime, hcalTime, time;
 
   if(!fTower) return;
 
-  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
-
-//  energy = gRandom->Gaus(fTowerEnergy, sigma);
-//  if(energy < 0.0) energy = 0.0;
-
-  energy = LogNormal(fTowerEnergy, sigma);
-  time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime;
+  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerECalEnergy);
+
+//  ecalEnergy = gRandom->Gaus(fTowerECalEnergy, ecalSigma);
+//  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
+
+  ecalEnergy = LogNormal(fTowerECalEnergy, ecalSigma);
+  ecalTime = (fTowerECalWeightTime < 1.0E-09 ) ? 0 : fTowerECalTime/fTowerECalWeightTime;
+
+  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerHCalEnergy);
+
+//  hcalEnergy = gRandom->Gaus(fTowerHCalEnergy, hcalSigma);
+//  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
+
+  hcalEnergy = LogNormal(fTowerHCalEnergy, hcalSigma);
+  hcalTime = (fTowerHCalWeightTime < 1.0E-09 ) ? 0 : fTowerHCalTime/fTowerHCalWeightTime;
+
+  energy = ecalEnergy + hcalEnergy;
+  time = (TMath::Sqrt(ecalEnergy)*ecalTime + TMath::Sqrt(hcalEnergy)*hcalTime)/(TMath::Sqrt(ecalEnergy) + TMath::Sqrt(hcalEnergy)); 
+
+//  eta = fTowerEta;
+//  phi = fTowerPhi;
 
   eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
@@ -392,5 +445,7 @@
   fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
   fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
- 
+  fTower->Eem = ecalEnergy;
+  fTower->Ehad = hcalEnergy;
+
   fTower->Edges[0] = fTowerEdges[0];
   fTower->Edges[1] = fTowerEdges[1];
@@ -400,21 +455,63 @@
 
   // fill calorimeter towers
-  if(energy > 0.0) fTowerOutputArray->Add(fTower);
+  if(energy > 0.0)
+  {
+    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
+    {
+      fPhotonOutputArray->Add(fTower);
+    }
+   
+    fTowerOutputArray->Add(fTower);
+  }
+
+  // fill energy flow candidates
+
+  // save all the tracks as energy flow tracks
+  fItTowerTrackArray->Reset();
+  while((track = static_cast<Candidate*>(fItTowerTrackArray->Next())))
+  {
+    fEFlowTrackOutputArray->Add(track);
+  }
+
+  ecalEnergy -= fTrackECalEnergy;
+  if(ecalEnergy < 0.0) ecalEnergy = 0.0;
+
+  hcalEnergy -= fTrackHCalEnergy;
+  if(hcalEnergy < 0.0) hcalEnergy = 0.0;
+
+  energy = ecalEnergy + hcalEnergy;
 
   
-  // fill energy flow candidates
-  energy -= fTrackEnergy;
-  if(energy < 0.0) energy = 0.0;
-   
-  // save energy excess as an energy flow tower
-  if(energy > 0.0)
+  // save ECAL and/or HCAL energy excess as an energy flow tower
+  if(ecalEnergy > 0.0)
   {
     // create new photon tower
     tower = static_cast<Candidate*>(fTower->Clone());
-    pt = energy / TMath::CosH(eta);
-
-    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
-    fEFlowTowerOutputArray->Add(tower);
-  }
+
+    pt = ecalEnergy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalEnergy);
+    tower->Eem = ecalEnergy;
+    tower->Ehad = 0;
+
+    fEFlowPhotonOutputArray->Add(tower);
+  }
+
+  if(hcalEnergy > 0.0)
+  {
+    // create new neutral hadron tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+
+    pt = hcalEnergy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalEnergy);
+    tower->Eem = 0;
+    tower->Ehad = hcalEnergy;
+
+    fEFlowNeutralHadronOutputArray->Add(tower);
+  }
+
+
+
 
 }
Index: modules/Calorimeter.h
===================================================================
--- modules/Calorimeter.h	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ modules/Calorimeter.h	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
@@ -38,17 +38,17 @@
 private:
 
-  typedef std::map< Long64_t, Double_t > TFractionMap; //!
+  typedef std::map< Long64_t, std::pair< Double_t, Double_t > > TFractionMap; //!
   typedef std::map< Double_t, std::set< Double_t > > TBinMap; //!
 
   Candidate *fTower;
   Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
-  Double_t fTowerEnergy;
-  Double_t fTrackEnergy;
+  Double_t fTowerECalEnergy, fTowerHCalEnergy;
+  Double_t fTrackECalEnergy, fTrackHCalEnergy;
   
-  Double_t fTowerTime;
-  Double_t fTrackTime;
+  Double_t fTowerECalTime, fTowerHCalTime;
+  Double_t fTrackECalTime, fTrackHCalTime;
    
-  Double_t fTowerWeightTime;
-  Double_t fTrackWeightTime;
+  Double_t fTowerECalWeightTime, fTowerHCalWeightTime;
+  Double_t fTrackECalWeightTime, fTrackHCalWeightTime;
   
   Int_t fTowerTrackHits, fTowerPhotonHits;
@@ -62,10 +62,13 @@
   std::vector < Long64_t > fTowerHits;
 
-  std::vector < Double_t > fTowerFractions;
-  
-  std::vector < Double_t > fTrackFractions;
- 
-  DelphesFormula *fResolutionFormula; //!
- 
+  std::vector < Double_t > fTowerECalFractions;
+  std::vector < Double_t > fTowerHCalFractions;
+
+  std::vector < Double_t > fTrackECalFractions;
+  std::vector < Double_t > fTrackHCalFractions;
+
+  DelphesFormula *fECalResolutionFormula; //!
+  DelphesFormula *fHCalResolutionFormula; //!
+
   TIterator *fItParticleInputArray; //!
   TIterator *fItTrackInputArray; //!
@@ -75,7 +78,10 @@
 
   TObjArray *fTowerOutputArray; //!
+  TObjArray *fPhotonOutputArray; //!
  
-  TObjArray *fEFlowTowerOutputArray; //!
- 
+  TObjArray *fEFlowTrackOutputArray; //!
+  TObjArray *fEFlowPhotonOutputArray; //!
+  TObjArray *fEFlowNeutralHadronOutputArray; //!
+
   TObjArray *fTowerTrackArray; //!
   TIterator *fItTowerTrackArray; //!
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision 9687203133c32a58144435a449f6d077a04b38c7)
+++ modules/ModulesLinkDef.h	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
@@ -21,4 +21,5 @@
 #include "modules/ImpactParameterSmearing.h"
 #include "modules/TimeSmearing.h"
+#include "modules/SimpleCalorimeter.h"
 #include "modules/Calorimeter.h"
 #include "modules/Isolation.h"
@@ -57,4 +58,5 @@
 #pragma link C++ class ImpactParameterSmearing+;
 #pragma link C++ class TimeSmearing+;
+#pragma link C++ class SimpleCalorimeter+;
 #pragma link C++ class Calorimeter+;
 #pragma link C++ class Isolation+;
Index: modules/SimpleCalorimeter.cc
===================================================================
--- modules/SimpleCalorimeter.cc	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
+++ modules/SimpleCalorimeter.cc	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
@@ -0,0 +1,439 @@
+
+/** \class SimpleCalorimeter
+ *
+ *  Fills SimpleCalorimeter towers, performs SimpleCalorimeter resolution smearing,
+ *  and creates energy flow objects (tracks, photons, and neutral hadrons).
+ *
+ *  $Date: 2014-04-16 15:29:31 +0200 (Wed, 16 Apr 2014) $
+ *  $Revision: 1364 $
+ *
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/SimpleCalorimeter.h"
+
+#include "classes/DelphesClasses.h"
+#include "classes/DelphesFactory.h"
+#include "classes/DelphesFormula.h"
+
+#include "ExRootAnalysis/ExRootResult.h"
+#include "ExRootAnalysis/ExRootFilter.h"
+#include "ExRootAnalysis/ExRootClassifier.h"
+
+#include "TMath.h"
+#include "TString.h"
+#include "TFormula.h"
+#include "TRandom3.h"
+#include "TObjArray.h"
+#include "TDatabasePDG.h"
+#include "TLorentzVector.h"
+
+#include <algorithm>
+#include <stdexcept>
+#include <iostream>
+#include <sstream>
+
+using namespace std;
+
+//------------------------------------------------------------------------------
+
+SimpleCalorimeter::SimpleCalorimeter() :
+  fResolutionFormula(0),
+  fItParticleInputArray(0), fItTrackInputArray(0),
+  fTowerTrackArray(0), fItTowerTrackArray(0)
+{
+  fResolutionFormula = new DelphesFormula;
+ 
+  fTowerTrackArray = new TObjArray;
+  fItTowerTrackArray = fTowerTrackArray->MakeIterator();
+}
+
+//------------------------------------------------------------------------------
+
+SimpleCalorimeter::~SimpleCalorimeter()
+{
+  if(fResolutionFormula) delete fResolutionFormula;
+ 
+  if(fTowerTrackArray) delete fTowerTrackArray;
+  if(fItTowerTrackArray) delete fItTowerTrackArray;
+}
+
+//------------------------------------------------------------------------------
+
+void SimpleCalorimeter::Init()
+{
+  ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
+  Double_t fraction;
+  TBinMap::iterator itEtaBin;
+  set< Double_t >::iterator itPhiBin;
+  vector< Double_t > *phiBins;
+
+  // read eta and phi bins
+  param = GetParam("EtaPhiBins");
+  size = param.GetSize();
+  fBinMap.clear();
+  fEtaBins.clear();
+  fPhiBins.clear();
+  for(i = 0; i < size/2; ++i)
+  {
+    paramEtaBins = param[i*2];
+    sizeEtaBins = paramEtaBins.GetSize();
+    paramPhiBins = param[i*2 + 1];
+    sizePhiBins = paramPhiBins.GetSize();
+
+    for(j = 0; j < sizeEtaBins; ++j)
+    {
+      for(k = 0; k < sizePhiBins; ++k)
+      {
+        fBinMap[paramEtaBins[j].GetDouble()].insert(paramPhiBins[k].GetDouble());
+      }
+    }
+  }
+
+  // for better performance we transform map of sets to parallel vectors:
+  // vector< double > and vector< vector< double >* >
+  for(itEtaBin = fBinMap.begin(); itEtaBin != fBinMap.end(); ++itEtaBin)
+  {
+    fEtaBins.push_back(itEtaBin->first);
+    phiBins = new vector< double >(itEtaBin->second.size());
+    fPhiBins.push_back(phiBins);
+    phiBins->clear();
+    for(itPhiBin = itEtaBin->second.begin(); itPhiBin != itEtaBin->second.end(); ++itPhiBin)
+    {
+      phiBins->push_back(*itPhiBin);
+    }
+  }
+
+  // read energy fractions for different particles
+  param = GetParam("EnergyFraction");
+  size = param.GetSize();
+
+  // set default energy fractions values
+  fFractionMap.clear();
+  fFractionMap[0] = 1.0;
+
+  for(i = 0; i < size/2; ++i)
+  {
+    paramFractions = param[i*2 + 1];
+    fraction = paramFractions[0].GetDouble();
+    fFractionMap[param[i*2].GetInt()] = fraction;
+  }
+/*
+  TFractionMap::iterator itFractionMap;
+  for(itFractionMap = fFractionMap.begin(); itFractionMap != fFractionMap.end(); ++itFractionMap)
+  {
+    cout << itFractionMap->first << "   " << itFractionMap->second.first  << "   " << itFractionMap->second.second << endl;
+  }
+*/
+  // read resolution formulas
+  fResolutionFormula->Compile(GetString("ResolutionFormula", "0"));
+ 
+  // import array with output from other modules
+  fParticleInputArray = ImportArray(GetString("ParticleInputArray", "ParticlePropagator/particles"));
+  fItParticleInputArray = fParticleInputArray->MakeIterator();
+
+  fTrackInputArray = ImportArray(GetString("TrackInputArray", "ParticlePropagator/tracks"));
+  fItTrackInputArray = fTrackInputArray->MakeIterator();
+
+  // create output arrays
+  fTowerOutputArray = ExportArray(GetString("TowerOutputArray", "towers"));
+  fEFlowTowerOutputArray = ExportArray(GetString("EFlowTowerOutputArray", "eflowTowers"));
+ 
+}
+
+//------------------------------------------------------------------------------
+
+void SimpleCalorimeter::Finish()
+{
+  vector< vector< Double_t >* >::iterator itPhiBin;
+  if(fItParticleInputArray) delete fItParticleInputArray;
+  if(fItTrackInputArray) delete fItTrackInputArray;
+  for(itPhiBin = fPhiBins.begin(); itPhiBin != fPhiBins.end(); ++itPhiBin)
+  {
+    delete *itPhiBin;
+  }
+}
+
+//------------------------------------------------------------------------------
+
+void SimpleCalorimeter::Process()
+{
+  Candidate *particle, *track;
+  TLorentzVector position, momentum;
+  Short_t etaBin, phiBin, flags;
+  Int_t number;
+  Long64_t towerHit, towerEtaPhi, hitEtaPhi;
+  Double_t fraction;
+  Double_t energy;
+  Int_t pdgCode;
+
+  TFractionMap::iterator itFractionMap;
+
+  vector< Double_t >::iterator itEtaBin;
+  vector< Double_t >::iterator itPhiBin;
+  vector< Double_t > *phiBins;
+
+  vector< Long64_t >::iterator itTowerHits;
+
+  DelphesFactory *factory = GetFactory();
+  fTowerHits.clear();
+  fTowerFractions.clear();
+  fTrackFractions.clear();
+ 
+  // loop over all particles
+  fItParticleInputArray->Reset();
+  number = -1;
+  while((particle = static_cast<Candidate*>(fItParticleInputArray->Next())))
+  {
+    const TLorentzVector &particlePosition = particle->Position;
+    ++number;
+
+    pdgCode = TMath::Abs(particle->PID);
+
+    itFractionMap = fFractionMap.find(pdgCode);
+    if(itFractionMap == fFractionMap.end())
+    {
+      itFractionMap = fFractionMap.find(0);
+    }
+
+    fraction = itFractionMap->second;
+    fTowerFractions.push_back(fraction);
+   
+    if(fraction < 1.0E-9) continue;
+
+    // find eta bin [1, fEtaBins.size - 1]
+    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), particlePosition.Eta());
+    if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
+    etaBin = distance(fEtaBins.begin(), itEtaBin);
+
+    // phi bins for given eta bin
+    phiBins = fPhiBins[etaBin];
+
+    // find phi bin [1, phiBins.size - 1]
+    itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), particlePosition.Phi());
+    if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
+    phiBin = distance(phiBins->begin(), itPhiBin);
+
+    flags = 0;
+    flags |= (pdgCode == 11 || pdgCode == 22) << 1;
+
+    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for particle number}
+    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
+
+    fTowerHits.push_back(towerHit);
+  }
+
+  // loop over all tracks
+  fItTrackInputArray->Reset();
+  number = -1;
+  while((track = static_cast<Candidate*>(fItTrackInputArray->Next())))
+  {
+    const TLorentzVector &trackPosition = track->Position;
+    ++number;
+
+    pdgCode = TMath::Abs(track->PID);
+
+    itFractionMap = fFractionMap.find(pdgCode);
+    if(itFractionMap == fFractionMap.end())
+    {
+      itFractionMap = fFractionMap.find(0);
+    }
+
+    fraction = itFractionMap->second;
+  
+    fTrackFractions.push_back(fraction);
+  
+    // find eta bin [1, fEtaBins.size - 1]
+    itEtaBin = lower_bound(fEtaBins.begin(), fEtaBins.end(), trackPosition.Eta());
+    if(itEtaBin == fEtaBins.begin() || itEtaBin == fEtaBins.end()) continue;
+    etaBin = distance(fEtaBins.begin(), itEtaBin);
+
+    // phi bins for given eta bin
+    phiBins = fPhiBins[etaBin];
+
+    // find phi bin [1, phiBins.size - 1]
+    itPhiBin = lower_bound(phiBins->begin(), phiBins->end(), trackPosition.Phi());
+    if(itPhiBin == phiBins->begin() || itPhiBin == phiBins->end()) continue;
+    phiBin = distance(phiBins->begin(), itPhiBin);
+
+    flags = 1;
+
+    // make tower hit {16-bits for eta bin number, 16-bits for phi bin number, 8-bits for flags, 24-bits for track number}
+    towerHit = (Long64_t(etaBin) << 48) | (Long64_t(phiBin) << 32) | (Long64_t(flags) << 24) | Long64_t(number);
+
+    fTowerHits.push_back(towerHit);
+  }
+
+  // all hits are sorted first by eta bin number, then by phi bin number,
+  // then by flags and then by particle or track number
+  sort(fTowerHits.begin(), fTowerHits.end());
+
+  // loop over all hits
+  towerEtaPhi = 0;
+  fTower = 0;
+  for(itTowerHits = fTowerHits.begin(); itTowerHits != fTowerHits.end(); ++itTowerHits)
+  {
+    towerHit = (*itTowerHits);
+    flags = (towerHit >> 24) & 0x00000000000000FFLL;
+    number = (towerHit) & 0x0000000000FFFFFFLL;
+    hitEtaPhi = towerHit >> 32;
+
+    if(towerEtaPhi != hitEtaPhi)
+    {
+      // switch to next tower
+      towerEtaPhi = hitEtaPhi;
+
+      // finalize previous tower
+      FinalizeTower();
+
+      // create new tower
+      fTower = factory->NewCandidate();
+
+      phiBin = (towerHit >> 32) & 0x000000000000FFFFLL;
+      etaBin = (towerHit >> 48) & 0x000000000000FFFFLL;
+
+      // phi bins for given eta bin
+      phiBins = fPhiBins[etaBin];
+
+      // calculate eta and phi of the tower's center
+      fTowerEta = 0.5*(fEtaBins[etaBin - 1] + fEtaBins[etaBin]);
+      fTowerPhi = 0.5*((*phiBins)[phiBin - 1] + (*phiBins)[phiBin]);
+
+      fTowerEdges[0] = fEtaBins[etaBin - 1];
+      fTowerEdges[1] = fEtaBins[etaBin];
+      fTowerEdges[2] = (*phiBins)[phiBin - 1];
+      fTowerEdges[3] = (*phiBins)[phiBin];
+
+      fTowerEnergy = 0.0;
+      fTrackEnergy = 0.0;
+      
+      fTowerTime = 0.0;
+      fTrackTime = 0.0;
+     
+      fTowerWeightTime = 0.0; 
+      
+      fTowerTrackHits = 0;
+      fTowerPhotonHits = 0;
+      
+      fTowerTrackArray->Clear();
+    }
+
+    // check for track hits
+    if(flags & 1)
+    {
+      ++fTowerTrackHits;
+
+      track = static_cast<Candidate*>(fTrackInputArray->At(number));
+      momentum = track->Momentum;
+      position = track->Position;
+ 
+      energy = momentum.E() * fTrackFractions[number];
+      
+      fTrackEnergy += energy;
+      
+      fTrackTime += TMath::Sqrt(energy)*position.T();
+      fTrackWeightTime += TMath::Sqrt(energy);
+    
+      fTowerTrackArray->Add(track);
+
+      continue;
+    }
+   
+    // check for photon and electron hits in current tower
+    if(flags & 2) ++fTowerPhotonHits;
+    
+    particle = static_cast<Candidate*>(fParticleInputArray->At(number));
+    momentum = particle->Momentum;
+    position = particle->Position;
+
+    // fill current tower
+    energy = momentum.E() * fTowerFractions[number];
+   
+    fTowerEnergy += energy;
+    
+    fTowerTime += TMath::Sqrt(energy)*position.T();
+    fTowerWeightTime += TMath::Sqrt(energy);
+    
+    fTower->AddCandidate(particle);
+  }
+
+  // finalize last tower
+  FinalizeTower();
+}
+
+//------------------------------------------------------------------------------
+
+void SimpleCalorimeter::FinalizeTower()
+{
+  Candidate *tower;
+  Double_t energy, pt, eta, phi;
+  Double_t sigma;
+  Double_t time;
+
+  if(!fTower) return;
+
+  sigma = fResolutionFormula->Eval(0.0, fTowerEta, 0.0, fTowerEnergy);
+
+//  energy = gRandom->Gaus(fTowerEnergy, sigma);
+//  if(energy < 0.0) energy = 0.0;
+
+  energy = LogNormal(fTowerEnergy, sigma);
+  time = (fTowerWeightTime < 1.0E-09 ) ? 0 : fTowerTime/fTowerWeightTime;
+
+  eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
+  phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
+
+  pt = energy / TMath::CosH(eta);
+
+ // fTower->Position.SetXYZT(-time, 0.0, 0.0, time);
+  fTower->Position.SetPtEtaPhiE(1.0, eta, phi, time);
+  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+ 
+  fTower->Edges[0] = fTowerEdges[0];
+  fTower->Edges[1] = fTowerEdges[1];
+  fTower->Edges[2] = fTowerEdges[2];
+  fTower->Edges[3] = fTowerEdges[3];
+
+
+  // fill SimpleCalorimeter towers
+  if(energy > 0.0) fTowerOutputArray->Add(fTower);
+
+  
+  // fill energy flow candidates
+  energy -= fTrackEnergy;
+  if(energy < 0.0) energy = 0.0;
+   
+  // save energy excess as an energy flow tower
+  if(energy > 0.0)
+  {
+    // create new photon tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt = energy / TMath::CosH(eta);
+
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+    fEFlowTowerOutputArray->Add(tower);
+  }
+
+}
+
+//------------------------------------------------------------------------------
+
+Double_t SimpleCalorimeter::LogNormal(Double_t mean, Double_t sigma)
+{
+  Double_t a, b;
+
+  if(mean > 0.0)
+  {
+    b = TMath::Sqrt(TMath::Log((1.0 + (sigma*sigma)/(mean*mean))));
+    a = TMath::Log(mean) - 0.5*b*b;
+
+    return TMath::Exp(a + b*gRandom->Gaus(0, 1));
+  }
+  else
+  {
+    return 0.0;
+  }
+}
Index: modules/SimpleCalorimeter.h
===================================================================
--- modules/SimpleCalorimeter.h	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
+++ modules/SimpleCalorimeter.h	(revision 2e6a81baac3061592fe8970d1f5767e88a116640)
@@ -0,0 +1,89 @@
+#ifndef SimpleCalorimeter_h
+#define SimpleCalorimeter_h
+
+/** \class SimpleCalorimeter
+ *
+ *  Fills SimpleCalorimeter towers, performs SimpleCalorimeter resolution smearing,
+ *  and creates energy flow objects (tracks, photons, and neutral hadrons).
+ *
+ *  $Date: 2014-04-16 15:29:31 +0200 (Wed, 16 Apr 2014) $
+ *  $Revision: 1364 $
+ *
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+#include <map>
+#include <set>
+#include <vector>
+
+class TObjArray;
+class DelphesFormula;
+class Candidate;
+
+class SimpleCalorimeter: public DelphesModule
+{
+public:
+
+  SimpleCalorimeter();
+  ~SimpleCalorimeter();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  typedef std::map< Long64_t, Double_t > TFractionMap; //!
+  typedef std::map< Double_t, std::set< Double_t > > TBinMap; //!
+
+  Candidate *fTower;
+  Double_t fTowerEta, fTowerPhi, fTowerEdges[4];
+  Double_t fTowerEnergy;
+  Double_t fTrackEnergy;
+  
+  Double_t fTowerTime;
+  Double_t fTrackTime;
+   
+  Double_t fTowerWeightTime;
+  Double_t fTrackWeightTime;
+  
+  Int_t fTowerTrackHits, fTowerPhotonHits;
+
+  TFractionMap fFractionMap; //!
+  TBinMap fBinMap; //!
+
+  std::vector < Double_t > fEtaBins;
+  std::vector < std::vector < Double_t >* > fPhiBins;
+
+  std::vector < Long64_t > fTowerHits;
+
+  std::vector < Double_t > fTowerFractions;
+  
+  std::vector < Double_t > fTrackFractions;
+ 
+  DelphesFormula *fResolutionFormula; //!
+ 
+  TIterator *fItParticleInputArray; //!
+  TIterator *fItTrackInputArray; //!
+
+  const TObjArray *fParticleInputArray; //!
+  const TObjArray *fTrackInputArray; //!
+
+  TObjArray *fTowerOutputArray; //!
+ 
+  TObjArray *fEFlowTowerOutputArray; //!
+ 
+  TObjArray *fTowerTrackArray; //!
+  TIterator *fItTowerTrackArray; //!
+
+  void FinalizeTower();
+  Double_t LogNormal(Double_t mean, Double_t sigma);
+
+  ClassDef(SimpleCalorimeter, 1)
+};
+
+#endif
