PuppiContainer Class Reference

#include <PuppiContainer.h>

Classes
class	PuppiUserInfo

Public Member Functions
void	initialize (const std::vector< RecoObj > &iRecoObjects)
std::vector < fastjet::PseudoJet > const &	pfParticles () const
const std::vector< double > &	puppiAlphas ()
const std::vector< double > &	puppiAlphasMed ()
const std::vector< double > &	puppiAlphasRMS ()
	PuppiContainer (const edm::ParameterSet &iConfig)
int	puppiNAlgos ()
std::vector < fastjet::PseudoJet > const &	puppiParticles () const
const std::vector< double > &	puppiRawAlphas ()
std::vector< double > const &	puppiWeights ()
std::vector < fastjet::PseudoJet > const &	pvParticles () const
void	setNPV (int iNPV)
	~PuppiContainer ()

Protected Member Functions
double	getChi2FromdZ (double iDZ)
int	getPuppiId (float iPt, float iEta)
void	getRawAlphas (int iOpt, std::vector< fastjet::PseudoJet > const &iConstits, std::vector< fastjet::PseudoJet > const &iParticles, std::vector< fastjet::PseudoJet > const &iChargeParticles)
void	getRMSAvg (int iOpt, std::vector< fastjet::PseudoJet > const &iConstits, std::vector< fastjet::PseudoJet > const &iParticles, std::vector< fastjet::PseudoJet > const &iChargeParticles)
double	goodVar (fastjet::PseudoJet const &iPart, std::vector< fastjet::PseudoJet > const &iParts, int iOpt, double iRCone)
double	var_within_R (int iId, const std::vector< fastjet::PseudoJet > &particles, const fastjet::PseudoJet &centre, double R)

Protected Attributes
std::vector< double >	fAlphaMed
std::vector< double >	fAlphaRMS
bool	fApplyCHS
std::vector< fastjet::PseudoJet >	fChargedPV
bool	fInvert
int	fNAlgos
double	fNeutralMinPt
double	fNeutralSlope
int	fNPV
std::vector< fastjet::PseudoJet >	fPFParticles
std::vector< fastjet::PseudoJet >	fPupParticles
std::vector< PuppiAlgo >	fPuppiAlgo
bool	fPuppiDiagnostics
double	fPuppiWeightCut
double	fPVFrac
std::vector< double >	fRawAlphas
std::vector< RecoObj >	fRecoParticles
bool	fUseExp
std::vector< double >	fVals
std::vector< double >	fWeights

Detailed Description

Definition at line 13 of file PuppiContainer.h.

Constructor & Destructor Documentation

PuppiContainer::PuppiContainer

(

const edm::ParameterSet &

iConfig

)

Definition at line 14 of file PuppiContainer.cc.

References edm::ParameterSet::getParameter().

                                                             {
    fPuppiDiagnostics = iConfig.getParameter<bool>("puppiDiagnostics");
    fApplyCHS        = iConfig.getParameter<bool>("applyCHS");
    fInvert          = iConfig.getParameter<bool>("invertPuppi");    
    fUseExp          = iConfig.getParameter<bool>("useExp");
    fPuppiWeightCut  = iConfig.getParameter<double>("MinPuppiWeight");
    std::vector<edm::ParameterSet> lAlgos = iConfig.getParameter<std::vector<edm::ParameterSet> >("algos");
    fNAlgos = lAlgos.size();
    for(unsigned int i0 = 0; i0 < lAlgos.size(); i0++) {
        PuppiAlgo pPuppiConfig(lAlgos[i0]);
        fPuppiAlgo.push_back(pPuppiConfig);
    }
}

PuppiContainer::~PuppiContainer

(

)

Definition at line 74 of file PuppiContainer.cc.

{}

Member Function Documentation

double PuppiContainer::getChi2FromdZ ( double  iDZ )

protected

Definition at line 185 of file PuppiContainer.cc.

