Skip to content

Commit

Permalink
Merge pull request #79 from amarini/topic_simclustering_cmssw90
Browse files Browse the repository at this point in the history
Sim Clustering
  • Loading branch information
jbsauvan authored Feb 13, 2017
2 parents d51ae80 + ea0a3cd commit 3f43338
Show file tree
Hide file tree
Showing 19 changed files with 901 additions and 60 deletions.
85 changes: 85 additions & 0 deletions DataFormats/L1THGCal/interface/ClusterShapes.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,85 @@
#ifndef CLUSTER_SHAPES_H
#define CLUSTER_SHAPES_H
#include <vector>
#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);}


// --- this is what I want out:
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

// -----------------------------------
// Local Variables:
// mode:c++
// indent-tabs-mode:nil
// tab-width:4
// c-basic-offset:4
// End:
// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
6 changes: 5 additions & 1 deletion DataFormats/L1THGCal/interface/HGCalCluster.h
Original file line number Diff line number Diff line change
Expand Up @@ -3,9 +3,11 @@

#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 +38,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
146 changes: 146 additions & 0 deletions DataFormats/L1THGCal/src/ClusterShapes.cc
Original file line number Diff line number Diff line change
@@ -0,0 +1,146 @@
#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 )
{
float phi = phi_0;
isPhi0=true;
while (phi < - M_PI) phi += 2*M_PI;
while (phi > M_PI) phi -= 2*M_PI;
return spp_0;
}
else
{
float phi = phi_1 ;
isPhi0=false;
while (phi < - M_PI) phi += 2*M_PI;
while (phi > M_PI) phi -= 2*M_PI;
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));
}

// -----------------------------------
// Local Variables:
// mode:c++
// indent-tabs-mode:nil
// tab-width:4
// c-basic-offset:4
// End:
// vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
1 change: 1 addition & 0 deletions DataFormats/L1THGCal/src/HGCalCluster.cc
Original file line number Diff line number Diff line change
Expand Up @@ -31,3 +31,4 @@ bool HGCalCluster::operator<(const HGCalCluster& cl) const
}
return res;
}

4 changes: 4 additions & 0 deletions DataFormats/L1THGCal/src/classes.h
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
Expand Up @@ -29,4 +29,6 @@
<class name="l1t::HGCalTriggerCellBxCollection"/>
<class name="edm::Wrapper<l1t::HGCalTriggerCellBxCollection>"/>

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

</lcgdict>
26 changes: 22 additions & 4 deletions L1Trigger/L1THGCal/interface/HGCalTriggerBackendAlgorithmBase.h
Original file line number Diff line number Diff line change
Expand Up @@ -29,10 +29,19 @@

class HGCalTriggerBackendAlgorithmBase {
public:
/*
HGCalTriggerBackendAlgorithmBase(const edm::ParameterSet& conf) :
geometry_(nullptr),
name_(conf.getParameter<std::string>("AlgorithmName"))
{}
*/

// 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 +51,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 +73,15 @@ namespace HGCalTriggerBackend {
template<typename FECODEC>
class Algorithm : public HGCalTriggerBackendAlgorithmBase {
public:
/*
Algorithm(const edm::ParameterSet& conf) :
HGCalTriggerBackendAlgorithmBase(conf),
codec_(conf.getParameterSet("FECodec")){ }
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 +94,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
Original file line number Diff line number Diff line change
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
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ struct HGCalBestChoiceDataPayload
{
payload.fill(0);
}

};


Expand Down
1 change: 1 addition & 0 deletions L1Trigger/L1THGCal/plugins/BuildFile.xml
Original file line number Diff line number Diff line change
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
Loading

0 comments on commit 3f43338

Please sign in to comment.