MVAJetPuId.cc

Go to the documentation of this file.

#include "RecoJets/JetProducers/interface/MVAJetPuId.h"
#include "DataFormats/JetReco/interface/PFJet.h"
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "DataFormats/JetReco/interface/Jet.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
#include "DataFormats/Math/interface/deltaR.h"
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
#include "TMatrixDSym.h"
#include "TMatrixDSymEigen.h"
#include "DataFormats/JetReco/interface/PileupJetIdentifier.h"

const float large_val = std::numeric_limits<float>::max();


MVAJetPuId::MVAJetPuId(const edm::ParameterSet & ps) 
{
    impactParTkThreshod_ = 1.;

    tmvaWeights_         = edm::FileInPath(ps.getParameter<std::string>("tmvaWeights")).fullPath(); 
    tmvaMethod_          = ps.getParameter<std::string>("tmvaMethod");
    tmvaVariables_       = ps.getParameter<std::vector<std::string> >("tmvaVariables");
    tmvaSpectators_      = ps.getParameter<std::vector<std::string> >("tmvaSpectators");
    version_             = ps.getParameter<int>("version");
    reader_              = nullptr;
    edm::ParameterSet jetConfig = ps.getParameter<edm::ParameterSet>("JetIdParams");
    for(int i0 = 0; i0 < NWPs; i0++) { 
        std::string lCutType                            = "Tight";
        if(i0 == PileupJetIdentifier::kMedium) lCutType = "Medium";
        if(i0 == PileupJetIdentifier::kLoose)  lCutType = "Loose";
        for(int i1 = 0; i1 < 1; i1++) {
            std::vector<double> pt010  = jetConfig.getParameter<std::vector<double> >(("Pt010_" +lCutType).c_str());
            std::vector<double> pt1020 = jetConfig.getParameter<std::vector<double> >(("Pt1020_"+lCutType).c_str());
            std::vector<double> pt2030 = jetConfig.getParameter<std::vector<double> >(("Pt2030_"+lCutType).c_str());
            std::vector<double> pt3050 = jetConfig.getParameter<std::vector<double> >(("Pt3050_"+lCutType).c_str());
            for(int i2 = 0; i2 < NPts; i2++) mvacut_[i0][0][i2] = pt010 [i2];
            for(int i2 = 0; i2 < NPts; i2++) mvacut_[i0][1][i2] = pt1020[i2];
            for(int i2 = 0; i2 < NPts; i2++) mvacut_[i0][2][i2] = pt2030[i2];
            for(int i2 = 0; i2 < NPts; i2++) mvacut_[i0][3][i2] = pt3050[i2];

        }
    }
    setup();
}





MVAJetPuId::MVAJetPuId(int version,
        const std::string & tmvaWeights, 
        const std::string & tmvaMethod, 
        Float_t impactParTkThreshod,
        const std::vector<std::string> & tmvaVariables
        ) 
{
    impactParTkThreshod_ = impactParTkThreshod;
    tmvaWeights_         = tmvaWeights;
    tmvaMethod_          = tmvaMethod;
    tmvaVariables_       = tmvaVariables;
    version_             = version;

    reader_              = nullptr;

    setup();
}



