HLTanalyzers/plugins/HLTMCtruth.cc

0001 #include <iostream>
0002 #include <sstream>
0003 #include <istream>
0004 #include <fstream>
0005 #include <iomanip>
0006 #include <string>
0007 #include <cmath>
0008 #include <functional>
0009 #include <cstdlib>
0010 #include <cstring>
0011
0012 #include "HLTMCtruth.h"
0013 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0014
0015 HLTMCtruth::HLTMCtruth() {
0016   //set parameter defaults
0017   _Monte = false;
0018   _Gen = false;
0019   _Debug = false;
0020 }
0021
0022 /*  Setup the analysis to put the branch-variables into the tree. */
0023 void HLTMCtruth::setup(const edm::ParameterSet& pSet, TTree* HltTree) {
0024   edm::ParameterSet myMCParams = pSet.getParameter<edm::ParameterSet>("RunParameters");
0025   std::vector<std::string> parameterNames = myMCParams.getParameterNames();
0026
0027   for (auto& parameterName : parameterNames) {
0028     if (parameterName == "Monte")
0029       _Monte = myMCParams.getParameter<bool>(parameterName);
0030     else if (parameterName == "GenTracks")
0031       _Gen = myMCParams.getParameter<bool>(parameterName);
0032     else if (parameterName == "Debug")
0033       _Debug = myMCParams.getParameter<bool>(parameterName);
0034   }
0035
0036   const int kMaxMcTruth = 10000;
0037   mcpid = new int[kMaxMcTruth];
0038   mcstatus = new int[kMaxMcTruth];
0039   mcvx = new float[kMaxMcTruth];
0040   mcvy = new float[kMaxMcTruth];
0041   mcvz = new float[kMaxMcTruth];
0042   mcpt = new float[kMaxMcTruth];
0043   mceta = new float[kMaxMcTruth];
0044   mcphi = new float[kMaxMcTruth];
0045
0046   // MCtruth-specific branches of the tree
0047   HltTree->Branch("NMCpart", &nmcpart, "NMCpart/I");
0048   HltTree->Branch("MCpid", mcpid, "MCpid[NMCpart]/I");
0049   HltTree->Branch("MCstatus", mcstatus, "MCstatus[NMCpart]/I");
0050   HltTree->Branch("MCvtxX", mcvx, "MCvtxX[NMCpart]/F");
0051   HltTree->Branch("MCvtxY", mcvy, "MCvtxY[NMCpart]/F");
0052   HltTree->Branch("MCvtxZ", mcvz, "MCvtxZ[NMCpart]/F");
0053   HltTree->Branch("MCpt", mcpt, "MCpt[NMCpart]/F");
0054   HltTree->Branch("MCeta", mceta, "MCeta[NMCpart]/F");
0055   HltTree->Branch("MCphi", mcphi, "MCphi[NMCpart]/F");
0056   HltTree->Branch("MCPtHat", &pthatf, "MCPtHat/F");
0057   HltTree->Branch("MCWeight", &weightf, "MCWeight/F");
0058   HltTree->Branch("MCWeightSign", &weightsignf, "MCWeightSign/F");
0059   HltTree->Branch("MCmu3", &nmu3, "MCmu3/I");
0060   HltTree->Branch("MCel3", &nel3, "MCel3/I");
0061   HltTree->Branch("MCbb", &nbb, "MCbb/I");
0062   HltTree->Branch("MCab", &nab, "MCab/I");
0063   HltTree->Branch("MCWenu", &nwenu, "MCWenu/I");
0064   HltTree->Branch("MCWmunu", &nwmunu, "MCmunu/I");
0065   HltTree->Branch("MCZee", &nzee, "MCZee/I");
0066   HltTree->Branch("MCZmumu", &nzmumu, "MCZmumu/I");
0067   HltTree->Branch("MCptEleMax", &ptEleMax, "MCptEleMax/F");
0068   HltTree->Branch("MCptMuMax", &ptMuMax, "MCptMuMax/F");
0069   HltTree->Branch("NPUTrueBX0", &npubx0, "NPUTrueBX0/I");
0070   HltTree->Branch("NPUgenBX0", &npuvertbx0, "NPUgenBX0/I");
0071 }
0072
0073 /* **Analyze the event** */
0074 void HLTMCtruth::analyze(const edm::Handle<reco::CandidateView>& mctruth,
0075                          const double& pthat,
0076                          const double& weight,
0077                          const edm::Handle<std::vector<SimTrack> >& simTracks,
0078                          const edm::Handle<std::vector<SimVertex> >& simVertices,
0079                          const edm::Handle<std::vector<PileupSummaryInfo> >& PupInfo,
0080                          TTree* HltTree) {
0081   //std::cout << " Beginning HLTMCtruth " << std::endl;
0082
0083   if (_Monte) {
0084     int nmc = 0;
0085     int mu3 = 0;
0086     int el3 = 0;
0087     int mab = 0;
0088     int mbb = 0;
0089     int wel = 0;
0090     int wmu = 0;
0091     int zee = 0;
0092     int zmumu = 0;
0093
0094     ptEleMax = -999.0;
0095     ptMuMax = -999.0;
0096     pthatf = pthat;
0097     weightf = weight;
0098     weightsignf = (weight > 0) ? 1. : -1.;
0099     npubx0 = 0.0;
0100     npuvertbx0 = 0;
0101
0102     int npvtrue = 0;
0103     int npuvert = 0;
0104
0105     if (PupInfo.isValid()) {
0106       std::vector<PileupSummaryInfo>::const_iterator PVI;
0107       for (PVI = PupInfo->begin(); PVI != PupInfo->end(); ++PVI) {
0108         int BX = PVI->getBunchCrossing();
0109         npvtrue = PVI->getTrueNumInteractions();
0110         npuvert = PVI->getPU_NumInteractions();
0111
0112         if (BX == 0) {
0113           npubx0 += npvtrue;
0114           npuvertbx0 += npuvert;
0115         }
0116       }
0117     }
0118
0119     if ((simTracks.isValid()) && (simVertices.isValid())) {
0120       for (auto const& j : *simTracks) {
0121         int pdgid = j.type();
0122         if (abs(pdgid) != 13)
0123           continue;
0124         double pt = j.momentum().pt();
0125         if (pt < 5.0)
0126           continue;
0127         double eta = j.momentum().eta();
0128         if (std::abs(eta) > 2.5)
0129           continue;
0130         if (j.noVertex())
0131           continue;
0132         int vertIndex = j.vertIndex();
0133         double x = simVertices->at(vertIndex).position().x();
0134         double y = simVertices->at(vertIndex).position().y();
0135         double r = sqrt(x * x + y * y);
0136         if (r > 200.)
0137           continue;  // I think units are cm here
0138         double z = simVertices->at(vertIndex).position().z();
0139         if (std::abs(z) > 400.)
0140           continue;  // I think units are cm here
0141         mu3 += 1;
0142         break;
0143       }
0144     }
0145
0146     if (_Gen && mctruth.isValid()) {
0147       for (size_t i = 0; i < mctruth->size(); ++i) {
0148         const reco::Candidate& p = (*mctruth)[i];
0149
0150         mcpid[nmc] = p.pdgId();
0151         mcstatus[nmc] = p.status();
0152         mcpt[nmc] = p.pt();
0153         mceta[nmc] = p.eta();
0154         mcphi[nmc] = p.phi();
0155         mcvx[nmc] = p.vx();
0156         mcvy[nmc] = p.vy();
0157         mcvz[nmc] = p.vz();
0158
0159         if ((mcpid[nmc] == 24) || (mcpid[nmc] == -24)) {  // Checking W -> e/mu nu
0160           size_t idg = p.numberOfDaughters();
0161           for (size_t j = 0; j != idg; ++j) {
0162             const reco::Candidate& d = *p.daughter(j);
0163             if ((d.pdgId() == 11) || (d.pdgId() == -11)) {
0164               wel += 1;
0165             }
0166             if ((d.pdgId() == 13) || (d.pdgId() == -13)) {
0167               wmu += 1;
0168             }
0169             //      if ( (abs(d.pdgId())!=24) && ((mcpid[nmc])*(d.pdgId())>0) )
0170             //        {std::cout << "Wrong sign between mother-W and daughter !" << std::endl;}
0171           }
0172         }
0173         if (mcpid[nmc] == 23) {  // Checking Z -> 2 e/mu
0174           size_t idg = p.numberOfDaughters();
0175           for (size_t j = 0; j != idg; ++j) {
0176             const reco::Candidate& d = *p.daughter(j);
0177             if (d.pdgId() == 11) {
0178               zee += 1;
0179             }
0180             if (d.pdgId() == -11) {
0181               zee += 2;
0182             }
0183             if (d.pdgId() == 13) {
0184               zmumu += 1;
0185             }
0186             if (d.pdgId() == -13) {
0187               zmumu += 2;
0188             }
0189           }
0190         }
0191
0192         // Set-up flags, based on Pythia-generator information, for avoiding double-counting events when
0193         // using both pp->{e,mu}X AND QCD samples
0194         //  if (((mcpid[nmc]==13)||(mcpid[nmc]==-13))&&(mcpt[nmc]>2.5)) {mu3 += 1;} // Flag for muons with pT > 2.5 GeV/c
0195         if (((mcpid[nmc] == 11) || (mcpid[nmc] == -11)) && (mcpt[nmc] > 2.5)) {
0196           el3 += 1;
0197         }  // Flag for electrons with pT > 2.5 GeV/c
0198
0199         if (mcpid[nmc] == -5) {
0200           mab += 1;
0201         }  // Flag for bbar
0202         if (mcpid[nmc] == 5) {
0203           mbb += 1;
0204         }  // Flag for b
0205
0206         if ((mcpid[nmc] == 13) || (mcpid[nmc] == -13)) {
0207           if (p.pt() > ptMuMax) {
0208             ptMuMax = p.pt();
0209           }
0210         }  // Save max pt of generated Muons
0211         if ((mcpid[nmc] == 11) || (mcpid[nmc] == -11)) {
0212           if (p.pt() > ptEleMax)
0213             ptEleMax = p.pt();
0214         }  // Save max pt of generated Electrons
0215
0216         nmc++;
0217       }
0218     }
0219     //    else {std::cout << "%HLTMCtruth -- No MC truth information" << std::endl;}
0220
0221     nmcpart = nmc;
0222     nmu3 = mu3;
0223     nel3 = el3;
0224     nbb = mbb;
0225     nab = mab;
0226     nwenu = wel;
0227     nwmunu = wmu;
0228     if ((zee % 3) == 0) {
0229       nzee = zee / 3;
0230     }
0231     //     else {std::cout << "Z does not decay in e+ e- !" << std::endl;}
0232     if ((zmumu % 3) == 0) {
0233       nzmumu = zmumu / 3;
0234     }
0235     //     else {std::cout << "Z does not decay in mu+ mu- !" << std::endl;}
0236   }
0237 }