ClusterShapes.cc

Go to the documentation of this file.

#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_));
}