References funct::abs(), and alignCSCRings::e.

                                               {
    //We need to obtain prob of PU + (1-Prob of LV)
    // Prob(LV) = Gaus(dZ,sigma) where sigma = 1.5mm  (its really more like 1mm)
    //double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ),0.2)*2.; //*2 is to do it double sided
    //Take iDZ to be corrected by sigma already
    double lProbLV = ROOT::Math::normal_cdf_c(std::abs(iDZ),1.)*2.; //*2 is to do it double sided
    double lProbPU = 1-lProbLV;
    if(lProbPU <= 0) lProbPU = 1e-16;   //Quick Trick to through out infs
    if(lProbPU >= 0) lProbPU = 1-1e-16; //Ditto
    double lChi2PU = TMath::ChisquareQuantile(lProbPU,1);
    lChi2PU*=lChi2PU;
    return lChi2PU;
}

int PuppiContainer::getPuppiId ( float  iPt, float  iEta  )

protected

Definition at line 168 of file PuppiContainer.cc.

References funct::abs(), HLT_25ns14e33_v1_cff::etaMax, HLT_25ns14e33_v1_cff::etaMin, and PtMinSelector_cfg::ptMin.

                                                        {
    int lId = -1;
    for(int i0 = 0; i0 < fNAlgos; i0++) {
        int nEtaBinsPerAlgo = fPuppiAlgo[i0].etaBins();
        for (int i1 = 0; i1 < nEtaBinsPerAlgo; i1++){
            if ( (std::abs(iEta) > fPuppiAlgo[i0].etaMin(i1)) && (std::abs(iEta) < fPuppiAlgo[i0].etaMax(i1)) ){ 
                fPuppiAlgo[i0].fixAlgoEtaBin( i1 );
                if(iPt > fPuppiAlgo[i0].ptMin()){
                    lId = i0; 
                    break;
                }
            }
        }
    }
    //if(lId == -1) std::cerr << "Error : Full fiducial range is not defined " << std::endl;
    return lId;
}

void PuppiContainer::getRawAlphas ( int  iOpt, std::vector< fastjet::PseudoJet > const &  iConstits, std::vector< fastjet::PseudoJet > const &  iParticles, std::vector< fastjet::PseudoJet > const &  iChargeParticles  )

protected

Definition at line 149 of file PuppiContainer.cc.

References edm::isFinite(), and LogDebug.

                                                                                                                                                                                          {
    for(int j0 = 0; j0 < fNAlgos; j0++){
        for(unsigned int i0 = 0; i0 < iConstits.size(); i0++ ) {
            double pVal = -1;
            //Get the Puppi Sub Algo (given iteration)
            int  pAlgo    = fPuppiAlgo[j0].algoId   (iOpt);
            bool pCharged = fPuppiAlgo[j0].isCharged(iOpt);
            double pCone  = fPuppiAlgo[j0].coneSize (iOpt);
            //Compute the Puppi Metric
            if(!pCharged) pVal = goodVar(iConstits[i0],iParticles       ,pAlgo,pCone);
            if( pCharged) pVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
            fRawAlphas.push_back(pVal);
            if( ! edm::isFinite(pVal)) {
                LogDebug( "NotFound" )  << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt() << " -- " << iConstits[i0].eta() << endl;
                continue;
            }
        }
    }
}

void PuppiContainer::getRMSAvg ( int  iOpt, std::vector< fastjet::PseudoJet > const &  iConstits, std::vector< fastjet::PseudoJet > const &  iParticles, std::vector< fastjet::PseudoJet > const &  iChargeParticles  )

protected

Definition at line 109 of file PuppiContainer.cc.