void MVAJetPuId::setup()
{
    initVariables();


    tmvaVariables_.clear();
    tmvaVariables_.push_back( "rho"  );
    tmvaVariables_.push_back( "nParticles" );
    tmvaVariables_.push_back( "nCharged" );
    tmvaVariables_.push_back( "majW"   );
    tmvaVariables_.push_back( "minW"   );
    tmvaVariables_.push_back( "frac01"    );
    tmvaVariables_.push_back( "frac02" );
    tmvaVariables_.push_back( "frac03"   );
    tmvaVariables_.push_back( "frac04"  );
    tmvaVariables_.push_back( "ptD"  );
    tmvaVariables_.push_back( "beta"   );
    tmvaVariables_.push_back( "betaStar"  );
    tmvaVariables_.push_back( "dR2Mean"  );
    tmvaVariables_.push_back( "pull" );
    tmvaVariables_.push_back( "jetR");
    tmvaVariables_.push_back( "jetRchg");

    tmvaNames_["rho"] = "rho";
    tmvaNames_["nParticles"] = "nParticles";
    tmvaNames_["nCharged"] = "nCharged";
    tmvaNames_["majW"] = "majW";
    tmvaNames_["minW"] = "minW";
    tmvaNames_["frac01"] = "frac01";
    tmvaNames_["frac02"] = "frac02";
    tmvaNames_["frac03"] = "frac03";
    tmvaNames_["frac04"] = "frac04";
    tmvaNames_["ptD"] = "ptD";
    tmvaNames_["beta"] = "beta";
    tmvaNames_["betaStar"] = "betaStar";
    tmvaNames_["dR2Mean"] = "dR2Mean";
    tmvaNames_["pull"] = "pull";  
    tmvaNames_["jetR"] = "jetR";  
    tmvaNames_["jetRchg"] = "jetRchg";
}



MVAJetPuId::~MVAJetPuId() 
{
    if( reader_ ) {
        delete reader_;
    }
}


void Assign(const std::vector<float> & vec, float & a, float & b, float & c, float & d )
{
    size_t sz = vec.size();
    a = ( sz > 0 ? vec[0] : 0. );
    b = ( sz > 1 ? vec[1] : 0. );
    c = ( sz > 2 ? vec[2] : 0. );
    d = ( sz > 3 ? vec[3] : 0. );
}


void SetPtEtaPhi(const reco::Candidate & p, float & pt, float & eta, float &phi )
{
    pt  = p.pt();
    eta = p.eta();
    phi = p.phi();
}


void MVAJetPuId::bookReader()
{
    reader_ = new TMVA::Reader("!Color:Silent");
    assert( ! tmvaMethod_.empty() && !  tmvaWeights_.empty() );
    for(std::vector<std::string>::iterator it=tmvaVariables_.begin(); it!=tmvaVariables_.end(); ++it) {
        if(  tmvaNames_[*it].empty() ) { 
            tmvaNames_[*it] = *it;
        }
        reader_->AddVariable( *it, variables_[ tmvaNames_[*it] ].first );
    }
    for(std::vector<std::string>::iterator it=tmvaSpectators_.begin(); it!=tmvaSpectators_.end(); ++it) {
        if(  tmvaNames_[*it].empty() ) { 
            tmvaNames_[*it] = *it;
        }
        reader_->AddSpectator( *it, variables_[ tmvaNames_[*it] ].first );
    }
    reco::details::loadTMVAWeights(reader_,  tmvaMethod_, tmvaWeights_ ); 
}


void MVAJetPuId::set(const PileupJetIdentifier & id)
{
    internalId_ = id;
}


void MVAJetPuId::runMva()
{
    if( ! reader_ ) { bookReader();}
    if(fabs(internalId_.jetEta_) <  5.0) internalId_.mva_ = reader_->EvaluateMVA( tmvaMethod_.c_str() );
    if(fabs(internalId_.jetEta_) >= 5.0) internalId_.mva_ = -2.;
    internalId_.idFlag_ = computeIDflag(internalId_.mva_,internalId_.jetPt_,internalId_.jetEta_);
}

std::pair<int,int> MVAJetPuId::getJetIdKey(float jetPt, float jetEta)
{
    int ptId = 0;                                                                                                                                    
    if(jetPt > 10 && jetPt < 20) ptId = 1;                                                                                 
    if(jetPt >= 20 && jetPt < 30) ptId = 2;                                                                                 
    if(jetPt >= 30              ) ptId = 3;                                                                                          

    int etaId = 0;                                                                                                                                   
    if(fabs(jetEta) > 2.5  && fabs(jetEta) < 2.75) etaId = 1;                                                              
    if(fabs(jetEta) >= 2.75 && fabs(jetEta) < 3.0 ) etaId = 2;                                                              
    if(fabs(jetEta) >= 3.0  && fabs(jetEta) < 5.0 ) etaId = 3;                                 
    return std::pair<int,int>(ptId,etaId);
}


