Skip to content

Commit

Permalink
Merge pull request #17691 from PFCal-dev/hgc-tpg-integration-90X-170228
Browse files Browse the repository at this point in the history
HGCAL trigger sim clustering
  • Loading branch information
davidlange6 committed Mar 9, 2017
2 parents 588950a + 7dc0f7c commit d6d088f
Show file tree
Hide file tree
Showing 26 changed files with 966 additions and 87 deletions.
77 changes: 77 additions & 0 deletions DataFormats/L1THGCal/interface/ClusterShapes.h
@@ -0,0 +1,77 @@
#ifndef DataFormats_L1HGCal_ClusterShapes_H
#define DataFormats_L1HGCal_ClusterShapes_H
#include <cmath>

namespace l1t{
// this class is design to contain and compute
// efficiently the cluster shapes
// running only once on the cluster members.
class ClusterShapes{
private:
float sum_e_ = 0.0;
float sum_e2_ = 0.0;
float sum_logE_ = 0.0;
int n_=0.0;

float emax_ = 0.0;

float sum_w_ =0.0; // just here for clarity
float sum_eta_ = 0.0;
float sum_r_ = 0.0;
// i will discriminate using the rms in -pi,pi or in 0,pi
float sum_phi_0_ = 0.0; // computed in -pi,pi
float sum_phi_1_ = 0.0; // computed in 0, 2pi

float sum_eta2_=0.0;
float sum_r2_ = 0.0;
float sum_phi2_0_=0.0; //computed in -pi,pi
float sum_phi2_1_=0.0; //computed in 0,2pi

// off diagonal element of the tensor
float sum_eta_r_ =0.0;
float sum_r_phi_0_ = 0.0;
float sum_r_phi_1_ = 0.0;
float sum_eta_phi_0_ = 0.0;
float sum_eta_phi_1_ = 0.0;

// caching of informations
mutable bool isPhi0_ = true;
mutable bool modified_ = false; // check wheneever i need

public:
ClusterShapes(){}
ClusterShapes(float e, float eta, float phi, float r) { Init(e,eta,phi,r);}
~ClusterShapes(){}
ClusterShapes(const ClusterShapes&x)= default;
//init an empty cluster
void Init(float e, float eta, float phi, float r=0.);
inline void Add(float e, float eta, float phi, float r=0.0)
{ (*this) += ClusterShapes(e,eta,phi,r);}


// ---
float SigmaEtaEta() const ;
float SigmaPhiPhi() const ;
float SigmaRR() const ;
// ----
float Phi() const ;
float R() const ;
float Eta() const ;
inline int N()const {return n_;}
// --
float SigmaEtaR()const ;
float SigmaEtaPhi() const ;
float SigmaRPhi() const ;
// --
float LogEoverE() const { return sum_logE_/sum_e_;}
float eD() const { return std::sqrt(sum_e2_)/sum_e_;}

ClusterShapes operator+(const ClusterShapes &);
void operator+=(const ClusterShapes &);
ClusterShapes& operator=(const ClusterShapes &) = default;
};

}; // end namespace

#endif

5 changes: 4 additions & 1 deletion DataFormats/L1THGCal/interface/HGCalCluster.h
Expand Up @@ -3,9 +3,10 @@

#include "DataFormats/L1Trigger/interface/L1Candidate.h"
#include "DataFormats/L1Trigger/interface/BXVector.h"
#include "DataFormats/L1THGCal/interface/ClusterShapes.h"