References eta, edm::isFinite(), LogDebug, and EnergyCorrector::pt.

                                                                                                                                                                                       {
    for(unsigned int i0 = 0; i0 < iConstits.size(); i0++ ) {
        double pVal = -1;
        //Calculate the Puppi Algo to use
        int  pPupId   = getPuppiId(iConstits[i0].pt(),iConstits[i0].eta());
        if(pPupId == -1 || fPuppiAlgo[pPupId].numAlgos() <= iOpt){
            fVals.push_back(-1);
            continue;
        }
        //Get the Puppi Sub Algo (given iteration)
        int  pAlgo    = fPuppiAlgo[pPupId].algoId   (iOpt);
        bool pCharged = fPuppiAlgo[pPupId].isCharged(iOpt);
        double pCone  = fPuppiAlgo[pPupId].coneSize (iOpt);
        //Compute the Puppi Metric
        if(!pCharged) pVal = goodVar(iConstits[i0],iParticles       ,pAlgo,pCone);
        if( pCharged) pVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
        fVals.push_back(pVal);
        //if(std::isnan(pVal) || std::isinf(pVal)) cerr << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt() << " -- " << iConstits[i0].eta() << endl;
        if( ! edm::isFinite(pVal)) {
            LogDebug( "NotFound" )  << "====> Value is Nan " << pVal << " == " << iConstits[i0].pt() << " -- " << iConstits[i0].eta() << endl;
            continue;
        }
        
        // // fPuppiAlgo[pPupId].add(iConstits[i0],pVal,iOpt);
        //code added by Nhan, now instead for every algorithm give it all the particles
        for(int i1 = 0; i1 < fNAlgos; i1++){
            pAlgo    = fPuppiAlgo[i1].algoId   (iOpt);
            pCharged = fPuppiAlgo[i1].isCharged(iOpt);
            pCone    = fPuppiAlgo[i1].coneSize (iOpt);
            double curVal = -1; 
            if(!pCharged) curVal = goodVar(iConstits[i0],iParticles       ,pAlgo,pCone);
            if( pCharged) curVal = goodVar(iConstits[i0],iChargedParticles,pAlgo,pCone);
            //std::cout << "i1 = " << i1 << ", curVal = " << curVal << ", eta = " << iConstits[i0].eta() << ", pupID = " << pPupId << std::endl;
            fPuppiAlgo[i1].add(iConstits[i0],curVal,iOpt);
        }

    }
    for(int i0 = 0; i0 < fNAlgos; i0++) fPuppiAlgo[i0].computeMedRMS(iOpt,fPVFrac);
}

double PuppiContainer::goodVar ( fastjet::PseudoJet const &  iPart, std::vector< fastjet::PseudoJet > const &  iParts, int  iOpt, double  iRCone  )

protected

Definition at line 76 of file PuppiContainer.cc.

                                                                                                                {
    double lPup = 0;
    lPup = var_within_R(iOpt,iParts,iPart,iRCone);
    return lPup;
}

void PuppiContainer::initialize

(

const std::vector< RecoObj > &

iRecoObjects

)

Definition at line 28 of file PuppiContainer.cc.

References funct::abs(), i, edm::isFinite(), and or.

                                                                      {
    //Clear everything
    fRecoParticles.resize(0);
    fPFParticles  .resize(0);
    fChargedPV    .resize(0);
    fPupParticles .resize(0);
    fWeights      .resize(0);
    fVals.resize(0);
    fRawAlphas.resize(0);
    fAlphaMed     .resize(0);
    fAlphaRMS     .resize(0);
    //fChargedNoPV.resize(0);
    //Link to the RecoObjects
    fPVFrac = 0.;
    fNPV    = 1.;
    fRecoParticles = iRecoObjects;
    for (unsigned int i = 0; i < fRecoParticles.size(); i++){
        fastjet::PseudoJet curPseudoJet;
        auto fRecoParticle = fRecoParticles[i];
        // float nom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt)*(cosh(fRecoParticle.eta))*(cosh(fRecoParticle.eta))) + (fRecoParticle.pt)*sinh(fRecoParticle.eta);//hacked
        // float denom = sqrt((fRecoParticle.m)*(fRecoParticle.m) + (fRecoParticle.pt)*(fRecoParticle.pt));//hacked
        // float rapidity = log(nom/denom);//hacked
        if (edm::isFinite(fRecoParticle.rapidity)){
            curPseudoJet.reset_PtYPhiM(fRecoParticle.pt,fRecoParticle.rapidity,fRecoParticle.phi,fRecoParticle.m);//hacked
        } else {        
            curPseudoJet.reset_PtYPhiM(0, 99., 0, 0);//skipping may have been a better choice     
        }                   
        //curPseudoJet.reset_PtYPhiM(fRecoParticle.pt,fRecoParticle.eta,fRecoParticle.phi,fRecoParticle.m);
        int puppi_register = 0;
        if(fRecoParticle.id == 0 or fRecoParticle.charge == 0)  puppi_register = 0; // zero is neutral hadron
        if(fRecoParticle.id == 1 and fRecoParticle.charge != 0) puppi_register = fRecoParticle.charge; // from PV use the
        if(fRecoParticle.id == 2 and fRecoParticle.charge != 0) puppi_register = fRecoParticle.charge+5; // from NPV use the charge as key +5 as key
        curPseudoJet.set_user_info( new PuppiUserInfo( puppi_register ) );
        // fill vector of pseudojets for internal references
        fPFParticles.push_back(curPseudoJet);
        //Take Charged particles associated to PV
        if(std::abs(fRecoParticle.id) == 1) fChargedPV.push_back(curPseudoJet);
        if(std::abs(fRecoParticle.id) >= 1 ) fPVFrac+=1.;
        //if((fRecoParticle.id == 0) && (inParticles[i].id == 2))  _genParticles.push_back( curPseudoJet);
        //if(fRecoParticle.id <= 2 && !(inParticles[i].pt < fNeutralMinE && fRecoParticle.id < 2)) _pfchsParticles.push_back(curPseudoJet);
        //if(fRecoParticle.id == 3) _chargedNoPV.push_back(curPseudoJet);
        // if(fNPV < fRecoParticle.vtxId) fNPV = fRecoParticle.vtxId;
    }
    if (fPVFrac != 0) fPVFrac = double(fChargedPV.size())/fPVFrac;
    else fPVFrac = 0;
}

