Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:36 

0001 #include <fstream>
0002 #include <iomanip>
0003 #include <iostream>
0004 
0005 // DataFormat //
0006 
0007 #include "DataFormats/Common/interface/Handle.h"
0008 #include "DataFormats/Provenance/interface/Provenance.h"
0009 
0010 #include "DataFormats/JetReco/interface/Jet.h"
0011 #include "DataFormats/JetReco/interface/JetTracksAssociation.h"
0012 
0013 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0014 #include "DataFormats/JetReco/interface/CaloJet.h"
0015 
0016 #include "DataFormats/JetReco/interface/JPTJetCollection.h"
0017 #include "DataFormats/JetReco/interface/JPTJet.h"
0018 
0019 #include "DataFormats/JetReco/interface/JetID.h"
0020 
0021 #include "DataFormats/Math/interface/deltaR.h"
0022 #include "DataFormats/Math/interface/LorentzVector.h"
0023 #include "DataFormats/Math/interface/Vector3D.h"
0024 
0025 // FWCore //
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 #include "FWCore/Framework/interface/ESHandle.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 
0032 // TrackingTools //
0033 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0034 
0035 // Constants, Math //
0036 #include "CLHEP/Units/GlobalPhysicalConstants.h"
0037 #include "Math/GenVector/VectorUtil.h"
0038 #include "Math/GenVector/PxPyPzE4D.h"
0039 
0040 #include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
0041 
0042 #include <vector>
0043 
0044 ///////////////////////////////////////////////////////////////////////////////////////
0045 
0046 #include "RecoJets/JetProducers/interface/PileupJPTJetIdAlgo.h"
0047 
0048 ////////////////////////////////////////////////////////////////////////////////////////
0049 using namespace std;
0050 namespace cms {
0051 
0052   PileupJPTJetIdAlgo::PileupJPTJetIdAlgo(const edm::ParameterSet& iConfig) {
0053     verbosity = iConfig.getParameter<int>("Verbosity");
0054     tmvaWeights_ = edm::FileInPath(iConfig.getParameter<std::string>("tmvaWeightsCentral")).fullPath();
0055     tmvaWeightsF_ = edm::FileInPath(iConfig.getParameter<std::string>("tmvaWeightsForward")).fullPath();
0056     tmvaMethod_ = iConfig.getParameter<std::string>("tmvaMethod");
0057   }
0058 
0059   PileupJPTJetIdAlgo::~PileupJPTJetIdAlgo() {
0060     //  std::cout<<" JetPlusTrack destructor "<<std::endl;
0061   }
0062 
0063   void PileupJPTJetIdAlgo::bookMVAReader() {
0064     // Read TMVA tree
0065     //    std::string mvaw     = "RecoJets/JetProducers/data/TMVAClassification_BDTG.weights.xml";
0066     //    std::string mvawF     = "RecoJets/JetProducers/data/TMVAClassification_BDTG.weights_F.xml";
0067     //    tmvaWeights_         = edm::FileInPath(mvaw).fullPath();
0068     //    tmvaWeightsF_         = edm::FileInPath(mvawF).fullPath();
0069     //    tmvaMethod_          = "BDTG method";
0070 
0071     if (verbosity > 0)
0072       std::cout << " TMVA method " << tmvaMethod_.c_str() << " " << tmvaWeights_.c_str() << std::endl;
0073     reader_ = new TMVA::Reader("!Color:!Silent");
0074 
0075     reader_->AddVariable("Nvtx", &Nvtx);
0076     reader_->AddVariable("PtJ", &PtJ);
0077     reader_->AddVariable("EtaJ", &EtaJ);
0078     reader_->AddVariable("Beta", &Beta);
0079     reader_->AddVariable("MultCalo", &MultCalo);
0080     reader_->AddVariable("dAxis1c", &dAxis1c);
0081     reader_->AddVariable("dAxis2c", &dAxis2c);
0082     reader_->AddVariable("MultTr", &MultTr);
0083     reader_->AddVariable("dAxis1t", &dAxis1t);
0084     reader_->AddVariable("dAxis2t", &dAxis2t);
0085 
0086     reco::details::loadTMVAWeights(reader_, tmvaMethod_, tmvaWeights_);
0087 
0088     //std::cout<<" Method booked "<<std::endl;
0089 
0090     readerF_ = new TMVA::Reader("!Color:!Silent");
0091 
0092     readerF_->AddVariable("Nvtx", &Nvtx);
0093     readerF_->AddVariable("PtJ", &PtJ);
0094     readerF_->AddVariable("EtaJ", &EtaJ);
0095     readerF_->AddVariable("MultCalo", &MultCalo);
0096     readerF_->AddVariable("dAxis1c", &dAxis1c);
0097     readerF_->AddVariable("dAxis2c", &dAxis2c);
0098 
0099     reco::details::loadTMVAWeights(readerF_, tmvaMethod_, tmvaWeightsF_);
0100 
0101     // std::cout<<" Method booked F "<<std::endl;
0102   }
0103 
0104   float PileupJPTJetIdAlgo::fillJPTBlock(const reco::JPTJet* jet) {
0105     if (verbosity > 0) {
0106       std::cout << "================================    JET LOOP  ====================  " << std::endl;
0107       std::cout << "================    jetPt   =  " << (*jet).pt() << std::endl;
0108       std::cout << "================    jetEta   =  " << (*jet).eta() << std::endl;
0109     }
0110 
0111     const edm::RefToBase<reco::Jet> jptjetRef = jet->getCaloJetRef();
0112     reco::CaloJet const* rawcalojet = dynamic_cast<reco::CaloJet const*>(&*jptjetRef);
0113 
0114     int ncalotowers = 0.;
0115     double sumpt = 0.;
0116     double dphi2 = 0.;
0117     double deta2 = 0.;
0118     double dphi1 = 0.;
0119     double dphideta = 0.;
0120     double deta1 = 0.;
0121     double EE = 0.;
0122     double HE = 0.;
0123     double EELong = 0.;
0124     double EEShort = 0.;
0125 
0126     std::vector<CaloTowerPtr> calotwrs = (*rawcalojet).getCaloConstituents();
0127 
0128     if (verbosity > 0) {
0129       std::cout << "=======    CaloTowerPtr DONE   ==  " << std::endl;
0130     }
0131 
0132     for (std::vector<CaloTowerPtr>::const_iterator icalot = calotwrs.begin(); icalot != calotwrs.end(); icalot++) {
0133       ncalotowers++;
0134 
0135       double deta = (*jet).eta() - (*icalot)->eta();
0136       double dphi = (*jet).phi() - (*icalot)->phi();
0137 
0138       if (dphi > M_PI)
0139         dphi = dphi - 2. * M_PI;
0140       if (dphi < -1. * M_PI)
0141         dphi = dphi + 2. * M_PI;
0142 
0143       if (verbosity > 0)
0144         std::cout << " CaloTower jet eta " << (*jet).eta() << " tower eta " << (*icalot)->eta() << " jet phi "
0145                   << (*jet).phi() << " tower phi " << (*icalot)->phi() << " dphi " << dphi << " " << (*icalot)->pt()
0146                   << " ieta " << (*icalot)->ieta() << " " << abs((*icalot)->ieta()) << std::endl;
0147 
0148       if (abs((*icalot)->ieta()) < 30)
0149         EE = EE + (*icalot)->emEnergy();
0150       if (abs((*icalot)->ieta()) < 30)
0151         HE = HE + (*icalot)->hadEnergy();
0152       if (abs((*icalot)->ieta()) > 29)
0153         EELong = EELong + (*icalot)->emEnergy();
0154       if (abs((*icalot)->ieta()) > 29)
0155         EEShort = EEShort + (*icalot)->hadEnergy();
0156 
0157       sumpt = sumpt + (*icalot)->pt();
0158       dphi2 = dphi2 + dphi * dphi * (*icalot)->pt();
0159       deta2 = deta2 + deta * deta * (*icalot)->pt();
0160       dphi1 = dphi1 + dphi * (*icalot)->pt();
0161       deta1 = deta1 + deta * (*icalot)->pt();
0162       dphideta = dphideta + dphi * deta * (*icalot)->pt();
0163 
0164     }  // calojet constituents
0165 
0166     if (sumpt > 0.) {
0167       deta1 = deta1 / sumpt;
0168       dphi1 = dphi1 / sumpt;
0169       deta2 = deta2 / sumpt;
0170       dphi2 = dphi2 / sumpt;
0171       dphideta = dphideta / sumpt;
0172     }
0173 
0174     // W.r.t. principal axis
0175 
0176     double detavar = deta2 - deta1 * deta1;
0177     double dphivar = dphi2 - dphi1 * dphi1;
0178     double dphidetacov = dphideta - deta1 * dphi1;
0179 
0180     double det = (detavar - dphivar) * (detavar - dphivar) + 4 * dphidetacov * dphidetacov;
0181     det = sqrt(det);
0182     double x1 = (detavar + dphivar + det) / 2.;
0183     double x2 = (detavar + dphivar - det) / 2.;
0184 
0185     if (verbosity > 0)
0186       std::cout << " ncalo " << ncalotowers << " deta2 " << deta2 << " dphi2 " << dphi2 << " deta1 " << deta1
0187                 << " dphi1 " << dphi1 << " detavar " << detavar << " dphivar " << dphivar << " dphidetacov "
0188                 << dphidetacov << " sqrt(det) " << sqrt(det) << " x1 " << x1 << " x2 " << x2 << std::endl;
0189 
0190     // For jets with |eta|<2 take also tracks shape
0191     int ntracks = 0.;
0192     double sumpttr = 0.;
0193     double dphitr2 = 0.;
0194     double detatr2 = 0.;
0195     double dphitr1 = 0.;
0196     double detatr1 = 0.;
0197     double dphidetatr = 0.;
0198 
0199     const reco::TrackRefVector pioninin = (*jet).getPionsInVertexInCalo();
0200 
0201     for (reco::TrackRefVector::const_iterator it = pioninin.begin(); it != pioninin.end(); it++) {
0202       if ((*it)->pt() > 0.5 && ((*it)->ptError() / (*it)->pt()) < 0.05) {
0203         ntracks++;
0204         sumpttr = sumpttr + (*it)->pt();
0205         double deta = (*jet).eta() - (*it)->eta();
0206         double dphi = (*jet).phi() - (*it)->phi();
0207 
0208         if (dphi > M_PI)
0209           dphi = dphi - 2. * M_PI;
0210         if (dphi < -1. * M_PI)
0211           dphi = dphi + 2. * M_PI;
0212 
0213         dphitr2 = dphitr2 + dphi * dphi * (*it)->pt();
0214         detatr2 = detatr2 + deta * deta * (*it)->pt();
0215         dphitr1 = dphitr1 + dphi * (*it)->pt();
0216         detatr1 = detatr1 + deta * (*it)->pt();
0217         dphidetatr = dphidetatr + dphi * deta * (*it)->pt();
0218         if (verbosity > 0)
0219           std::cout << " Tracks-in-in " << (*it)->pt() << " " << (*it)->eta() << " " << (*it)->phi() << " in jet "
0220                     << (*jet).eta() << " " << (*jet).phi() << " jet pt " << (*jet).pt() << std::endl;
0221       }
0222     }  // pioninin
0223 
0224     const reco::TrackRefVector pioninout = (*jet).getPionsInVertexOutCalo();
0225 
0226     for (reco::TrackRefVector::const_iterator it = pioninout.begin(); it != pioninout.end(); it++) {
0227       if ((*it)->pt() > 0.5 && ((*it)->ptError() / (*it)->pt()) < 0.05) {
0228         ntracks++;
0229         sumpttr = sumpttr + (*it)->pt();
0230         double deta = (*jet).eta() - (*it)->eta();
0231         double dphi = (*jet).phi() - (*it)->phi();
0232 
0233         if (dphi > M_PI)
0234           dphi = dphi - 2. * M_PI;
0235         if (dphi < -1. * M_PI)
0236           dphi = dphi + 2. * M_PI;
0237 
0238         dphitr2 = dphitr2 + dphi * dphi * (*it)->pt();
0239         detatr2 = detatr2 + deta * deta * (*it)->pt();
0240         dphitr1 = dphitr1 + dphi * (*it)->pt();
0241         detatr1 = detatr1 + deta * (*it)->pt();
0242         dphidetatr = dphidetatr + dphi * deta * (*it)->pt();
0243         if (verbosity > 0)
0244           std::cout << " Tracks-in-in " << (*it)->pt() << " " << (*it)->eta() << " " << (*it)->phi() << " in jet "
0245                     << (*jet).eta() << " " << (*jet).phi() << " jet pt " << (*jet).pt() << std::endl;
0246       }
0247     }  // pioninout
0248 
0249     if (verbosity > 0)
0250       std::cout << " Number of tracks in-in and in-out " << pioninin.size() << " " << pioninout.size() << std::endl;
0251 
0252     if (sumpttr > 0.) {
0253       detatr1 = detatr1 / sumpttr;
0254       dphitr1 = dphitr1 / sumpttr;
0255       detatr2 = detatr2 / sumpttr;
0256       dphitr2 = dphitr2 / sumpttr;
0257       dphidetatr = dphidetatr / sumpttr;
0258     }
0259 
0260     // W.r.t. principal axis
0261 
0262     double detavart = detatr2 - detatr1 * detatr1;
0263     double dphivart = dphitr2 - dphitr1 * dphitr1;
0264     double dphidetacovt = dphidetatr - detatr1 * dphitr1;
0265 
0266     double dettr = (detavart - dphivart) * (detavart - dphivart) + 4 * dphidetacovt * dphidetacovt;
0267     dettr = sqrt(dettr);
0268     double x1tr = (detavart + dphivart + dettr) / 2.;
0269     double x2tr = (detavart + dphivart - dettr) / 2.;
0270 
0271     if (verbosity > 0)
0272       std::cout << " ntracks " << ntracks << " detatr2 " << detatr2 << " dphitr2 " << dphitr2 << " detatr1 " << detatr1
0273                 << " dphitr1 " << dphitr1 << " detavart " << detavart << " dphivart " << dphivart << " dphidetacovt "
0274                 << dphidetacovt << " sqrt(det) " << sqrt(dettr) << " x1tr " << x1tr << " x2tr " << x2tr << std::endl;
0275 
0276     // Multivariate analisis
0277     PtJ = (*jet).pt();
0278     EtaJ = (*jet).eta();
0279     Beta = (*jet).getSpecific().Zch;
0280     dAxis2c = x2;
0281     dAxis1c = x1;
0282     dAxis2t = x2tr;
0283     dAxis1t = x1tr;
0284     MultCalo = ncalotowers;
0285     MultTr = ntracks;
0286 
0287     float mva_ = 1.;
0288     if (fabs(EtaJ) < 2.6) {
0289       mva_ = reader_->EvaluateMVA("BDTG method");
0290       // std::cout<<" MVA analysis "<<mva_<<std::endl;
0291     } else {
0292       mva_ = readerF_->EvaluateMVA("BDTG method");
0293       // std::cout<<" MVA analysis Forward "<<mva_<<std::endl;
0294     }
0295     if (verbosity > 0)
0296       std::cout << "=======  Computed MVA =  " << mva_ << std::endl;
0297     return mva_;
0298   }  // fillJPTBlock
0299 }  // namespace cms
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!