int MVAJetPuId::computeIDflag(float mva, float jetPt, float jetEta)
{
    std::pair<int,int> jetIdKey = getJetIdKey(jetPt,jetEta);
    return computeIDflag(mva,jetIdKey.first,jetIdKey.second);
}

int MVAJetPuId::computeIDflag(float mva,int ptId,int etaId)
{
    int idFlag(0);
    if(mva > mvacut_[PileupJetIdentifier::kTight ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kTight;
    if(mva > mvacut_[PileupJetIdentifier::kMedium][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kMedium;
    if(mva > mvacut_[PileupJetIdentifier::kLoose ][ptId][etaId]) idFlag += 1 << PileupJetIdentifier::kLoose;
    return idFlag;
}

PileupJetIdentifier MVAJetPuId::computeMva()
{
    runMva();
    return PileupJetIdentifier(internalId_);
}

PileupJetIdentifier MVAJetPuId::computeIdVariables(const reco::Jet * jet, float jec, const reco::Vertex * vtx,
        const reco::VertexCollection & allvtx, double rho,
        bool calculateMva) 
{

    typedef std::vector <reco::PFCandidatePtr> constituents_type;
    typedef std::vector <reco::PFCandidatePtr>::iterator constituents_iterator;

    resetVariables();

    const reco::PFJet * pfjet = dynamic_cast<const reco::PFJet *>(jet);

    if( jec < 0. ) {
        jec = 1.;
    }

    constituents_type constituents = pfjet->getPFConstituents();

    reco::PFCandidatePtr lLead, lSecond, lLeadNeut, lLeadEm, lLeadCh, lTrail;
    std::vector<float> frac, fracCh, fracEm, fracNeut;
    constexpr int ncones = 4;
    std::array<float, ncones> cones{ {0.1,0.2,0.3,0.4}};
    std::array<float *, ncones> coneFracs{ {&internalId_.frac01_,&internalId_.frac02_,&internalId_.frac03_,&internalId_.frac04_}};
    TMatrixDSym covMatrix(2); covMatrix = 0.;

    reco::TrackRef impactTrack;
    float jetPt = jet->pt() / jec; // use uncorrected pt for shape variables
    float sumPt = 0., sumPt2 = 0., sumTkPt = 0.,sumPtCh=0,sumPtNe = 0; float sum_deta =0 ; float sum_dphi =0 ; float sum_deta2 =0 ; float sum_detadphi =0 ; float sum_dphi2=0;
    SetPtEtaPhi(*jet,internalId_.jetPt_,internalId_.jetEta_,internalId_.jetPhi_); // use corrected pt for jet kinematics
    internalId_.jetM_ = jet->mass(); 
    internalId_.rho_ = rho;//allvtx.size();
    for(constituents_iterator it=constituents.begin(); it!=constituents.end(); ++it) {
        reco::PFCandidatePtr & icand = *it;
        float candPt = icand->pt();
        float candPtFrac = candPt/jetPt;
        float candDr   = reco::deltaR(**it,*jet);
        float candDeta = fabs( (*it)->eta() - jet->eta() );
        float candDphi = reco::deltaPhi(**it,*jet);
        float candPtDr = candPt * candDr;
        size_t icone = std::lower_bound(&cones[0],&cones[ncones],candDr) - &cones[0];
        float weight2 = candPt * candPt;


        if( lLead.isNull() || candPt > lLead->pt() ) {
            lSecond = lLead;
            lLead = icand; 
        } else if( (lSecond.isNull() || candPt > lSecond->pt()) && (candPt < lLead->pt()) ) {
            lSecond = icand;
        }

        //internalId_.dRMean_     += candPtDr;
        internalId_.dR2Mean_    += candPtDr*candPtDr;

        internalId_.ptD_ += candPt*candPt;
        sumPt += candPt;
        sumPt2 += candPt*candPt;
        sum_deta     += candDeta*weight2;
        sum_dphi     += candDphi*weight2;
        sum_deta2    += candDeta*candDeta*weight2;
        sum_detadphi += candDeta*candDphi*weight2;
        sum_dphi2    += candDphi*candDphi*weight2;
        //Teta         += candPt * candDR * candDeta;
        //Tphi         += candPt * candDR * candDphi;



        frac.push_back(candPtFrac);
        if( icone < ncones ) { *coneFracs[icone] += candPt; }


        if( icand->particleId() == reco::PFCandidate::h0 ) {
            if (lLeadNeut.isNull() || candPt > lLeadNeut->pt()) { lLeadNeut = icand; }
            internalId_.dRMeanNeut_ += candPtDr;
            fracNeut.push_back(candPtFrac);
            sumPtNe               += candPt;
        }

        if( icand->particleId() == reco::PFCandidate::gamma ) {
            if(lLeadEm.isNull() || candPt > lLeadEm->pt())  { lLeadEm = icand; }
            internalId_.dRMeanEm_ += candPtDr;
            fracEm.push_back(candPtFrac);
            sumPtNe               += candPt;
        }

        if(  icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
            if (lLeadCh.isNull() || candPt > lLeadCh->pt()) { lLeadCh = icand; }
            //internalId_.jetRchg_  += candPtDr;
            fracCh.push_back(candPtFrac);
            sumPtCh                += candPt;
        }

        if(  icand->trackRef().isNonnull() && icand->trackRef().isAvailable() ) {
            float tkpt = icand->trackRef()->pt(); 
            sumTkPt += tkpt;
            bool inVtx0 = find( vtx->tracks_begin(), vtx->tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != vtx->tracks_end();
            bool inAnyOther = false;

            double dZ0 = fabs(icand->trackRef()->dz(vtx->position()));
            double dZ = dZ0;
            for(reco::VertexCollection::const_iterator  vi=allvtx.begin(); vi!=allvtx.end(); ++vi ) {
                const reco::Vertex & iv = *vi;
                if( iv.isFake() || iv.ndof() < 4 ) { continue; }

                bool isVtx0  = (iv.position() - vtx->position()).r() < 0.02;

                if( ! isVtx0 && ! inAnyOther ) {
                    inAnyOther = find( iv.tracks_begin(), iv.tracks_end(), reco::TrackBaseRef(icand->trackRef()) ) != iv.tracks_end();
                }

                dZ = std::min(dZ,fabs(icand->trackRef()->dz(iv.position())));
            }
            if( inVtx0 && ! inAnyOther ) {
                internalId_.betaClassic_ += tkpt;
            } else if( ! inVtx0 && inAnyOther ) {
                internalId_.betaStarClassic_ += tkpt;
            }

            if( dZ0 < 0.2 ) {
                internalId_.beta_ += tkpt;
            } else if( dZ < 0.2 ) {
                internalId_.betaStar_ += tkpt;
            }
        } 

        if( lTrail.isNull() || candPt < lTrail->pt() ) {
            lTrail = icand; 
        }
    }

    assert( lLead.isNonnull() );
    if ( lSecond.isNull() )   { lSecond   = lTrail; }
    if ( lLeadNeut.isNull() ) { lLeadNeut = lTrail; }
    if ( lLeadEm.isNull() )   { lLeadEm   = lTrail; }
    if ( lLeadCh.isNull() )   { lLeadCh   = lTrail; }
    impactTrack = lLeadCh->trackRef();

    internalId_.nCharged_    = pfjet->chargedMultiplicity();
    internalId_.nNeutrals_   = pfjet->neutralMultiplicity();
    internalId_.chgEMfrac_   = pfjet->chargedEmEnergy()    /jet->energy();
    internalId_.neuEMfrac_   = pfjet->neutralEmEnergy()    /jet->energy();
    internalId_.chgHadrfrac_ = pfjet->chargedHadronEnergy()/jet->energy();
    internalId_.neuHadrfrac_ = pfjet->neutralHadronEnergy()/jet->energy();

    if( impactTrack.isNonnull() && impactTrack.isAvailable() ) {
        internalId_.d0_ = fabs(impactTrack->dxy(vtx->position()));
        internalId_.dZ_ = fabs(impactTrack->dz(vtx->position()));
    }else{ 
        internalId_.nParticles_ = constituents.size(); 
        SetPtEtaPhi(*lLead,internalId_.leadPt_,internalId_.leadEta_,internalId_.leadPhi_);                 
        SetPtEtaPhi(*lSecond,internalId_.secondPt_,internalId_.secondEta_,internalId_.secondPhi_);        
        SetPtEtaPhi(*lLeadNeut,internalId_.leadNeutPt_,internalId_.leadNeutEta_,internalId_.leadNeutPhi_); 
        SetPtEtaPhi(*lLeadEm,internalId_.leadEmPt_,internalId_.leadEmEta_,internalId_.leadEmPhi_);        
        SetPtEtaPhi(*lLeadCh,internalId_.leadChPt_,internalId_.leadChEta_,internalId_.leadChPhi_);         
        std::sort(frac.begin(),frac.end(),std::greater<float>());
        std::sort(fracCh.begin(),fracCh.end(),std::greater<float>());
        std::sort(fracEm.begin(),fracEm.end(),std::greater<float>());
        std::sort(fracNeut.begin(),fracNeut.end(),std::greater<float>());
        Assign(frac,    internalId_.leadFrac_,    internalId_.secondFrac_,    internalId_.thirdFrac_,    internalId_.fourthFrac_);

        //covMatrix(0,0) /= sumPt2;
        //covMatrix(0,1) /= sumPt2;
        //covMatrix(1,1) /= sumPt2;
        //covMatrix(1,0)  = covMatrix(0,1);
        //internalId_.etaW_ = sqrt(covMatrix(0,0));
        //internalId_.phiW_ = sqrt(covMatrix(1,1));
        //internalId_.jetW_ = 0.5*(internalId_.etaW_+internalId_.phiW_);
        //TVectorD eigVals(2); eigVals = TMatrixDSymEigen(covMatrix).GetEigenValues();
        //  
        if( internalId_.majW_ < internalId_.minW_ ) { std::swap(internalId_.majW_,internalId_.minW_); }

        //internalId_.dRLeadCent_ = reco::deltaR(*jet,*lLead);
        if( lSecond.isNonnull() ) { internalId_.dRLead2nd_  = reco::deltaR(*jet,*lSecond); }
        internalId_.dRMeanNeut_ /= jetPt;
        internalId_.dRMeanEm_   /= jetPt;
        //internalId_.jetRchg_   /= jetPt;
        internalId_.dR2Mean_    /= sumPt2;
        for(size_t ic=0; ic<ncones; ++ic){
            *coneFracs[ic]     /= jetPt;
        }

        double ptMean = sumPt/internalId_.nParticles_;
        double ptRMS  = 0;
        for(unsigned int i0 = 0; i0 < frac.size(); i0++) {ptRMS+=(frac[i0]-ptMean)*(frac[i0]-ptMean);}
        ptRMS/=internalId_.nParticles_;
        ptRMS=sqrt(ptRMS);
        internalId_.jetRchg_ = internalId_.leadChPt_/sumPt;
        internalId_.jetR_ = internalId_.leadPt_/sumPt;

        internalId_.ptMean_  = ptMean;
        internalId_.ptRMS_   = ptRMS/jetPt;
        internalId_.pt2A_    = sqrt( internalId_.ptD_     /internalId_.nParticles_)/jetPt;
        internalId_.ptD_     = sqrt( internalId_.ptD_)    / sumPt;
        internalId_.sumPt_   = sumPt;
        internalId_.sumChPt_ = sumPtCh;
        internalId_.sumNePt_ = sumPtNe;
        if (sumPt > 0) {
            internalId_.beta_     /= sumPt;
            internalId_.betaStar_ /= sumPt;
        } else {
            assert( internalId_.beta_ == 0. && internalId_.betaStar_ == 0.);
        }
        float ave_deta = sum_deta/sumPt2;
        float ave_dphi = sum_dphi/sumPt2;
        float ave_deta2 = sum_deta2/sumPt2;
        float ave_dphi2 = sum_dphi2/sumPt2;
        float a = ave_deta2-ave_deta*ave_deta;
        float b = ave_dphi2-ave_dphi*ave_dphi;
        float c = -(sum_detadphi/sumPt2-ave_deta*ave_dphi);
        float axis1=0; float axis2=0;
        if((((a-b)*(a-b)+4*c*c))>0) {
        float delta = sqrt(((a-b)*(a-b)+4*c*c));
        if (a+b+delta > 0) {
            axis1 = sqrt(0.5*(a+b+delta));
        }
        if (a+b-delta > 0) {
            axis2 = sqrt(0.5*(a+b-delta));
        }
        }
        else {
           axis1=-1;
           axis2=-1;
        }
        internalId_.majW_ = axis1; //sqrt(fabs(eigVals(0)));
        internalId_.minW_ = axis2;//sqrt(fabs(eigVals(1)));
        //compute Pull

        float ddetaR_sum(0.0), ddphiR_sum(0.0);
        for(int i=0; i<internalId_.nParticles_; ++i) {
            reco::PFCandidatePtr part = pfjet->getPFConstituent(i);
            float weight =part->pt()*part->pt() ;
            float deta = part->eta() - jet->eta();
            float dphi = reco::deltaPhi(*part, *jet);
            float ddeta, ddphi, ddR;
            ddeta = deta - ave_deta ;
            ddphi = reco::deltaPhi(dphi,ave_dphi);//2*atan(tan((dphi - ave_dphi)/2.)) ;
            ddR = sqrt(ddeta*ddeta + ddphi*ddphi);
            ddetaR_sum += ddR*ddeta*weight;
            ddphiR_sum += ddR*ddphi*weight;
        }if (sumPt2 > 0) {
            float ddetaR_ave = ddetaR_sum/sumPt2;
            float ddphiR_ave = ddphiR_sum/sumPt2;
            internalId_.dRMean_ = sqrt(ddetaR_ave*ddetaR_ave+ddphiR_ave*ddphiR_ave);
        }



    }

    if( calculateMva ) {
        runMva();
    }

    return PileupJetIdentifier(internalId_);
}

std::string MVAJetPuId::dumpVariables() const
{
    std::stringstream out;
    for(variables_list_t::const_iterator it=variables_.begin(); 
            it!=variables_.end(); ++it ) {
        out << std::setw(15) << it->first << std::setw(3) << "=" 
            << std::setw(5) << *it->second.first 
            << " (" << std::setw(5) << it->second.second << ")" << std::endl;
    }
    return out.str();
}

void MVAJetPuId::resetVariables()
{
    internalId_.idFlag_    = 0;
    for(variables_list_t::iterator it=variables_.begin(); 
            it!=variables_.end(); ++it ) {
        *it->second.first = it->second.second;
    }
}

#define INIT_VARIABLE(NAME,TMVANAME,VAL)    \
    internalId_.NAME ## _ = VAL; \
variables_[ # NAME   ] = std::make_pair(& internalId_.NAME ## _, VAL);

void MVAJetPuId::initVariables()
{
    internalId_.idFlag_    = 0;
    INIT_VARIABLE(mva        , "", -100.);

    INIT_VARIABLE(jetPt      , "jetPt", 0.);
    INIT_VARIABLE(jetEta     , "jetEta", large_val);
    INIT_VARIABLE(jetPhi     , "", large_val);
    INIT_VARIABLE(jetM       , "", 0.);
    INIT_VARIABLE(nCharged   , "nCharged", 0.);
    INIT_VARIABLE(nNeutrals  , "", 0.);

    INIT_VARIABLE(chgEMfrac  , "", 0.);
    INIT_VARIABLE(neuEMfrac  , "", 0.);
    INIT_VARIABLE(chgHadrfrac, "", 0.);
    INIT_VARIABLE(neuHadrfrac, "", 0.);

    INIT_VARIABLE(d0         , ""    , -1000.);   
    INIT_VARIABLE(dZ         , ""    , -1000.);  
    INIT_VARIABLE(nParticles , "nParticles"  , 0.);  

    INIT_VARIABLE(leadPt     , ""    , 0.);  
    INIT_VARIABLE(leadEta    , ""   , large_val);  
    INIT_VARIABLE(leadPhi    , ""   , large_val);  
    INIT_VARIABLE(secondPt   , ""    , 0.);  
    INIT_VARIABLE(secondEta  , ""   , large_val);  
    INIT_VARIABLE(secondPhi  , ""   , large_val);  
    INIT_VARIABLE(leadNeutPt , ""    , 0.);  
    INIT_VARIABLE(leadNeutEta, ""   , large_val);  

    INIT_VARIABLE(jetR , "jetR"   , 0.);  
    INIT_VARIABLE(pull     , "pull"    , 0.);  
    INIT_VARIABLE(jetRchg   , "jetRchg"   , 0.);  
    INIT_VARIABLE(dR2Mean    , "dR2Mean"         , 0.);  

    INIT_VARIABLE(ptD        , "ptD", 0.);
    INIT_VARIABLE(ptMean     , "", 0.);
    INIT_VARIABLE(ptRMS      , "", 0.);
    INIT_VARIABLE(pt2A       , "", 0.);
    INIT_VARIABLE(ptDCh      , "", 0.);
    INIT_VARIABLE(ptDNe      , "", 0.);
    INIT_VARIABLE(sumPt      , "", 0.);
    INIT_VARIABLE(sumChPt    , "", 0.);
    INIT_VARIABLE(sumNePt    , "", 0.);
    INIT_VARIABLE(secondFrac  ,"" ,0.);  
    INIT_VARIABLE(thirdFrac   ,"" ,0.);  
    INIT_VARIABLE(fourthFrac  ,"" ,0.);  
    INIT_VARIABLE(leadChFrac    ,"" ,0.);  
    INIT_VARIABLE(secondChFrac  ,"" ,0.);  
    INIT_VARIABLE(thirdChFrac   ,"" ,0.);  
    INIT_VARIABLE(fourthChFrac  ,"" ,0.);  
    INIT_VARIABLE(leadNeutFrac    ,"" ,0.);  
    INIT_VARIABLE(secondNeutFrac  ,"" ,0.);  
    INIT_VARIABLE(thirdNeutFrac   ,"" ,0.);  
    INIT_VARIABLE(fourthNeutFrac  ,"" ,0.);  
    INIT_VARIABLE(leadEmFrac    ,"" ,0.);  
    INIT_VARIABLE(secondEmFrac  ,"" ,0.);  
    INIT_VARIABLE(thirdEmFrac   ,"" ,0.);  
    INIT_VARIABLE(fourthEmFrac  ,"" ,0.);  
    INIT_VARIABLE(jetW  ,"" ,1.);  
    INIT_VARIABLE(etaW  ,"" ,1.);  
    INIT_VARIABLE(phiW  ,"" ,1.);  
    INIT_VARIABLE(majW  ,"majW" ,1.);  
    INIT_VARIABLE(minW  ,"minW" ,1.);  
    INIT_VARIABLE(frac01    ,"frac01" ,0.);  
    INIT_VARIABLE(frac02    ,"frac02" ,0.);  
    INIT_VARIABLE(frac03    ,"frac03" ,0.);  
    INIT_VARIABLE(frac04    ,"frac04" ,0.);  

    INIT_VARIABLE(beta   ,"beta" ,0.);  
    INIT_VARIABLE(betaStar   ,"betaStar" ,0.);  
    INIT_VARIABLE(betaClassic   ,"betaClassic" ,0.);  
    INIT_VARIABLE(betaStarClassic   ,"betaStarClassic" ,0.);  
    INIT_VARIABLE(rho   ,"rho" ,0.);  

}
#undef INIT_VARIABLE