std::vector<fastjet::PseudoJet> const& PuppiContainer::pfParticles ( ) const

inline

Definition at line 42 of file PuppiContainer.h.

References fPFParticles.

{ return fPFParticles; }    

const std::vector<double>& PuppiContainer::puppiAlphas ( )

inline

Definition at line 46 of file PuppiContainer.h.

References fVals.

{ return fVals; }

const std::vector<double>& PuppiContainer::puppiAlphasMed ( )

inline

Definition at line 48 of file PuppiContainer.h.

References fAlphaMed.

{return fAlphaMed;}

const std::vector<double>& PuppiContainer::puppiAlphasRMS ( )

inline

Definition at line 49 of file PuppiContainer.h.

References fAlphaRMS.

{return fAlphaRMS;}

int PuppiContainer::puppiNAlgos ( )

inline

Definition at line 51 of file PuppiContainer.h.

References fNAlgos.

{ return fNAlgos; }

std::vector<fastjet::PseudoJet> const& PuppiContainer::puppiParticles ( ) const

inline

Definition at line 52 of file PuppiContainer.h.

References fPupParticles.

{ return fPupParticles;}

const std::vector<double>& PuppiContainer::puppiRawAlphas ( )

inline

Definition at line 45 of file PuppiContainer.h.

References fRawAlphas.

{ return fRawAlphas; }

std::vector< double > const & PuppiContainer::puppiWeights

(

)

Definition at line 198 of file PuppiContainer.cc.