namespace l1t {

class HGCalCluster : public L1Candidate {
public:
HGCalCluster(){}
Expand Down Expand Up @@ -36,6 +37,8 @@ namespace l1t {

uint32_t hOverE() const {return hOverE_;}

ClusterShapes shapes ;

bool operator<(const HGCalCluster& cl) const;
bool operator>(const HGCalCluster& cl) const {return cl<*this;};
bool operator<=(const HGCalCluster& cl) const {return !(cl>*this);};
Expand Down
132 changes: 132 additions & 0 deletions DataFormats/L1THGCal/src/ClusterShapes.cc
@@ -0,0 +1,132 @@
#include "DataFormats/L1THGCal/interface/ClusterShapes.h"
#include <cmath>

using namespace l1t;

ClusterShapes ClusterShapes::operator+(const ClusterShapes& x)
{
ClusterShapes cs(*this); // copy constructor
cs += x;
return cs;
}


void ClusterShapes::operator +=(const ClusterShapes &x){

sum_e_ += x.sum_e_;
sum_e2_ += x.sum_e2_;
sum_logE_ += x.sum_logE_;
n_ += x.n_;

sum_w_ += x.sum_w_;

emax_ = (emax_> x.emax_) ? emax_: x.emax_;

// mid-point
sum_eta_ += x.sum_eta_;
sum_phi_0_ += x.sum_phi_0_; //
sum_phi_1_ += x.sum_phi_1_; //
sum_r_ += x.sum_r_;

// square
sum_eta2_ += x.sum_eta2_;
sum_phi2_0_ += x.sum_phi2_0_;
sum_phi2_1_ += x.sum_phi2_1_;
sum_r2_ += x.sum_r2_;

// off diagonal
sum_eta_r_ += x.sum_eta_r_ ;
sum_r_phi_0_ += x.sum_r_phi_0_ ;
sum_r_phi_1_ += x.sum_r_phi_1_ ;
sum_eta_phi_0_ += x.sum_eta_phi_0_;
sum_eta_phi_1_ += x.sum_eta_phi_1_;

}


// -------------- CLUSTER SHAPES ---------------
void ClusterShapes::Init(float e ,float eta, float phi, float r){
if (e<=0 ) return;
sum_e_ = e;
sum_e2_ = e*e;
sum_logE_ = std::log(e);

float w = e;

n_=1;

sum_w_ = w;

sum_phi_0_ = w *( phi );
sum_phi_1_ = w* (phi + M_PI);
sum_r_ = w * r;
sum_eta_ = w * eta;

//--
sum_r2_ += w * (r*r);
sum_eta2_ += w * (eta*eta);
sum_phi2_0_ += w * (phi*phi);
sum_phi2_1_ += w * (phi+M_PI)*(phi+M_PI);

// -- off diagonal
sum_eta_r_ += w * (r*eta);
sum_r_phi_0_ += w* (r *phi);
sum_r_phi_1_ += w* r *(phi + M_PI);
sum_eta_phi_0_ += w* (eta *phi);
sum_eta_phi_1_ += w* eta * (phi+M_PI);

}
// ------
float ClusterShapes::Eta()const { return sum_eta_/sum_w_;}
float ClusterShapes::R() const { return sum_r_/sum_w_;}

float ClusterShapes::SigmaEtaEta()const {return sum_eta2_/sum_w_ - Eta()*Eta();}

float ClusterShapes::SigmaRR()const { return sum_r2_/sum_w_ - R() *R();}


float ClusterShapes::SigmaPhiPhi()const {
float phi_0 = (sum_phi_0_ / sum_w_);
float phi_1 = (sum_phi_1_ / sum_w_);
float spp_0 = sum_phi2_0_ / sum_w_ - phi_0*phi_0;
float spp_1 = sum_phi2_1_ / sum_w_ - phi_1*phi_1;

if (spp_0 < spp_1 )
{
isPhi0_=true;
return spp_0;
}
else
{
isPhi0_=false;
return spp_1;
}
}

float ClusterShapes::Phi()const {
SigmaPhiPhi(); //update phi
if (isPhi0_) return (sum_phi_0_ / sum_w_);
else return (sum_phi_1_ / sum_w_);
}


// off - diagonal
float ClusterShapes::SigmaEtaR() const { return -(sum_eta_r_ / sum_w_ - Eta() *R()) ;}

float ClusterShapes::SigmaEtaPhi()const {
SigmaPhiPhi() ; // decide which phi use, update phi

if (isPhi0_)
return -(sum_eta_phi_0_ /sum_w_ - Eta()*(sum_phi_0_ / sum_w_));
else
return -(sum_eta_phi_1_ / sum_w_ - Eta()*(sum_phi_1_ / sum_w_));
}

float ClusterShapes::SigmaRPhi()const {
SigmaPhiPhi() ; // decide which phi use, update phi
if (isPhi0_)
return -(sum_r_phi_0_ / sum_w_ - R() *(sum_phi_0_ / sum_w_));
else
return -(sum_r_phi_1_ / sum_w_ - R() * (sum_phi_1_ / sum_w_));
}

4 changes: 4 additions & 0 deletions DataFormats/L1THGCal/src/classes.h
Expand Up @@ -5,6 +5,8 @@
#include "DataFormats/L1THGCal/interface/HGCalTower.h"
#include "DataFormats/L1THGCal/interface/HGCalTriggerCell.h"

#include "DataFormats/L1THGCal/interface/ClusterShapes.h"

namespace DataFormats {
namespace L1THGCal {
l1t::HGCFETriggerDigi hgcfetd;
Expand All @@ -27,5 +29,7 @@ namespace DataFormats {
edm::Wrapper<l1t::HGCalTowerBxCollection> w_hgcalTowerBxColl;
edm::Wrapper<l1t::HGCalClusterBxCollection> w_hgcalClusterBxColl;
edm::Wrapper<l1t::HGCalTriggerCellBxCollection> w_hgcalTriggerCellBxColl;

l1t::ClusterShapes clusterShapes;
}
}
2 changes: 2 additions & 0 deletions DataFormats/L1THGCal/src/classes_def.xml
Expand Up @@ -29,4 +29,6 @@
<class name="l1t::HGCalTriggerCellBxCollection"/>
<class name="edm::Wrapper<l1t::HGCalTriggerCellBxCollection>"/>

<class name="l1t::ClusterShapes"/>

</lcgdict>
17 changes: 11 additions & 6 deletions L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h
Expand Up @@ -29,10 +29,12 @@

class HGCalTriggerBackendAlgorithmBase {
public:
HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf) :
// Allow HGCalTriggerBackend to be passed a consume collector
HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf, edm::ConsumesCollector &cc) :
geometry_(nullptr),
name_(conf.getParameter<std::string>("AlgorithmName"))
{}

virtual ~HGCalTriggerBackendAlgorithmBase() {}

const std::string& name() const { return name_; }
Expand All @@ -42,7 +44,10 @@ class HGCalTriggerBackendAlgorithmBase {
//runs the trigger algorithm, storing internally the results
virtual void setProduces(edm::EDProducer& prod) const = 0;

virtual void run(const l1t::HGCFETriggerDigiCollection& coll, const edm::EventSetup& es) = 0;
virtual void run(const l1t::HGCFETriggerDigiCollection& coll,
const edm::EventSetup& es,
const edm::Event &e
) = 0;

virtual void putInEvent(edm::Event& evt) = 0;

Expand All @@ -61,9 +66,9 @@ namespace HGCalTriggerBackend {
template<typename FECODEC>
class Algorithm : public HGCalTriggerBackendAlgorithmBase {
public:
Algorithm(const edm::ParameterSet& conf) :
HGCalTriggerBackendAlgorithmBase(conf),
codec_(conf.getParameterSet("FECodec")){ }
Algorithm(const edm::ParameterSet& conf, edm::ConsumesCollector &cc ) :
HGCalTriggerBackendAlgorithmBase(conf, cc),
codec_(conf.getParameterSet("FECodec")){ }

virtual void setGeometry(const HGCalTriggerGeometryBase* const geom) override final {
HGCalTriggerBackendAlgorithmBase::setGeometry(geom);
Expand All @@ -76,6 +81,6 @@ namespace HGCalTriggerBackend {
}

#include "FWCore/PluginManager/interface/PluginFactory.h"
typedef edmplugin::PluginFactory< HGCalTriggerBackendAlgorithmBase* (const edm::ParameterSet&) > HGCalTriggerBackendAlgorithmFactory;
typedef edmplugin::PluginFactory< HGCalTriggerBackendAlgorithmBase* (const edm::ParameterSet&,edm::ConsumesCollector & ) > HGCalTriggerBackendAlgorithmFactory;

#endif
7 changes: 5 additions & 2 deletions L1Trigger/L1THGCal/interface/HGCalTriggerBackendProcessor.h
Expand Up @@ -31,13 +31,16 @@ class HGCalTriggerBackendProcessor {
public:
typedef std::unique_ptr<HGCalTriggerBackendAlgorithmBase> algo_ptr;

HGCalTriggerBackendProcessor(const edm::ParameterSet& conf);
HGCalTriggerBackendProcessor(const edm::ParameterSet& conf, edm::ConsumesCollector&&cc);

void setGeometry(const HGCalTriggerGeometryBase* const geom);

void setProduces(edm::EDProducer& prod) const;

void run(const l1t::HGCFETriggerDigiCollection& coll, const edm::EventSetup& es);
void run(const l1t::HGCFETriggerDigiCollection& coll,
const edm::EventSetup& es,
const edm::Event&e
);

void putInEvent(edm::Event& evt);

Expand Down
1 change: 1 addition & 0 deletions L1Trigger/L1THGCal/plugins/BuildFile.xml
Expand Up @@ -21,6 +21,7 @@
<use name="rootcore"/>
<use name="L1Trigger/L1THGCal"/>
<use name="DataFormats/L1THGCal"/>
<use name="SimDataFormats/CaloTest"/>
<flags EDM_PLUGIN="1"/>
</library>

Expand Down
7 changes: 3 additions & 4 deletions L1Trigger/L1THGCal/plugins/HGCalTriggerDigiFEReproducer.cc
Expand Up @@ -38,7 +38,8 @@ DEFINE_FWK_MODULE(HGCalTriggerDigiFEReproducer);

/*****************************************************************/
HGCalTriggerDigiFEReproducer::HGCalTriggerDigiFEReproducer(const edm::ParameterSet& conf):
inputdigi_(consumes<l1t::HGCFETriggerDigiCollection>(conf.getParameter<edm::InputTag>("feDigis")))
inputdigi_(consumes<l1t::HGCFETriggerDigiCollection>(conf.getParameter<edm::InputTag>("feDigis"))),
backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"), consumesCollector()))
/*****************************************************************/
{
//setup FE codec
Expand All @@ -49,8 +50,6 @@ HGCalTriggerDigiFEReproducer::HGCalTriggerDigiFEReproducer(const edm::ParameterS
codec_->unSetDataPayload();

produces<l1t::HGCFETriggerDigiCollection>();
//setup BE processor
backEndProcessor_ = std::make_unique<HGCalTriggerBackendProcessor>(conf.getParameterSet("BEConfiguration"));
// register backend processor products
backEndProcessor_->setProduces(*this);
}
Expand Down Expand Up @@ -98,7 +97,7 @@ void HGCalTriggerDigiFEReproducer::produce(edm::Event& e, const edm::EventSetup&
auto fe_digis_coll = *fe_digis_handle;

//now we run the emulation of the back-end processor
backEndProcessor_->run(fe_digis_coll, es);
backEndProcessor_->run(fe_digis_coll,es,e);
backEndProcessor_->putInEvent(e);
backEndProcessor_->reset();
}
13 changes: 5 additions & 8 deletions L1Trigger/L1THGCal/plugins/HGCalTriggerDigiProducer.cc
Expand Up @@ -38,11 +38,10 @@ DEFINE_FWK_MODULE(HGCalTriggerDigiProducer);
HGCalTriggerDigiProducer::
HGCalTriggerDigiProducer(const edm::ParameterSet& conf):
inputee_(consumes<HGCEEDigiCollection>(conf.getParameter<edm::InputTag>("eeDigis"))),
inputfh_(consumes<HGCHEDigiCollection>(conf.getParameter<edm::InputTag>("fhDigis")))
//inputbh_(consumes<HGCHEDigiCollection>(conf.getParameter<edm::InputTag>("bhDigis")))
inputfh_(consumes<HGCHEDigiCollection>(conf.getParameter<edm::InputTag>("fhDigis"))),
//inputbh_(consumes<HGCHEDigiCollection>(conf.getParameter<edm::InputTag>("bhDigis"))),
backEndProcessor_(new HGCalTriggerBackendProcessor(conf.getParameterSet("BEConfiguration"),consumesCollector()) )
{


//setup FE codec
const edm::ParameterSet& feCodecConfig = conf.getParameterSet("FECodec");
const std::string& feCodecName = feCodecConfig.getParameter<std::string>("CodecName");
Expand All @@ -51,8 +50,6 @@ HGCalTriggerDigiProducer(const edm::ParameterSet& conf):
codec_->unSetDataPayload();

produces<l1t::HGCFETriggerDigiCollection>();
//setup BE processor
backEndProcessor_ = std::make_unique<HGCalTriggerBackendProcessor>(conf.getParameterSet("BEConfiguration"));
// register backend processor products
backEndProcessor_->setProduces(*this);
}
Expand Down Expand Up @@ -130,7 +127,7 @@ void HGCalTriggerDigiProducer::produce(edm::Event& e, const edm::EventSetup& es)
auto fe_digis_coll = *fe_digis_handle;

//now we run the emulation of the back-end processor
backEndProcessor_->run(fe_digis_coll, es);
backEndProcessor_->reset();
backEndProcessor_->run(fe_digis_coll,es,e);
backEndProcessor_->putInEvent(e);
backEndProcessor_->reset();
}

0 comments on commit d6d088f

Please sign in to comment.