Index: cards/delphes_card_IDEA.tcl
===================================================================
--- cards/delphes_card_IDEA.tcl	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
+++ cards/delphes_card_IDEA.tcl	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
@@ -0,0 +1,621 @@
+set MaxEvents 10
+
+#######################################
+# Order of execution of various modules
+#######################################
+
+set ExecutionPath {
+  ParticlePropagator
+
+  ChargedHadronTrackingEfficiency
+  ElectronTrackingEfficiency
+  MuonTrackingEfficiency
+
+  ChargedHadronMomentumSmearing
+  ElectronMomentumSmearing
+  MuonMomentumSmearing
+
+  TrackMerger
+  Calorimeter
+  EFlowMerger
+
+  PhotonEfficiency
+  PhotonIsolation
+
+  ElectronFilter
+  ElectronEfficiency
+  ElectronIsolation
+
+  MuonEfficiency
+  MuonIsolation
+
+  MissingET
+
+  NeutrinoFilter
+  GenJetFinder
+  GenMissingET
+  
+  FastJetFinder
+
+  JetEnergyScale
+
+  JetFlavorAssociation
+
+  BTagging
+  TauTagging
+
+  UniqueObjectFinder
+
+  ScalarHT
+
+  TreeWriter
+}
+
+#################################
+# Propagate particles in cylinder
+#################################
+
+module ParticlePropagator ParticlePropagator {
+  set InputArray Delphes/stableParticles
+
+  set OutputArray stableParticles
+  set ChargedHadronOutputArray chargedHadrons
+  set ElectronOutputArray electrons
+  set MuonOutputArray muons
+
+  # radius of the magnetic field coverage, in m
+  set Radius 1.29
+  # half-length of the magnetic field coverage, in m
+  set HalfLength 3.00
+
+  # magnetic field
+  set Bz 3.8
+}
+
+####################################
+# Charged hadron tracking efficiency
+####################################
+
+module Efficiency ChargedHadronTrackingEfficiency {
+  set InputArray ParticlePropagator/chargedHadrons
+  set OutputArray chargedHadrons
+
+  # add EfficiencyFormula {efficiency formula as a function of eta and pt}
+
+  # tracking efficiency formula for charged hadrons
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.60) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.85) +
+                         (abs(eta) > 2.5)                                                  * (0.00)}
+}
+
+##############################
+# Electron tracking efficiency
+##############################
+
+module Efficiency ElectronTrackingEfficiency {
+  set InputArray ParticlePropagator/electrons
+  set OutputArray electrons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+
+  # tracking efficiency formula for electrons
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.73) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0   && pt <= 1.0e2) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e2)                * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.50) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0   && pt <= 1.0e2) * (0.83) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0e2)                * (0.90) +
+                         (abs(eta) > 2.5)                                                  * (0.00)}
+}
+
+##########################
+# Muon tracking efficiency
+##########################
+
+module Efficiency MuonTrackingEfficiency {
+  set InputArray ParticlePropagator/muons
+  set OutputArray muons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+
+  # tracking efficiency formula for muons
+  set EfficiencyFormula {                                                    (pt <= 0.1)   * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 0.1   && pt <= 1.0)   * (0.75) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0)                  * (0.99) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1   && pt <= 1.0)   * (0.70) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 1.0)                  * (0.98) +
+                         (abs(eta) > 2.5)                                                  * (0.00)}
+}
+
+########################################
+# Momentum resolution for charged tracks
+########################################
+
+module MomentumSmearing ChargedHadronMomentumSmearing {
+  set InputArray ChargedHadronTrackingEfficiency/chargedHadrons
+  set OutputArray chargedHadrons
+
+  # set ResolutionFormula {resolution formula as a function of eta and pt}
+
+  # resolution formula for charged hadrons
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*2.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*5.5e-4^2)}
+}
+
+###################################
+# Momentum resolution for electrons
+###################################
+
+module MomentumSmearing ElectronMomentumSmearing {
+  set InputArray ElectronTrackingEfficiency/electrons
+  set OutputArray electrons
+
+  # set ResolutionFormula {resolution formula as a function of eta and energy}
+
+  # resolution formula for electrons
+  # based on arXiv:1405.6569
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.03^2 + pt^2*1.3e-3^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.05^2 + pt^2*1.7e-3^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.15^2 + pt^2*3.1e-3^2)}
+}
+
+###############################
+# Momentum resolution for muons
+###############################
+
+module MomentumSmearing MuonMomentumSmearing {
+  set InputArray MuonTrackingEfficiency/muons
+  set OutputArray muons
+
+  # set ResolutionFormula {resolution formula as a function of eta and pt}
+
+  # resolution formula for muons
+  set ResolutionFormula {                  (abs(eta) <= 0.5) * (pt > 0.1) * sqrt(0.01^2 + pt^2*1.0e-4^2) +
+                         (abs(eta) > 0.5 && abs(eta) <= 1.5) * (pt > 0.1) * sqrt(0.015^2 + pt^2*1.5e-4^2) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 0.1) * sqrt(0.025^2 + pt^2*3.5e-4^2)}
+}
+
+##############
+# Track merger
+##############
+
+module Merger TrackMerger {
+# add InputArray InputArray
+  add InputArray ChargedHadronMomentumSmearing/chargedHadrons
+  add InputArray ElectronMomentumSmearing/electrons
+  add InputArray MuonMomentumSmearing/muons
+  set OutputArray tracks
+}
+
+#############
+# Calorimeter
+#############
+
+module DualReadoutCalorimeter Calorimeter {
+  set ParticleInputArray ParticlePropagator/stableParticles
+  set TrackInputArray TrackMerger/tracks
+
+  set TowerOutputArray towers
+  set PhotonOutputArray photons
+
+  set EFlowTrackOutputArray eflowTracks
+  set EFlowPhotonOutputArray eflowPhotons
+  set EFlowNeutralHadronOutputArray eflowNeutralHadrons
+
+  set ECalEnergyMin 0.5
+  set HCalEnergyMin 1.0
+  set EnergyMin 0.0
+
+  set ECalEnergySignificanceMin 1.0
+  set HCalEnergySignificanceMin 1.0
+  set EnergySignificanceMin 1.0
+
+  set SmearTowerCenter true
+
+  set pi [expr {acos(-1)}]
+
+  # lists of the edges of each tower in eta and phi
+  # each list starts with the lower edge of the first tower
+  # the list ends with the higher edged of the last tower
+
+  # 5 degrees towers
+  set PhiBins {}
+  for {set i -36} {$i <= 36} {incr i} {
+    add PhiBins [expr {$i * $pi/36.0}]
+  }
+  foreach eta {-1.566 -1.479 -1.392 -1.305 -1.218 -1.131 -1.044 -0.957 -0.87 -0.783 -0.696 -0.609 -0.522 -0.435 -0.348 -0.261 -0.174 -0.087 0 0.087 0.174 0.261 0.348 0.435 0.522 0.609 0.696 0.783 0.87 0.957 1.044 1.131 1.218 1.305 1.392 1.479 1.566 1.653} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # 10 degrees towers
+  set PhiBins {}
+  for {set i -18} {$i <= 18} {incr i} {
+    add PhiBins [expr {$i * $pi/18.0}]
+  }
+  foreach eta {-4.35 -4.175 -4 -3.825 -3.65 -3.475 -3.3 -3.125 -2.95 -2.868 -2.65 -2.5 -2.322 -2.172 -2.043 -1.93 -1.83 -1.74 -1.653 1.74 1.83 1.93 2.043 2.172 2.322 2.5 2.65 2.868 2.95 3.125 3.3 3.475 3.65 3.825 4 4.175 4.35 4.525} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # 20 degrees towers
+  set PhiBins {}
+  for {set i -9} {$i <= 9} {incr i} {
+    add PhiBins [expr {$i * $pi/9.0}]
+  }
+  foreach eta {-5 -4.7 -4.525 4.7 5} {
+    add EtaPhiBins $eta $PhiBins
+  }
+
+  # default energy fractions {abs(PDG code)} {Fecal Fhcal}
+  add EnergyFraction {0} {0.0 1.0}
+  # energy fractions for e, gamma and pi0
+  add EnergyFraction {11} {1.0 0.0}
+  add EnergyFraction {22} {1.0 0.0}
+  add EnergyFraction {111} {1.0 0.0}
+  # energy fractions for muon, neutrinos and neutralinos
+  add EnergyFraction {12} {0.0 0.0}
+  add EnergyFraction {13} {0.0 0.0}
+  add EnergyFraction {14} {0.0 0.0}
+  add EnergyFraction {16} {0.0 0.0}
+  add EnergyFraction {1000022} {0.0 0.0}
+  add EnergyFraction {1000023} {0.0 0.0}
+  add EnergyFraction {1000025} {0.0 0.0}
+  add EnergyFraction {1000035} {0.0 0.0}
+  add EnergyFraction {1000045} {0.0 0.0}
+  # energy fractions for K0short and Lambda
+  add EnergyFraction {310} {0.3 0.7}
+  add EnergyFraction {3122} {0.3 0.7}
+
+  # set ECalResolutionFormula {resolution formula as a function of eta and energy}
+  # Eta shape from arXiv:1306.2016, Energy shape from arXiv:1502.02701
+  set ECalResolutionFormula {                  (abs(eta) <= 1.5) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 1.5 && abs(eta) <= 2.5) * sqrt(energy^2*0.008^2 + energy*0.11^2 + 0.40^2) +
+                             (abs(eta) > 2.5 && abs(eta) <= 5.0) * sqrt(energy^2*0.107^2 + energy*2.08^2)}
+
+  # set HCalResolutionFormula {resolution formula as a function of eta and energy}
+  set HCalResolutionFormula {                  (abs(eta) <= 3.0) * sqrt(energy^2*0.050^2 + energy*1.50^2) +
+                             (abs(eta) > 3.0 && abs(eta) <= 5.0) * sqrt(energy^2*0.130^2 + energy*2.70^2)}
+}
+
+####################
+# Energy flow merger
+####################
+
+module Merger EFlowMerger {
+# add InputArray InputArray
+  add InputArray Calorimeter/eflowTracks
+  add InputArray Calorimeter/eflowPhotons
+  add InputArray Calorimeter/eflowNeutralHadrons
+  set OutputArray eflow
+}
+
+###################
+# Photon efficiency
+###################
+
+module Efficiency PhotonEfficiency {
+  set InputArray Calorimeter/eflowPhotons
+  set OutputArray photons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+
+  # efficiency formula for photons
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
+                         (abs(eta) > 2.5)                                   * (0.00)}
+}
+
+##################
+# Photon isolation
+##################
+
+module Isolation PhotonIsolation {
+  set CandidateInputArray PhotonEfficiency/photons
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray photons
+
+  set DeltaRMax 0.5
+
+  set PTMin 0.5
+
+  set PTRatioMax 0.12
+}
+
+#################
+# Electron filter
+#################
+
+module PdgCodeFilter ElectronFilter {
+  set InputArray Calorimeter/eflowTracks
+  set OutputArray electrons
+  set Invert true
+  add PdgCode {11}
+  add PdgCode {-11}
+}
+
+#####################
+# Electron efficiency
+#####################
+
+module Efficiency ElectronEfficiency {
+  set InputArray ElectronFilter/electrons
+  set OutputArray electrons
+
+  # set EfficiencyFormula {efficiency formula as a function of eta and pt}
+
+  # efficiency formula for electrons
+  set EfficiencyFormula {                                      (pt <= 10.0) * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0)  * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.5) * (pt > 10.0)  * (0.85) +
+                         (abs(eta) > 2.5)                                   * (0.00)}
+}
+
+####################
+# Electron isolation
+####################
+
+module Isolation ElectronIsolation {
+  set CandidateInputArray ElectronEfficiency/electrons
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray electrons
+
+  set DeltaRMax 0.5
+
+  set PTMin 0.5
+
+  set PTRatioMax 0.12
+}
+
+#################
+# Muon efficiency
+#################
+
+module Efficiency MuonEfficiency {
+  set InputArray MuonMomentumSmearing/muons
+  set OutputArray muons
+
+  # set EfficiencyFormula {efficiency as a function of eta and pt}
+
+  # efficiency formula for muons
+  set EfficiencyFormula {                                      (pt <= 10.0)               * (0.00) +
+                                           (abs(eta) <= 1.5) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                                           (abs(eta) <= 1.5) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 10.0 && pt <= 1.0e3) * (0.95) +
+                         (abs(eta) > 1.5 && abs(eta) <= 2.4) * (pt > 1.0e3)               * (0.95 * exp(0.5 - pt*5.0e-4)) +
+                         (abs(eta) > 2.4)                                                 * (0.00)}
+}
+
+################
+# Muon isolation
+################
+
+module Isolation MuonIsolation {
+  set CandidateInputArray MuonEfficiency/muons
+  set IsolationInputArray EFlowMerger/eflow
+
+  set OutputArray muons
+
+  set DeltaRMax 0.5
+
+  set PTMin 0.5
+
+  set PTRatioMax 0.25
+}
+
+###################
+# Missing ET merger
+###################
+
+module Merger MissingET {
+# add InputArray InputArray
+  add InputArray EFlowMerger/eflow
+  set MomentumOutputArray momentum
+}
+
+##################
+# Scalar HT merger
+##################
+
+module Merger ScalarHT {
+# add InputArray InputArray
+  add InputArray UniqueObjectFinder/jets
+  add InputArray UniqueObjectFinder/electrons
+  add InputArray UniqueObjectFinder/photons
+  add InputArray UniqueObjectFinder/muons
+  set EnergyOutputArray energy
+}
+
+
+#####################
+# Neutrino Filter
+#####################
+
+module PdgCodeFilter NeutrinoFilter {
+
+  set InputArray Delphes/stableParticles
+  set OutputArray filteredParticles
+
+  set PTMin 0.0
+
+  add PdgCode {12}
+  add PdgCode {14}
+  add PdgCode {16}
+  add PdgCode {-12}
+  add PdgCode {-14}
+  add PdgCode {-16}
+
+}
+
+
+#####################
+# MC truth jet finder
+#####################
+
+module FastJetFinder GenJetFinder {
+  set InputArray NeutrinoFilter/filteredParticles
+
+  set OutputArray jets
+
+  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 6
+  set ParameterR 0.5
+
+  set JetPTMin 20.0
+}
+
+#########################
+# Gen Missing ET merger
+########################
+
+module Merger GenMissingET {
+# add InputArray InputArray
+  add InputArray NeutrinoFilter/filteredParticles
+  set MomentumOutputArray momentum
+}
+
+
+
+############
+# Jet finder
+############
+
+module FastJetFinder FastJetFinder {
+#  set InputArray Calorimeter/towers
+  set InputArray EFlowMerger/eflow
+
+  set OutputArray jets
+
+  # algorithm: 1 CDFJetClu, 2 MidPoint, 3 SIScone, 4 kt, 5 Cambridge/Aachen, 6 antikt
+  set JetAlgorithm 6
+  set ParameterR 0.5
+
+  set JetPTMin 20.0
+}
+
+##################
+# Jet Energy Scale
+##################
+
+module EnergyScale JetEnergyScale {
+  set InputArray FastJetFinder/jets
+  set OutputArray jets
+
+  # scale formula for jets
+  set ScaleFormula {sqrt( (2.5 - 0.15*(abs(eta)))^2 / pt + 1.0 )}
+}
+
+########################
+# Jet Flavor Association
+########################
+
+module JetFlavorAssociation JetFlavorAssociation {
+
+  set PartonInputArray Delphes/partons
+  set ParticleInputArray Delphes/allParticles
+  set ParticleLHEFInputArray Delphes/allParticlesLHEF
+  set JetInputArray JetEnergyScale/jets
+
+  set DeltaR 0.5
+  set PartonPTMin 1.0
+  set PartonEtaMax 2.5
+
+}
+
+###########
+# b-tagging
+###########
+
+module BTagging BTagging {
+  set JetInputArray JetEnergyScale/jets
+
+  set BitNumber 0
+
+  # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
+  # PDG code = the highest PDG code of a quark or gluon inside DeltaR cone around jet axis
+  # gluon's PDG code has the lowest priority
+
+  # based on arXiv:1211.4462
+  
+  # default efficiency formula (misidentification rate)
+  add EfficiencyFormula {0} {0.01+0.000038*pt}
+
+  # efficiency formula for c-jets (misidentification rate)
+  add EfficiencyFormula {4} {0.25*tanh(0.018*pt)*(1/(1+ 0.0013*pt))}
+
+  # efficiency formula for b-jets
+  add EfficiencyFormula {5} {0.85*tanh(0.0025*pt)*(25.0/(1+0.063*pt))}
+}
+
+#############
+# tau-tagging
+#############
+
+module TauTagging TauTagging {
+  set ParticleInputArray Delphes/allParticles
+  set PartonInputArray Delphes/partons
+  set JetInputArray JetEnergyScale/jets
+
+  set DeltaR 0.5
+
+  set TauPTMin 1.0
+
+  set TauEtaMax 2.5
+
+  # add EfficiencyFormula {abs(PDG code)} {efficiency formula as a function of eta and pt}
+
+  # default efficiency formula (misidentification rate)
+  add EfficiencyFormula {0} {0.01}
+  # efficiency formula for tau-jets
+  add EfficiencyFormula {15} {0.6}
+}
+
+#####################################################
+# Find uniquely identified photons/electrons/tau/jets
+#####################################################
+
+module UniqueObjectFinder UniqueObjectFinder {
+# earlier arrays take precedence over later ones
+# add InputArray InputArray OutputArray
+  add InputArray PhotonIsolation/photons photons
+  add InputArray ElectronIsolation/electrons electrons
+  add InputArray MuonIsolation/muons muons
+  add InputArray JetEnergyScale/jets jets
+}
+
+##################
+# ROOT tree writer
+##################
+
+# tracks, towers and eflow objects are not stored by default in the output.
+# if needed (for jet constituent or other studies), uncomment the relevant
+# "add Branch ..." lines.
+
+module TreeWriter TreeWriter {
+# add Branch InputArray BranchName BranchClass
+  add Branch Delphes/allParticles Particle GenParticle
+
+  add Branch TrackMerger/tracks Track Track
+  add Branch Calorimeter/towers Tower Tower
+
+  add Branch Calorimeter/eflowTracks EFlowTrack Track
+  add Branch Calorimeter/eflowPhotons EFlowPhoton Tower
+  add Branch Calorimeter/eflowNeutralHadrons EFlowNeutralHadron Tower
+
+  add Branch GenJetFinder/jets GenJet Jet
+  add Branch GenMissingET/momentum GenMissingET MissingET
+
+  add Branch UniqueObjectFinder/jets Jet Jet
+  add Branch UniqueObjectFinder/electrons Electron Electron
+  add Branch UniqueObjectFinder/photons Photon Photon
+  add Branch UniqueObjectFinder/muons Muon Muon
+  add Branch MissingET/momentum MissingET MissingET
+  add Branch ScalarHT/energy ScalarHT ScalarHT
+}
Index: modules/DualReadoutCalorimeter.cc
===================================================================
--- modules/DualReadoutCalorimeter.cc	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
+++ modules/DualReadoutCalorimeter.cc	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
@@ -0,0 +1,716 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+
+/** \class DualReadoutCalorimeter
+ *
+ *  Fills DualReadoutCalorimeter towers, performs DualReadoutCalorimeter resolution smearing,
+ *  and creates energy flow objects (tracks, photons, and neutral hadrons).
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "modules/DualReadoutCalorimeter.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;
+
+//------------------------------------------------------------------------------
+
+DualReadoutCalorimeter::DualReadoutCalorimeter() :
+  fECalResolutionFormula(0), fHCalResolutionFormula(0),
+  fItParticleInputArray(0), fItTrackInputArray(0)
+{
+  
+  fECalResolutionFormula = new DelphesFormula;
+  fHCalResolutionFormula = new DelphesFormula;
+
+  fECalTowerTrackArray = new TObjArray;
+  fItECalTowerTrackArray = fECalTowerTrackArray->MakeIterator();
+
+  fHCalTowerTrackArray = new TObjArray;
+  fItHCalTowerTrackArray = fHCalTowerTrackArray->MakeIterator();
+  
+}
+
+//------------------------------------------------------------------------------
+
+DualReadoutCalorimeter::~DualReadoutCalorimeter()
+{
+  
+  if(fECalResolutionFormula) delete fECalResolutionFormula;
+  if(fHCalResolutionFormula) delete fHCalResolutionFormula;
+
+  if(fECalTowerTrackArray) delete fECalTowerTrackArray;
+  if(fItECalTowerTrackArray) delete fItECalTowerTrackArray;
+
+  if(fHCalTowerTrackArray) delete fHCalTowerTrackArray;
+  if(fItHCalTowerTrackArray) delete fItHCalTowerTrackArray;
+ 
+}
+
+//------------------------------------------------------------------------------
+
+void DualReadoutCalorimeter::Init()
+{
+  ExRootConfParam param, paramEtaBins, paramPhiBins, paramFractions;
+  Long_t i, j, k, size, sizeEtaBins, sizePhiBins;
+  Double_t ecalFraction, hcalFraction;
+  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] = make_pair(0.0, 1.0);
+
+  for(i = 0; i < size/2; ++i)
+  {
+    paramFractions = param[i*2 + 1];
+
+    ecalFraction = paramFractions[0].GetDouble();
+    hcalFraction = paramFractions[1].GetDouble();
+
+    fFractionMap[param[i*2].GetInt()] = make_pair(ecalFraction, hcalFraction);
+  }
+
+  // read min E value for timing measurement in ECAL
+  fTimingEnergyMin = GetDouble("TimingEnergyMin",4.);
+  // For timing
+  // So far this flag needs to be false
+  // Curved extrapolation not supported
+  fElectronsFromTrack = false;
+
+  // read min E value for towers to be saved
+  fECalEnergyMin = GetDouble("ECalEnergyMin", 0.0);
+  fHCalEnergyMin = GetDouble("HCalEnergyMin", 0.0);
+
+  fECalEnergySignificanceMin = GetDouble("ECalEnergySignificanceMin", 0.0);
+  fHCalEnergySignificanceMin = GetDouble("HCalEnergySignificanceMin", 0.0);
+
+  // switch on or off the dithering of the center of DualReadoutCalorimeter towers
+  fSmearTowerCenter = GetBool("SmearTowerCenter", true);
+
+  // read resolution formulas
+  fECalResolutionFormula->Compile(GetString("ECalResolutionFormula", "0"));
+  fHCalResolutionFormula->Compile(GetString("HCalResolutionFormula", "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"));
+  fPhotonOutputArray = ExportArray(GetString("PhotonOutputArray", "photons"));
+
+  fEFlowTrackOutputArray = ExportArray(GetString("EFlowTrackOutputArray", "eflowTracks"));
+  fEFlowPhotonOutputArray = ExportArray(GetString("EFlowPhotonOutputArray", "eflowPhotons"));
+  fEFlowNeutralHadronOutputArray = ExportArray(GetString("EFlowNeutralHadronOutputArray", "eflowNeutralHadrons"));
+}
+
+//------------------------------------------------------------------------------
+
+void DualReadoutCalorimeter::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 DualReadoutCalorimeter::Process()
+{
+  Candidate *particle, *track;
+  TLorentzVector position, momentum;
+  Short_t etaBin, phiBin, flags;
+  Int_t number;
+  Long64_t towerHit, towerEtaPhi, hitEtaPhi;
+  Double_t ecalFraction, hcalFraction;
+  Double_t ecalEnergy, hcalEnergy;
+  Double_t ecalSigma, hcalSigma;
+  Double_t energyGuess;
+  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();
+  fECalTowerFractions.clear();
+  fHCalTowerFractions.clear();
+  fECalTrackFractions.clear();
+  fHCalTrackFractions.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);
+    }
+
+    ecalFraction = itFractionMap->second.first;
+    hcalFraction = itFractionMap->second.second;
+
+    fECalTowerFractions.push_back(ecalFraction);
+    fHCalTowerFractions.push_back(hcalFraction);
+
+    if(ecalFraction < 1.0E-9 && hcalFraction < 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);
+    }
+
+    ecalFraction = itFractionMap->second.first;
+    hcalFraction = itFractionMap->second.second;
+
+    fECalTrackFractions.push_back(ecalFraction);
+    fHCalTrackFractions.push_back(hcalFraction);
+
+    // 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];
+
+      fECalTowerEnergy = 0.0;
+      fHCalTowerEnergy = 0.0;
+
+      fECalTrackEnergy = 0.0;
+      fHCalTrackEnergy = 0.0;
+      
+      fECalTrackSigma = 0.0;
+      fHCalTrackSigma = 0.0;
+      
+      fTowerTrackHits = 0;
+      fTowerPhotonHits = 0;
+
+      fECalTowerTrackArray->Clear();
+      fHCalTowerTrackArray->Clear();
+    
+    }
+
+    // check for track hits
+    if(flags & 1)
+    {
+      ++fTowerTrackHits;
+
+      track = static_cast<Candidate*>(fTrackInputArray->At(number));
+      momentum = track->Momentum;
+      position = track->Position;
+
+      ecalEnergy = momentum.E() * fECalTrackFractions[number];
+      hcalEnergy = momentum.E() * fHCalTrackFractions[number];
+
+      if(ecalEnergy > fTimingEnergyMin && fTower)
+      {
+        if(fElectronsFromTrack)
+        {
+          fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, track->Position.T()));
+        }
+      }
+
+      if(fECalTrackFractions[number] > 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
+      {
+        fECalTrackEnergy += ecalEnergy;
+        ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());        
+        if(ecalSigma/momentum.E() < track->TrackResolution) energyGuess = ecalEnergy;        
+        else energyGuess = momentum.E();
+
+        fECalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fECalTowerTrackArray->Add(track);
+      }
+     
+      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] > 1.0E-9)
+      {
+        fHCalTrackEnergy += hcalEnergy;
+        hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, momentum.E());
+        if(hcalSigma/momentum.E() < track->TrackResolution) energyGuess = hcalEnergy;
+        else energyGuess = momentum.E();
+
+        fHCalTrackSigma += (track->TrackResolution)*energyGuess*(track->TrackResolution)*energyGuess;
+        fHCalTowerTrackArray->Add(track);
+      }
+      
+      else if(fECalTrackFractions[number] < 1.0E-9 && fHCalTrackFractions[number] < 1.0E-9)
+      {
+        fEFlowTrackOutputArray->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
+    ecalEnergy = momentum.E() * fECalTowerFractions[number];
+    hcalEnergy = momentum.E() * fHCalTowerFractions[number];
+
+    fECalTowerEnergy += ecalEnergy;
+    fHCalTowerEnergy += hcalEnergy;
+
+    if(ecalEnergy > fTimingEnergyMin && fTower)
+    {
+      if (abs(particle->PID) != 11 || !fElectronsFromTrack)
+      {
+        fTower->ECalEnergyTimePairs.push_back(make_pair<Float_t, Float_t>(ecalEnergy, particle->Position.T()));
+      }
+    }
+
+    fTower->AddCandidate(particle);
+  }
+
+  // finalize last tower
+  FinalizeTower();
+}
+
+//------------------------------------------------------------------------------
+
+void DualReadoutCalorimeter::FinalizeTower()
+{
+
+
+
+
+  Candidate *track, *tower, *mother;
+  Double_t energy, pt, eta, phi;
+  Double_t ecalEnergy, hcalEnergy;
+  Double_t ecalNeutralEnergy, hcalNeutralEnergy;
+  
+  Double_t ecalSigma, hcalSigma, sigma;
+  Double_t ecalNeutralSigma, hcalNeutralSigma;
+
+  Double_t weightTrack, weightCalo, bestEnergyEstimate, rescaleFactor;
+  
+  TLorentzVector momentum;
+  TFractionMap::iterator itFractionMap;
+
+  Float_t weight, sumWeightedTime, sumWeight;
+
+  if(!fTower) return;
+
+
+  //if (fHCalTowerEnergy < 30 && fECalTowerEnergy < 30) return;
+  //cout<<"----------- New tower ---------"<<endl;
+
+
+  // here we change behaviour w.r.t to standard calorimeter. Start with worse case scenario. If fHCalTowerEnergy > 0, assume total energy smeared by HCAL resolution.
+  // For example, if overlapping charged pions and photons take hadronic resolution as combined measurement
+
+  // if no hadronic fraction at all, then use ECAL resolution
+
+  //cout<<"fECalTowerEnergy: "<<fECalTowerEnergy<<", fHCalTowerEnergy: "<<fHCalTowerEnergy<<endl;
+
+  if (fHCalTowerEnergy <= fHCalEnergyMin)
+  {
+    energy = fECalTowerEnergy;
+    sigma  = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
+    //cout<<"em: "<<energy<<","<<sigma<<endl;
+  }
+
+  // if hadronic fraction > 0, use HCAL resolution 
+  else
+  {
+    energy = fECalTowerEnergy + fHCalTowerEnergy;
+    sigma  = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, energy);
+    //cout<<"had: "<<energy<<","<<sigma<<endl;
+  }
+
+  energy = LogNormal(energy, sigma);
+  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
+  
+  if(energy < fEnergyMin || energy < fEnergySignificanceMin*sigma) energy = 0.0;
+
+  // for now keep this the same
+  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fECalTowerEnergy);
+  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, fHCalTowerEnergy);
+
+  ecalEnergy = LogNormal(fECalTowerEnergy, ecalSigma);
+  hcalEnergy = LogNormal(fHCalTowerEnergy, hcalSigma);
+
+  ecalSigma = fECalResolutionFormula->Eval(0.0, fTowerEta, 0.0, ecalEnergy);
+  hcalSigma = fHCalResolutionFormula->Eval(0.0, fTowerEta, 0.0, hcalEnergy);
+
+  if(ecalEnergy < fECalEnergyMin || ecalEnergy < fECalEnergySignificanceMin*ecalSigma) ecalEnergy = 0.0;
+  if(hcalEnergy < fHCalEnergyMin || hcalEnergy < fHCalEnergySignificanceMin*hcalSigma) hcalEnergy = 0.0;
+  //cout<<energy<<","<<ecalEnergy<<","<<hcalEnergy<<endl;
+
+
+  if(fSmearTowerCenter)
+  {
+    eta = gRandom->Uniform(fTowerEdges[0], fTowerEdges[1]);
+    phi = gRandom->Uniform(fTowerEdges[2], fTowerEdges[3]);
+  }
+  else
+  {
+    eta = fTowerEta;
+    phi = fTowerPhi;
+  }
+
+  pt = energy / TMath::CosH(eta);
+
+  // Time calculation for tower
+  fTower->NTimeHits = 0;
+  sumWeightedTime = 0.0;
+  sumWeight = 0.0;
+
+  for(size_t i = 0; i < fTower->ECalEnergyTimePairs.size(); ++i)
+  {
+    weight = TMath::Sqrt(fTower->ECalEnergyTimePairs[i].first);
+    sumWeightedTime += weight * fTower->ECalEnergyTimePairs[i].second;
+    sumWeight += weight;
+    fTower->NTimeHits++;
+  }
+
+  if(sumWeight > 0.0)
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, sumWeightedTime/sumWeight);
+  }
+  else
+  {
+    fTower->Position.SetPtEtaPhiE(1.0, eta, phi, 999999.9);
+  }
+
+
+  fTower->Momentum.SetPtEtaPhiE(pt, eta, phi, energy);
+  fTower->Eem = ecalEnergy;
+  fTower->Ehad = hcalEnergy;
+
+  fTower->Edges[0] = fTowerEdges[0];
+  fTower->Edges[1] = fTowerEdges[1];
+  fTower->Edges[2] = fTowerEdges[2];
+  fTower->Edges[3] = fTowerEdges[3];
+
+  if(energy > 0.0)
+  {
+    if(fTowerPhotonHits > 0 && fTowerTrackHits == 0)
+    {
+      fPhotonOutputArray->Add(fTower);
+    }
+
+    fTowerOutputArray->Add(fTower);
+  }
+  
+  // fill energy flow candidates
+  fECalTrackSigma = TMath::Sqrt(fECalTrackSigma);
+  fHCalTrackSigma = TMath::Sqrt(fHCalTrackSigma);
+
+  //compute neutral excesses
+  ecalNeutralEnergy = max( (ecalEnergy - fECalTrackEnergy) , 0.0);
+  hcalNeutralEnergy = max( (hcalEnergy - fHCalTrackEnergy) , 0.0);
+  
+  ecalNeutralSigma = ecalNeutralEnergy / TMath::Sqrt(fECalTrackSigma*fECalTrackSigma + ecalSigma*ecalSigma);
+  hcalNeutralSigma = hcalNeutralEnergy / TMath::Sqrt(fHCalTrackSigma*fHCalTrackSigma + hcalSigma*hcalSigma);
+  
+   // if ecal neutral excess is significant, simply create neutral EflowPhoton tower and clone each track into eflowtrack
+  if(ecalNeutralEnergy > fECalEnergyMin && ecalNeutralSigma > fECalEnergySignificanceMin)
+  {
+    // create new photon tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  ecalNeutralEnergy / TMath::CosH(eta);
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, ecalNeutralEnergy);
+    tower->Eem = ecalNeutralEnergy;
+    tower->Ehad = 0.0;
+    tower->PID = 22;
+    
+    fEFlowPhotonOutputArray->Add(tower);
+   
+    //clone tracks
+    fItECalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  
+  }
+ 
+  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
+  else if(fECalTrackEnergy > 0.0)
+  {
+    weightTrack = (fECalTrackSigma > 0.0) ? 1 / (fECalTrackSigma*fECalTrackSigma) : 0.0;
+    weightCalo  = (ecalSigma > 0.0) ? 1 / (ecalSigma*ecalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fECalTrackEnergy + weightCalo*ecalEnergy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate/fECalTrackEnergy;
+
+    //rescale tracks
+    fItECalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItECalTowerTrackArray->Next())))
+    {  
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+
+
+  // if hcal neutral excess is significant, simply create neutral EflowNeutralHadron tower and clone each track into eflowtrack
+  if(hcalNeutralEnergy > fHCalEnergyMin && hcalNeutralSigma > fHCalEnergySignificanceMin)
+  {
+    // create new photon tower
+    tower = static_cast<Candidate*>(fTower->Clone());
+    pt =  hcalNeutralEnergy / TMath::CosH(eta);
+    
+    tower->Momentum.SetPtEtaPhiE(pt, eta, phi, hcalNeutralEnergy);
+    tower->Ehad = hcalNeutralEnergy;
+    tower->Eem = 0.0;
+    
+    fEFlowNeutralHadronOutputArray->Add(tower);
+   
+    //clone tracks
+    fItHCalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
+    {
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  
+  }
+ 
+  // if neutral excess is not significant, rescale eflow tracks, such that the total charged equals the best measurement given by the DualReadoutCalorimeter and tracking
+  else if(fHCalTrackEnergy > 0.0)
+  {
+    weightTrack = (fHCalTrackSigma > 0.0) ? 1 / (fHCalTrackSigma*fHCalTrackSigma) : 0.0;
+    weightCalo  = (hcalSigma > 0.0) ? 1 / (hcalSigma*hcalSigma) : 0.0;
+  
+    bestEnergyEstimate = (weightTrack*fHCalTrackEnergy + weightCalo*hcalEnergy) / (weightTrack + weightCalo); 
+    rescaleFactor = bestEnergyEstimate / fHCalTrackEnergy;
+
+    //rescale tracks
+    fItHCalTowerTrackArray->Reset();
+    while((track = static_cast<Candidate*>(fItHCalTowerTrackArray->Next())))
+    {  
+      mother = track;
+      track = static_cast<Candidate*>(track->Clone());
+      track->AddCandidate(mother);
+
+      track->Momentum *= rescaleFactor;
+
+      fEFlowTrackOutputArray->Add(track);
+    }
+  }
+  
+  
+}
+
+//------------------------------------------------------------------------------
+
+Double_t DualReadoutCalorimeter::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.0, 1.0));
+  }
+  else
+  {
+    return 0.0;
+  }
+}
Index: modules/DualReadoutCalorimeter.h
===================================================================
--- modules/DualReadoutCalorimeter.h	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
+++ modules/DualReadoutCalorimeter.h	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
@@ -0,0 +1,122 @@
+/*
+ *  Delphes: a framework for fast simulation of a generic collider experiment
+ *  Copyright (C) 2012-2014  Universite catholique de Louvain (UCL), Belgium
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU General Public License as published by
+ *  the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program is distributed in the hope that it will be useful,
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *  GNU General Public License for more details.
+ *
+ *  You should have received a copy of the GNU General Public License
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef DualReadoutCalorimeter_h
+#define DualReadoutCalorimeter_h
+
+/** \class DualReadoutCalorimeter
+ *
+ *  Fills DualReadoutCalorimeter towers, performs DualReadoutCalorimeter resolution smearing,
+ *  and creates energy flow objects (tracks, photons, and neutral hadrons).
+ *
+ *  \author P. Demin - UCL, Louvain-la-Neuve
+ *
+ */
+
+#include "classes/DelphesModule.h"
+
+#include <map>
+#include <set>
+#include <vector>
+
+class TObjArray;
+class DelphesFormula;
+class Candidate;
+
+class DualReadoutCalorimeter: public DelphesModule
+{
+public:
+
+  DualReadoutCalorimeter();
+  ~DualReadoutCalorimeter();
+
+  void Init();
+  void Process();
+  void Finish();
+
+private:
+
+  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 fECalTowerEnergy, fHCalTowerEnergy;
+  Double_t fECalTrackEnergy, fHCalTrackEnergy;
+
+  Double_t fTimingEnergyMin;
+  Bool_t fElectronsFromTrack;
+
+  Int_t fTowerTrackHits, fTowerPhotonHits;
+
+  Double_t fECalEnergyMin;
+  Double_t fHCalEnergyMin;
+  Double_t fEnergyMin;
+
+  Double_t fECalEnergySignificanceMin;
+  Double_t fHCalEnergySignificanceMin;
+  Double_t fEnergySignificanceMin;
+
+  Double_t fECalTrackSigma;
+  Double_t fHCalTrackSigma;
+
+  Bool_t fSmearTowerCenter;
+
+  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 > fECalTowerFractions;
+  std::vector < Double_t > fHCalTowerFractions;
+
+  std::vector < Double_t > fECalTrackFractions;
+  std::vector < Double_t > fHCalTrackFractions;
+
+  DelphesFormula *fECalResolutionFormula; //!
+  DelphesFormula *fHCalResolutionFormula; //!
+
+  TIterator *fItParticleInputArray; //!
+  TIterator *fItTrackInputArray; //!
+
+  const TObjArray *fParticleInputArray; //!
+  const TObjArray *fTrackInputArray; //!
+
+  TObjArray *fTowerOutputArray; //!
+  TObjArray *fPhotonOutputArray; //!
+
+  TObjArray *fEFlowTrackOutputArray; //!
+  TObjArray *fEFlowPhotonOutputArray; //!
+  TObjArray *fEFlowNeutralHadronOutputArray; //!
+
+  TObjArray *fECalTowerTrackArray; //!
+  TIterator *fItECalTowerTrackArray; //!
+
+  TObjArray *fHCalTowerTrackArray; //!
+  TIterator *fItHCalTowerTrackArray; //!
+
+  void FinalizeTower();
+  Double_t LogNormal(Double_t mean, Double_t sigma);
+
+  ClassDef(DualReadoutCalorimeter, 1)
+};
+
+#endif
Index: modules/ModulesLinkDef.h
===================================================================
--- modules/ModulesLinkDef.h	(revision d20cefb7bcdeff55de7ffee4708dbdf3a74c9d96)
+++ modules/ModulesLinkDef.h	(revision ef8a06df96aad24050b236455e8bf873bf2f0b32)
@@ -41,4 +41,5 @@
 #include "modules/DenseTrackFilter.h"
 #include "modules/Calorimeter.h"
+#include "modules/DualReadoutCalorimeter.h"
 #include "modules/OldCalorimeter.h"
 #include "modules/Isolation.h"
@@ -94,4 +95,5 @@
 #pragma link C++ class DenseTrackFilter+;
 #pragma link C++ class Calorimeter+;
+#pragma link C++ class DualReadoutCalorimeter+;
 #pragma link C++ class OldCalorimeter+;
 #pragma link C++ class Isolation+;