References alignCSCRings::e, eta, edm::isFinite(), LogDebug, bookConverter::max, EnergyCorrector::pt, reset(), and plotscripts::rms().

                                                       {
    fPupParticles .resize(0);
    fWeights      .resize(0);
    fVals         .resize(0);
    for(int i0 = 0; i0 < fNAlgos; i0++) fPuppiAlgo[i0].reset();
    
    int lNMaxAlgo = 1;
    for(int i0 = 0; i0 < fNAlgos; i0++) lNMaxAlgo = std::max(fPuppiAlgo[i0].numAlgos(),lNMaxAlgo);
    //Run through all compute mean and RMS
    int lNParticles    = fRecoParticles.size();
    for(int i0 = 0; i0 < lNMaxAlgo; i0++) {
        getRMSAvg(i0,fPFParticles,fPFParticles,fChargedPV);
    }
    if (fPuppiDiagnostics) getRawAlphas(0,fPFParticles,fPFParticles,fChargedPV);

    std::vector<double> pVals;
    for(int i0 = 0; i0 < lNParticles; i0++) {
        //Refresh
        pVals.clear();
        double pWeight = 1;
        //Get the Puppi Id and if ill defined move on
        int  pPupId   = getPuppiId(fRecoParticles[i0].pt,fRecoParticles[i0].eta);
        if(pPupId == -1) {
            fWeights .push_back(pWeight);
            fAlphaMed.push_back(-10);
            fAlphaRMS.push_back(-10);            
            continue;
        }
        // fill the p-values
        double pChi2   = 0;
        if(fUseExp){
            //Compute an Experimental Puppi Weight with delta Z info (very simple example)
            pChi2 = getChi2FromdZ(fRecoParticles[i0].dZ);
            //Now make sure Neutrals are not set
            if(fRecoParticles[i0].pfType > 3) pChi2 = 0;
        }
        //Fill and compute the PuppiWeight
        int lNAlgos = fPuppiAlgo[pPupId].numAlgos();
        for(int i1 = 0; i1 < lNAlgos; i1++) pVals.push_back(fVals[lNParticles*i1+i0]);

        pWeight = fPuppiAlgo[pPupId].compute(pVals,pChi2);
        //Apply the CHS weights
        if(fRecoParticles[i0].id == 1 && fApplyCHS ) pWeight = 1;
        if(fRecoParticles[i0].id == 2 && fApplyCHS ) pWeight = 0;
        //Basic Weight Checks
        if( ! edm::isFinite(pWeight)) {
            pWeight = 0.0;
            LogDebug("PuppiWeightError") << "====> Weight is nan : " << pWeight << " : pt " << fRecoParticles[i0].pt << " -- eta : " << fRecoParticles[i0].eta << " -- Value" << fVals[i0] << " -- id :  " << fRecoParticles[i0].id << " --  NAlgos: " << lNAlgos << std::endl;
        }
        //Basic Cuts
        if(pWeight                         < fPuppiWeightCut) pWeight = 0;  //==> Elminate the low Weight stuff
        if(pWeight*fPFParticles[i0].pt()   < fPuppiAlgo[pPupId].neutralPt(fNPV) && fRecoParticles[i0].id == 0 ) pWeight = 0;  //threshold cut on the neutral Pt
        if(fInvert) pWeight = 1.-pWeight;
        //std::cout << "fRecoParticles[i0].pt = " <<  fRecoParticles[i0].pt << ", fRecoParticles[i0].charge = " << fRecoParticles[i0].charge << ", fRecoParticles[i0].id = " << fRecoParticles[i0].id << ", weight = " << pWeight << std::endl;

        fWeights .push_back(pWeight);
        fAlphaMed.push_back(fPuppiAlgo[pPupId].median());
        fAlphaRMS.push_back(fPuppiAlgo[pPupId].rms());        
        //Now get rid of the thrown out weights for the particle collection

        // leave these lines in, in case want to move eventually to having no 1-to-1 correspondence between puppi and pf cands
        // if( std::abs(pWeight) < std::numeric_limits<double>::denorm_min() ) continue; // this line seems not to work like it's supposed to...
        // if(std::abs(pWeight) <= 0. ) continue; 
        
        //Produce
        PseudoJet curjet( pWeight*fPFParticles[i0].px(), pWeight*fPFParticles[i0].py(), pWeight*fPFParticles[i0].pz(), pWeight*fPFParticles[i0].e() );
        curjet.set_user_index(i0);
        fPupParticles.push_back(curjet);
    }
    return fWeights;
}

std::vector<fastjet::PseudoJet> const& PuppiContainer::pvParticles ( ) const

inline

Definition at line 43 of file PuppiContainer.h.

References fChargedPV.

{ return fChargedPV; }        

void PuppiContainer::setNPV ( int  iNPV )

inline

Definition at line 40 of file PuppiContainer.h.

References fNPV.

{ fNPV = iNPV; }

double PuppiContainer::var_within_R ( int  iId, const std::vector< fastjet::PseudoJet > &  particles, const fastjet::PseudoJet &  centre, double  R  )

protected

Definition at line 81 of file PuppiContainer.cc.

References funct::abs(), i, cmsBatch::log, M_PI, phi(), EgammaValidation_Wenu_cff::sel, and MetTreeProducer::var().

                                                                                                                  {
    if(iId == -1) return 1;
    fastjet::Selector sel = fastjet::SelectorCircle(R);
    sel.set_reference(centre);
    vector<PseudoJet> near_particles = sel(particles);
    double var = 0;
    //double lSumPt = 0;
    //if(iId == 1) for(unsigned int i=0; i<near_particles.size(); i++) lSumPt += near_particles[i].pt();
    for(unsigned int i=0; i<near_particles.size(); i++){
        double pDEta = near_particles[i].eta()-centre.eta();
        double pDPhi = std::abs(near_particles[i].phi()-centre.phi());
        if(pDPhi > 2.*M_PI-pDPhi) pDPhi =  2.*M_PI-pDPhi;
        double pDR2 = pDEta*pDEta+pDPhi*pDPhi;
        if(std::abs(pDR2)  <  0.0001) continue;
        if(iId == 0) var += (near_particles[i].pt()/pDR2);
        if(iId == 1) var += near_particles[i].pt();
        if(iId == 2) var += (1./pDR2);
        if(iId == 3) var += (1./pDR2);
        if(iId == 4) var += near_particles[i].pt();
        if(iId == 5) var += (near_particles[i].pt() * near_particles[i].pt()/pDR2);
    }
    if(iId == 1) var += centre.pt(); //Sum in a cone
    if(iId == 0 && var != 0) var = log(var);
    if(iId == 3 && var != 0) var = log(var);
    if(iId == 5 && var != 0) var = log(var);
    return var;
}

Member Data Documentation

std::vector<double> PuppiContainer::fAlphaMed

protected

Definition at line 70 of file PuppiContainer.h.

Referenced by puppiAlphasMed().

std::vector<double> PuppiContainer::fAlphaRMS

protected

Definition at line 71 of file PuppiContainer.h.

Referenced by puppiAlphasRMS().

bool PuppiContainer::fApplyCHS

protected

Definition at line 73 of file PuppiContainer.h.

std::vector<fastjet::PseudoJet> PuppiContainer::fChargedPV

protected

Definition at line 65 of file PuppiContainer.h.

Referenced by pvParticles().

bool PuppiContainer::fInvert

protected

Definition at line 74 of file PuppiContainer.h.

int PuppiContainer::fNAlgos

protected

Definition at line 79 of file PuppiContainer.h.

Referenced by puppiNAlgos().

double PuppiContainer::fNeutralMinPt

protected

Definition at line 76 of file PuppiContainer.h.

double PuppiContainer::fNeutralSlope

protected

Definition at line 77 of file PuppiContainer.h.

int PuppiContainer::fNPV

protected

Definition at line 80 of file PuppiContainer.h.

Referenced by setNPV().

std::vector<fastjet::PseudoJet> PuppiContainer::fPFParticles

protected

Definition at line 64 of file PuppiContainer.h.

Referenced by pfParticles().

std::vector<fastjet::PseudoJet> PuppiContainer::fPupParticles

protected

Definition at line 66 of file PuppiContainer.h.

Referenced by puppiParticles().

std::vector<PuppiAlgo> PuppiContainer::fPuppiAlgo

protected

Definition at line 82 of file PuppiContainer.h.

bool PuppiContainer::fPuppiDiagnostics

protected

Definition at line 62 of file PuppiContainer.h.

double PuppiContainer::fPuppiWeightCut

protected

Definition at line 78 of file PuppiContainer.h.

double PuppiContainer::fPVFrac

protected

Definition at line 81 of file PuppiContainer.h.

std::vector<double> PuppiContainer::fRawAlphas

protected

Definition at line 69 of file PuppiContainer.h.

Referenced by puppiRawAlphas().

std::vector<RecoObj> PuppiContainer::fRecoParticles

protected

Definition at line 63 of file PuppiContainer.h.

bool PuppiContainer::fUseExp

protected

Definition at line 75 of file PuppiContainer.h.

std::vector<double> PuppiContainer::fVals

protected

Definition at line 68 of file PuppiContainer.h.

Referenced by puppiAlphas().

std::vector<double> PuppiContainer::fWeights

protected

Definition at line 67 of file PuppiContainer.h.