[ViewVC] Diff of: cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents):
Revision 1.1 by peiffer, Mon Apr 2 12:40:12 2012 UTC vs.
Revision 1.16 by peiffer, Mon Apr 30 09:01:16 2012 UTC

+//
+//
-<
+#include "UserCode/UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
->
+#include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
+  doMuons = iConfig.getParameter<bool>("doMuons");
+  doTaus = iConfig.getParameter<bool>("doTaus");
+  doJets = iConfig.getParameter<bool>("doJets");
-+
+  doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");
+  doPhotons = iConfig.getParameter<bool>("doPhotons");
+  doMET = iConfig.getParameter<bool>("doMET");
+  doGenInfo = iConfig.getParameter<bool>("doGenInfo");
-+
+  doAllGenParticles = iConfig.getParameter<bool>("doAllGenParticles");
-+
+  doLumiInfo = iConfig.getParameter<bool>("doLumiInfo");
+  doPV = iConfig.getParameter<bool>("doPV");
+  doTopJets = iConfig.getParameter<bool>("doTopJets");
-+
+  doTrigger = iConfig.getParameter<bool>("doTrigger");
+  // initialization of tree variables
+  tr->Branch("luminosityBlock",&luminosityBlock);
+  tr->Branch("isRealData",&isRealData);
+  tr->Branch("HBHENoiseFilterResult",&HBHENoiseFilterResult);
-<
+  tr->Branch("beamspot_x0",&beamspot_x0);
-<
+  tr->Branch("beamspot_y0",&beamspot_y0);
-<
+  tr->Branch("beamspot_z0",&beamspot_z0);
-<
->
+  if(doLumiInfo){
->
+    tr->Branch("intgRecLumi",&intgRecLumi);
->
+    tr->Branch("intgDelLumi",&intgDelLumi);
->
+  }
->
+  if(doPV){
->
+    tr->Branch("beamspot_x0",&beamspot_x0);
->
+    tr->Branch("beamspot_y0",&beamspot_y0);
->
+    tr->Branch("beamspot_z0",&beamspot_z0);
->
+  }
+  if(doElectrons){
+    electron_sources = iConfig.getParameter<std::vector<std::string> >("electron_sources");
+    for(size_t j=0; j< electron_sources.size(); ++j){
+      tr->Branch( topjet_sources[j].c_str(), "std::vector<TopJet>", &topjets[j]);
-+
+  if(doGenTopJets){
-+
+    gentopjet_sources = iConfig.getParameter<std::vector<std::string> >("gentopjet_sources");
-+
+    gentopjet_ptmin = iConfig.getParameter<double> ("gentopjet_ptmin");
-+
+    gentopjet_etamax = iConfig.getParameter<double> ("gentopjet_etamax");
-+
+    for(size_t j=0; j< gentopjet_sources.size(); ++j){
-+
+      tr->Branch( gentopjet_sources[j].c_str(), "std::vector<TopJet>", &gentopjets[j]);
-+
+    }
-+
+  }
+  if(doPhotons){
+    photon_sources = iConfig.getParameter<std::vector<std::string> >("photon_sources");
+    for(size_t j=0; j< photon_sources.size(); ++j){
+    tr->Branch("genInfo","GenInfo",&genInfo);
+    tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
-<
-<
+  trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
-<
+  //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
-<
+  tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
-<
+  tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
-<
+  tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
-<
+  tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
-<
->
+  if(doTrigger){
->
+    trigger_prefixes = iConfig.getParameter<std::vector<std::string> >("trigger_prefixes");
->
+    //tr->Branch("triggerResults","std::map<std::string, bool>",&triggerResults);
->
+    tr->Branch("triggerNames", "std::vector<std::string>", &triggerNames);
->
+    tr->Branch("triggerResults", "std::vector<bool>", &triggerResults);
->
+    //tr->Branch("L1_prescale", "std::vector<int>", &L1_prescale);
->
+    //tr->Branch("HLT_prescale", "std::vector<int>", &HLT_prescale);
->
+  }
+  newrun = true;
+   else HBHENoiseFilterResult = false;
-+
+   // ------------- primary vertices and beamspot  -------------
-+
-+
+   if(doPV){
-+
+     for(size_t j=0; j< pv_sources.size(); ++j){
-+
+       pvs[j].clear();
-+
-+
+       edm::Handle< std::vector<reco::Vertex> > pv_handle;
-+
+       iEvent.getByLabel(pv_sources[j], pv_handle);
-+
+       const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
-+
-+
+       for (unsigned int i = 0; i <  reco_pvs.size(); ++i) {
-+
+         reco::Vertex reco_pv = reco_pvs[i];
-+
-+
+         PrimaryVertex pv;
-+
+         pv.x =  reco_pv.x();
-+
+         pv.y =  reco_pv.y();
-+
+         pv.z =  reco_pv.z();
-+
+         pv.nTracks =  reco_pv.nTracks();
-+
+         //pv.isValid =  reco_pv.isValid();
-+
+         pv.chi2 =  reco_pv.chi2();
-+
+         pv.ndof =  reco_pv.ndof();
-+
-+
+         pvs[j].push_back(pv);
-+
+       }
-+
+     }
-+
-+
+     edm::Handle<reco::BeamSpot> beamSpot;
-+
+     iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
-+
+     const reco::BeamSpot & bsp = *beamSpot;
-+
-+
+     beamspot_x0 = bsp.x0();
-+
+     beamspot_y0 = bsp.y0();
-+
+     beamspot_z0 = bsp.z0();
-+
+   }
-+
-+
+ // ------------- generator info -------------
-+
-+
+   if(doGenInfo){
-+
+     genInfo.weights.clear();
-+
+     genInfo.binningValues.clear();
-+
+     genps.clear();
-+
-+
+     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
-+
+     iEvent.getByLabel("generator", genEventInfoProduct);
-+
+     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
-+
-+
+     genInfo.binningValues = genEventInfo.binningValues();
-+
+     genInfo.weights = genEventInfo.weights();
-+
+     genInfo.alphaQCD = genEventInfo.alphaQCD();
-+
+     genInfo.alphaQED = genEventInfo.alphaQED();
-+
+     genInfo.qScale = genEventInfo.qScale();
-+
-+
+     const gen::PdfInfo* pdf = genEventInfo.pdf();
-+
+     if(pdf){
-+
+       genInfo.pdf_id1=pdf->id.first;
-+
+       genInfo.pdf_id2=pdf->id.second;
-+
+       genInfo.pdf_x1=pdf->x.first;
-+
+       genInfo.pdf_x2=pdf->x.second;
-+
+       genInfo.pdf_scalePDF=pdf->scalePDF;
-+
+       genInfo.pdf_xPDF1=pdf->xPDF.first;
-+
+       genInfo.pdf_xPDF2=pdf->xPDF.second;
-+
+     }
-+
+     else{
-+
+       genInfo.pdf_id1=-999;
-+
+       genInfo.pdf_id2=-999;
-+
+       genInfo.pdf_x1=-999;
-+
+       genInfo.pdf_x2=-999;
-+
+       genInfo.pdf_scalePDF=-999;
-+
+       genInfo.pdf_xPDF1=-999;
-+
+       genInfo.pdf_xPDF2=-999;
-+
+     }
-+
-+
+     edm::Handle<std::vector<PileupSummaryInfo> > pus;
-+
+     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
-+
+     genInfo.pileup_NumInteractions_intime=0;
-+
+     genInfo.pileup_NumInteractions_ootbefore=0;
-+
+     genInfo.pileup_NumInteractions_ootafter=0;
-+
+     if(pus.isValid()){
-+
+       genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
-+
+       for(size_t i=0; i<pus->size(); ++i){
-+
+         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
-+
+            genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
-+
+         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
-+
+            genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
-+
+         }
-+
+         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
-+
+           genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
-+
+         }
-+
+       }
-+
+     }
-+
-+
+     edm::Handle<reco::GenParticleCollection> genPartColl;
-+
+     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
-+
+     int index=-1;
-+
+     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
-+
+       index++;
-+
-+
+       //write out only top quarks and status 3 particles (works fine only for MadGraph)
-+
+       if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
-+
+         GenParticle genp;
-+
+         genp.charge = iter->charge();;
-+
+         genp.pt = iter->p4().pt();
-+
+         genp.eta = iter->p4().eta();
-+
+         genp.phi = iter->p4().phi();
-+
+         genp.energy = iter->p4().E();
-+
+         genp.index =index;
-+
+         genp.status = iter->status();
-+
+         genp.pdgId = iter->pdgId();
-+
-+
+         genp.mother1=-1;
-+
+         genp.mother2=-1;
-+
+         genp.daughter1=-1;
-+
+         genp.daughter2=-1;
-+
-+
+         int nm=iter->numberOfMothers();
-+
+         int nd=iter->numberOfDaughters();
-+
-+
-+
+         if (nm>0) genp.mother1 = iter->motherRef(0).key();
-+
+         if (nm>1) genp.mother2 = iter->motherRef(1).key();
-+
+         if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
-+
+         if (nd>1) genp.daughter2 = iter->daughterRef(1).key();
-+
+         //std::cout << genp.index <<"  pdgId = " << genp.pdgId << "  mo1 = " << genp.mother1 << "  mo2 = " << genp.mother2 <<"  da1 = " << genp.daughter1 << "  da2 = " << genp.daughter2 <<std::endl;
-+
+         genps.push_back(genp);
-+
+       }
-+
+     }
-+
+   }
-+
+   // ------------- electrons -------------
+   if(doElectrons){
-+
-+
+     edm::Handle<reco::ConversionCollection> hConversions;
-+
+     iEvent.getByLabel("allConversions", hConversions);
-+
-+
+     edm::Handle<reco::BeamSpot> beamSpot;
-+
+     iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
-+
+     const reco::BeamSpot & bsp = *beamSpot;
-+
+     for(size_t j=0; j< electron_sources.size(); ++j){
+       eles[j].clear();
+       edm::Handle< std::vector<pat::Electron> > ele_handle;
+         ele.neutralHadronIso = pat_ele.neutralHadronIso();
+         ele.chargedHadronIso = pat_ele.chargedHadronIso();
+         ele.trackIso = pat_ele.trackIso();
-+
+         ele.photonIso = pat_ele.photonIso();
+         ele.puChargedHadronIso = pat_ele.puChargedHadronIso();
-+
+         ele.gsfTrack_trackerExpectedHitsInner_numberOfLostHits = pat_ele.gsfTrack()->trackerExpectedHitsInner().numberOfLostHits();
-+
+         ele.gsfTrack_px= pat_ele.gsfTrack()->px();
-+
+         ele.gsfTrack_py= pat_ele.gsfTrack()->py();
-+
+         ele.gsfTrack_pz= pat_ele.gsfTrack()->pz();
-+
+         ele.gsfTrack_vx= pat_ele.gsfTrack()->vx();
-+
+         ele.gsfTrack_vy= pat_ele.gsfTrack()->vy();
-+
+         ele.gsfTrack_vz= pat_ele.gsfTrack()->vz();
-+
+         ele.passconversionveto = !ConversionTools::hasMatchedConversion(pat_ele,hConversions,bsp.position());
-+
+         ele.dEtaIn=pat_ele.deltaEtaSuperClusterTrackAtVtx();
-+
+         ele.dPhiIn=pat_ele.deltaPhiSuperClusterTrackAtVtx();
-+
+         ele.sigmaIEtaIEta=pat_ele.sigmaIetaIeta();
-+
+         ele.HoverE=pat_ele.hadronicOverEm();
-+
+         ele.fbrem=pat_ele.fbrem();
-+
+         ele.EoverPIn=pat_ele.eSuperClusterOverP();
-+
+         ele.EcalEnergy=pat_ele.ecalEnergy();
-+
+         //ele.mvaTrigV0=pat_ele.electronID("mvaTrigV0");
-+
+         //ele.mvaNonTrigV0=pat_ele.electronID("mvaNonTrigV0");
+         eles[j].push_back(ele);
+         mu.neutralHadronIso = pat_mu.neutralHadronIso();
+         mu.chargedHadronIso = pat_mu.chargedHadronIso();
+         mu.trackIso = pat_mu.trackIso();
-+
+         mu.photonIso = pat_mu.photonIso();
+         mu.puChargedHadronIso = pat_mu.puChargedHadronIso();
+         mu.isGlobalMuon = pat_mu.isGlobalMuon();
+         mu.isStandAloneMuon = pat_mu.isStandAloneMuon();
+           mu.innerTrack_d0Error = innerTrack->d0Error();
+           mu.innerTrack_numberOfValidHits = innerTrack->numberOfValidHits();
+           mu.innerTrack_numberOfLostHits = innerTrack->numberOfLostHits();
-+
+           mu.innerTrack_trackerLayersWithMeasurement = innerTrack->hitPattern().trackerLayersWithMeasurement();
-+
+           mu.innerTrack_numberOfValidPixelHits = innerTrack->hitPattern().numberOfValidPixelHits();
+         else{
+           mu.innerTrack_chi2 = 0;
+           mu.innerTrack_d0Error = 0;
+           mu.innerTrack_numberOfValidHits = 0;
+           mu.innerTrack_numberOfLostHits = 0;
-+
+           mu.innerTrack_trackerLayersWithMeasurement = 0;
-+
+           mu.innerTrack_numberOfValidPixelHits = 0;
+         reco::TrackRef outerTrack = pat_mu.outerTrack();
+         if(!outerTrack.isNull()){
+         tau.eta =  pat_tau.eta();
+         tau.phi =  pat_tau.phi();
+         tau.energy =  pat_tau.energy();
-<
->
+         tau.decayModeFinding = pat_tau.tauID("decayModeFinding")>0.5;
->
+         tau.byVLooseCombinedIsolationDeltaBetaCorr  = pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5;
->
+         tau.byLooseCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5;
->
+         tau.byMediumCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5;
->
+         tau.byTightCombinedIsolationDeltaBetaCorr = pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5;
->
+         tau.againstElectronLoose  = pat_tau.tauID("againstElectronLoose")>0.5;
->
+         tau.againstElectronMedium = pat_tau.tauID("againstElectronMedium")>0.5;
->
+         tau.againstElectronTight = pat_tau.tauID("againstElectronTight")>0.5;
->
+         tau.againstElectronMVA  = pat_tau.tauID("againstElectronMVA")>0.5;
->
+         tau.againstMuonLoose = pat_tau.tauID("againstMuonLoose")>0.5;
->
+         tau.againstMuonMedium = pat_tau.tauID("againstMuonMedium")>0.5;
->
+         tau.againstMuonTight = pat_tau.tauID("againstMuonTight")>0.5;
->
->
+         reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
->
+         if(!leadPFCand.isNull()){
->
+           tau.leadPFCand_px = leadPFCand->px();
->
+           tau.leadPFCand_py = leadPFCand->py();
->
+           tau.leadPFCand_pz = leadPFCand->pz();
->
+         }
->
+         else{
->
+           tau.leadPFCand_px = 0;
->
+           tau.leadPFCand_py = 0;
->
+           tau.leadPFCand_pz = 0;
->
+         }
+         taus[j].push_back(tau);
+           jet.electronMultiplicity =pat_jet.electronMultiplicity();
+           jet.photonMultiplicity =pat_jet.photonMultiplicity();
-+
-+
+         jecUnc->setJetEta(pat_jet.eta());
-+
+         jecUnc->setJetPt(pat_jet.pt());
-+
+         jet.JEC_uncertainty = jecUnc->getUncertainty(true);
-+
+         jet.JEC_factor_raw = pat_jet.jecFactor("Uncorrected");
-+
+         jet.btag_simpleSecondaryVertexHighEff=pat_jet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
+         jet.btag_simpleSecondaryVertexHighPur=pat_jet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
+         jet.btag_combinedSecondaryVertex=pat_jet.bDiscriminator("combinedSecondaryVertexBJetTags");
+         jet.btag_jetBProbability=pat_jet.bDiscriminator("jetBProbabilityBJetTags");
+         jet.btag_jetProbability=pat_jet.bDiscriminator("jetProbabilityBJetTags");
-+
+         const reco::GenJet *genj = pat_jet.genJet();
-+
+         if(genj){
-+
+           jet.genjet_pt = genj->pt();
-+
+           jet.genjet_eta = genj->eta();
-+
+           jet.genjet_phi = genj->phi();
-+
+           jet.genjet_energy = genj->energy();
-+
+           if(doAllGenParticles){
-+
+             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
-+
+             for(unsigned int l = 0; l<jetgenps.size(); ++l){
-+
+               for(unsigned int k=0; k< genps.size(); ++k){
-+
+                 if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
-+
+                   jet.genparticles_indices.push_back(genps[k].index);
-+
+                 }
-+
+               }
-+
+             }
-+
+             if(jet.genparticles_indices.size()!= jetgenps.size())
-+
+               std::cout << "WARNING: Found only " << jet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
-+
+           }
-+
+         }
+         jets[j].push_back(jet);
+       iEvent.getByLabel(topjet_sources[j], pat_topjets);
+       for (unsigned int i = 0; i < pat_topjets->size(); i++) {
-+
+         const pat::Jet  pat_topjet =  * dynamic_cast<pat::Jet const *>(&pat_topjets->at(i));
+         if(pat_topjet.pt() < topjet_ptmin) continue;
+         if(fabs(pat_topjet.eta()) > topjet_etamax) continue;
+//       topjet.electronMultiplicity =pat_topjet.electronMultiplicity();
+//       topjet.photonMultiplicity =pat_topjet.photonMultiplicity();
-+
+         jecUnc->setJetEta(pat_topjet.eta());
-+
+         jecUnc->setJetPt(pat_topjet.pt());
-+
+         topjet.JEC_uncertainty = jecUnc->getUncertainty(true);
-+
+         topjet.JEC_factor_raw = pat_topjet.jecFactor("Uncorrected");
-+
+         topjet.btag_simpleSecondaryVertexHighEff=pat_topjet.bDiscriminator("simpleSecondaryVertexHighEffBJetTags");
+         topjet.btag_simpleSecondaryVertexHighPur=pat_topjet.bDiscriminator("simpleSecondaryVertexHighPurBJetTags");
+         topjet.btag_combinedSecondaryVertex=pat_topjet.bDiscriminator("combinedSecondaryVertexBJetTags");
+         topjet.btag_jetBProbability=pat_topjet.bDiscriminator("jetBProbabilityBJetTags");
+         topjet.btag_jetProbability=pat_topjet.bDiscriminator("jetProbabilityBJetTags");
-+
+         const reco::GenJet *genj = pat_topjet.genJet();
-+
+         if(genj){
-+
+           topjet.genjet_pt = genj->pt();
-+
+           topjet.genjet_eta = genj->eta();
-+
+           topjet.genjet_phi = genj->phi();
-+
+           topjet.genjet_energy = genj->energy();
-+
+           if(doAllGenParticles){
-+
+             std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
-+
+             for(unsigned int l = 0; l<jetgenps.size(); ++l){
-+
+               for(unsigned int k=0; k< genps.size(); ++k){
-+
+                 if(jetgenps[l]->pt() == genps[k].pt && jetgenps[l]->pdgId() == genps[k].pdgId){
-+
+                   topjet.genparticles_indices.push_back(genps[k].index);
-+
+                 }
-+
+               }
-+
+             }
-+
+             if(topjet.genparticles_indices.size()!= jetgenps.size())
-+
+               std::cout << "WARNING: Found only " << topjet.genparticles_indices.size() << " from " << jetgenps.size() << " gen particles of this topjet"<<std::endl;
-+
+           }
-+
+         }
-+
+         for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
+           Particle subjet_v4;
+           subjet_v4.pt = pat_topjet.daughter(k)->p4().pt();
-+
-+
+   // ------------- generator top jets -------------
-+
+   if(doGenTopJets){
-+
+     for(size_t j=0; j< gentopjet_sources.size(); ++j){
-+
-+
+       gentopjets[j].clear();
-+
-+
+       edm::Handle<reco::BasicJetCollection> reco_gentopjets;
-+
+       //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
-+
+       iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);
-+
-+
+       for (unsigned int i = 0; i < reco_gentopjets->size(); i++) {
-+
-+
+         const reco::BasicJet  reco_gentopjet =  reco_gentopjets->at(i);
-+
+         if(reco_gentopjet.pt() < gentopjet_ptmin) continue;
-+
+         if(fabs(reco_gentopjet.eta()) > gentopjet_etamax) continue;
-+
-+
+         TopJet gentopjet;
-+
+         gentopjet.charge = reco_gentopjet.charge();
-+
+         gentopjet.pt = reco_gentopjet.pt();
-+
+         gentopjet.eta = reco_gentopjet.eta();
-+
+         gentopjet.phi = reco_gentopjet.phi();
-+
+         gentopjet.energy = reco_gentopjet.energy();
-+
+         gentopjet.numberOfDaughters =reco_gentopjet.numberOfDaughters();
-+
-+
+         for (unsigned int k = 0; k < reco_gentopjet.numberOfDaughters(); k++) {
-+
+           Particle subjet_v4;
-+
+           subjet_v4.pt = reco_gentopjet.daughter(k)->p4().pt();
-+
+           subjet_v4.eta = reco_gentopjet.daughter(k)->p4().eta();
-+
+           subjet_v4.phi = reco_gentopjet.daughter(k)->p4().phi();
-+
+           subjet_v4.energy = reco_gentopjet.daughter(k)->p4().E();
-+
+           gentopjet.subjets.push_back(subjet_v4);
-+
+         }
-+
+         gentopjets[j].push_back(gentopjet);
-+
+       }
-+
+     }
-+
+   }
-+
+   // ------------- photons -------------
+   if(doPhotons){
+     for(size_t j=0; j< photon_sources.size(); ++j){
+   // ------------- trigger -------------
-<
-<
+   edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
-<
+   edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
-<
+   iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
-<
-<
+   const edm::Provenance *meta = dummy_TriggerEvent.provenance();
-<
+   std::string nameProcess = meta->processName();
-<
+   edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
-<
+   triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
-<
-<
+   edm::Handle<edm::TriggerResults> trigger;
-<
+   iEvent.getByLabel(triggerResultTag, trigger);
-<
+   const edm::TriggerResults& trig = *(trigger.product());
-<
-<
+   triggerResults.clear();
-<
+   triggerNames.clear();
-<
+   L1_prescale.clear();
-<
+   HLT_prescale.clear();
-<
-<
+   edm::Service<edm::service::TriggerNamesService> tns;
-<
+   std::vector<std::string> triggerNames_all;
-<
+   tns->getTrigPaths(trig,triggerNames_all);
-<
-<
+   if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
-<
+   for(unsigned int i=0; i<trig.size(); ++i){
-<
+     std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
-<
+     for(; it!=trigger_prefixes.end(); ++it){
-<
+       if(triggerNames_all[i].substr(0, it->size()) == *it)break;
-<
+     }
-<
+     if(it==trigger_prefixes.end()) continue;
-<
-<
+     //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
-<
+     triggerResults.push_back(trig.accept(i));
-<
+     if(newrun) triggerNames.push_back(triggerNames_all[i]);
-<
+     if(isRealData){
-<
+       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
-<
+       L1_prescale.push_back(pre.first);
-<
+       HLT_prescale.push_back(pre.second);
-<
+     }
-<
+   }
-<
+//    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
-<
+//      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
-<
+//    }
-<
+   newrun=false;
-<
-<
+   // ------------- generator info -------------
-<
-<
+   if(doGenInfo){
-<
+     genInfo.weights.clear();
-<
+     genInfo.binningValues.clear();
-<
+     genps.clear();
-<
-<
+     edm::Handle<GenEventInfoProduct> genEventInfoProduct;
-<
+     iEvent.getByLabel("generator", genEventInfoProduct);
-<
+     const GenEventInfoProduct& genEventInfo = *(genEventInfoProduct.product());
-<
-<
+     genInfo.binningValues = genEventInfo.binningValues();
-<
+     genInfo.weights = genEventInfo.weights();
-<
+     genInfo.alphaQCD = genEventInfo.alphaQCD();
-<
+     genInfo.alphaQED = genEventInfo.alphaQED();
-<
+     genInfo.qScale = genEventInfo.qScale();
->
+   if(doTrigger){
->
+     edm::InputTag triggerEvent = edm::InputTag("hltTriggerSummaryAOD");
->
+     edm::Handle< trigger::TriggerEvent > dummy_TriggerEvent;
->
+     iEvent.getByLabel( edm::InputTag(triggerEvent.label(), triggerEvent.instance()), dummy_TriggerEvent );
-<
+     const gen::PdfInfo* pdf = genEventInfo.pdf();
-<
+     if(pdf){
-<
+       genInfo.pdf_id1=pdf->id.first;
-<
+       genInfo.pdf_id2=pdf->id.second;
-<
+       genInfo.pdf_x1=pdf->x.first;
-<
+       genInfo.pdf_x2=pdf->x.second;
-<
+       genInfo.pdf_scalePDF=pdf->scalePDF;
-<
+       genInfo.pdf_xPDF1=pdf->xPDF.first;
-<
+       genInfo.pdf_xPDF2=pdf->xPDF.second;
-<
+     }
-<
+     else{
-<
+       genInfo.pdf_id1=-999;
-<
+       genInfo.pdf_id2=-999;
-<
+       genInfo.pdf_x1=-999;
-<
+       genInfo.pdf_x2=-999;
-<
+       genInfo.pdf_scalePDF=-999;
-<
+       genInfo.pdf_xPDF1=-999;
-<
+       genInfo.pdf_xPDF2=-999;
-<
+     }
-<
-<
+     edm::Handle<std::vector<PileupSummaryInfo> > pus;
-<
+     iEvent.getByLabel(edm::InputTag("addPileupInfo"), pus);
-<
+     genInfo.pileup_NumInteractions_intime=0;
-<
+     genInfo.pileup_NumInteractions_ootbefore=0;
-<
+     genInfo.pileup_NumInteractions_ootafter=0;
-<
+     if(pus.isValid()){
-<
+       genInfo.pileup_TrueNumInteractions = (float) pus->at(0).getTrueNumInteractions();
-<
+       for(size_t i=0; i<pus->size(); ++i){
-<
+         if(pus->at(i).getBunchCrossing() == 0) // intime pileup
-<
+            genInfo.pileup_NumInteractions_intime += pus->at(i).getPU_NumInteractions();
-<
+         else if(pus->at(i).getBunchCrossing() == -1){ // oot pileup before
-<
+            genInfo.pileup_NumInteractions_ootbefore += pus->at(i).getPU_NumInteractions();
-<
+         }
-<
+         else if(pus->at(i).getBunchCrossing() == +1){ // oot pileup before
-<
+           genInfo.pileup_NumInteractions_ootafter += pus->at(i).getPU_NumInteractions();
-<
+         }
-<
+       }
-<
+     }
-<
-<
+     edm::Handle<reco::GenParticleCollection> genPartColl;
-<
+     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
-<
+     int index=-1;
-<
+     for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
-<
+       index++;
-<
-<
+       //write out only top quarks and status 3 particles (works fine only for MadGraph)
-<
+       if(abs(iter->pdgId())==6 || iter->status()==3){
-<
+         GenParticle genp;
-<
+         genp.charge = iter->charge();;
-<
+         genp.pt = iter->p4().pt();
-<
+         genp.eta = iter->p4().eta();
-<
+         genp.phi = iter->p4().phi();
-<
+         genp.energy = iter->p4().E();
-<
+         genp.index =index;
-<
+         genp.status = iter->status();
-<
+         genp.pdgId = iter->pdgId();
-<
-<
+         genp.mother1=-1;
-<
+         genp.mother2=-1;
-<
+         genp.daughter1=-1;
-<
+         genp.daughter2=-1;
-<
-<
+         int nm=iter->numberOfMothers();
-<
+         int nd=iter->numberOfDaughters();
-<
-<
+         if (nm>0) genp.mother1 = iter->motherRef(0).key();
-<
+         if (nm>1) genp.mother2 = iter->motherRef(nm-1).key();
-<
+         if (nd>0) genp.daughter1 = iter->daughterRef(0).key();
-<
+         if (nd>1) genp.daughter2 = iter->daughterRef(nd-1).key();
-<
-<
+         genps.push_back(genp);
-<
+       }
-<
+     }
-<
->
+     const edm::Provenance *meta = dummy_TriggerEvent.provenance();
->
+     std::string nameProcess = meta->processName();
->
+     edm::InputTag triggerResultTag = edm::InputTag("TriggerResults");
->
+     triggerResultTag = edm::InputTag( triggerResultTag.label(), triggerResultTag.instance(), nameProcess );
->
->
+     edm::Handle<edm::TriggerResults> trigger;
->
+     iEvent.getByLabel(triggerResultTag, trigger);
->
+     const edm::TriggerResults& trig = *(trigger.product());
->
->
+     triggerResults.clear();
->
+     triggerNames.clear();
->
+//      L1_prescale.clear();
->
+//      HLT_prescale.clear();
->
->
+     edm::Service<edm::service::TriggerNamesService> tns;
->
+     std::vector<std::string> triggerNames_all;
->
+     tns->getTrigPaths(trig,triggerNames_all);
->
->
+     if (trig.size()!=triggerNames_all.size()) std::cout <<"ERROR: length of names and paths not the same: "<<triggerNames_all.size()<<","<<trig.size()<< std::endl;
->
+     for(unsigned int i=0; i<trig.size(); ++i){
->
+       std::vector<std::string>::const_iterator it = trigger_prefixes.begin();
->
+       for(; it!=trigger_prefixes.end(); ++it){
->
+         if(triggerNames_all[i].substr(0, it->size()) == *it)break;
->
+       }
->
+       if(it==trigger_prefixes.end()) continue;
->
->
+       //triggerResults.insert(std::pair<std::string, bool>(triggerNames[i],trig.accept(i)));
->
+       triggerResults.push_back(trig.accept(i));
->
+       if(newrun) triggerNames.push_back(triggerNames_all[i]);
->
+//        if(isRealData){
->
+//       std::pair<int, int> pre=hlt_cfg.prescaleValues(iEvent, iSetup, triggerNames_all[i]);
->
+//       L1_prescale.push_back(pre.first);
->
+//       HLT_prescale.push_back(pre.second);
->
+//       //std::cout <<  triggerNames_all[i] << " " << pre.first << " " <<pre.second << "   " << hlt_cfg.prescaleValue(iEvent, iSetup, triggerNames_all[i]) << std::endl;
->
+//        }
->
+     }
->
+     //    for(std::map<std::string, bool>::const_iterator iter = triggerResults.begin(); iter!=triggerResults.end(); iter++){
->
+     //      std::cout << (*iter).first << "   " << (*iter).second << std::endl;
->
+     //    }
->
+     newrun=false;
-–
+   // ------------- primary vertices and beamspot  -------------
-–
-–
+   if(doPV){
-–
+     for(size_t j=0; j< pv_sources.size(); ++j){
-–
+       pvs[j].clear();
-–
-–
+       edm::Handle< std::vector<reco::Vertex> > pv_handle;
-–
+       iEvent.getByLabel(pv_sources[j], pv_handle);
-–
+       const std::vector<reco::Vertex>& reco_pvs = *(pv_handle.product());
-–
-–
+       for (unsigned int i = 0; i <  reco_pvs.size(); ++i) {
-–
+         reco::Vertex reco_pv = reco_pvs[i];
-–
-–
+         PrimaryVertex pv;
-–
+         pv.x =  reco_pv.x();
-–
+         pv.y =  reco_pv.y();
-–
+         pv.z =  reco_pv.z();
-–
+         pv.nTracks =  reco_pv.nTracks();
-–
+         //pv.isValid =  reco_pv.isValid();
-–
+         pv.chi2 =  reco_pv.chi2();
-–
+         pv.ndof =  reco_pv.ndof();
-–
-–
+         pvs[j].push_back(pv);
-–
+       }
-–
+     }
-–
+   }
-–
-–
+   edm::Handle<reco::BeamSpot> beamSpot;
-–
+   iEvent.getByLabel(edm::InputTag("offlineBeamSpot"), beamSpot);
-–
+   const reco::BeamSpot & bsp = *beamSpot;
-–
-–
+   beamspot_x0 = bsp.x0();
-–
+   beamspot_y0 = bsp.y0();
-–
+   beamspot_z0 = bsp.z0();
+   tr->Fill();
-<
->
+   if(doLumiInfo)
->
+     previouslumiblockwasfilled=true;
+void
+NtupleWriter::beginJob()
-+
+  if(doLumiInfo){
-+
+    totalRecLumi=0;
-+
+    totalDelLumi=0;
-+
+    previouslumiblockwasfilled=false;
-+
+  }
+// ------------ method called once each job just after ending the event loop  ------------
+void
+NtupleWriter::beginRun(edm::Run const& iRun, edm::EventSetup const&  iSetup)
-<
+  bool setup_changed = false;
-<
+  hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
-<
+  newrun=true;
->
+  if(doTrigger){
->
+    bool setup_changed = false;
->
+    hlt_cfg.init(iRun, iSetup, "HLT", setup_changed);
->
+    newrun=true;
->
+  }
->
->
+  if(doJets || doTopJets){
->
+    edm::ESHandle<JetCorrectorParametersCollection> JetCorParColl;
->
+    iSetup.get<JetCorrectionsRecord>().get("AK5PF",JetCorParColl);
->
+    JetCorrectorParameters const & JetCorPar = (*JetCorParColl)["Uncertainty"];
->
+    jecUnc = new JetCorrectionUncertainty(JetCorPar);
->
+  }
+// ------------ method called when ending the processing of a run  ------------
+void
+NtupleWriter::endRun(edm::Run const&, edm::EventSetup const&)
-+
+  if(doLumiInfo)
-+
+    std::cout << "total integ. luminosity: " << totalDelLumi <<"(del) " << totalRecLumi << "(rec)" << std::endl;
+// ------------ method called when starting to processes a luminosity block  ------------
+void
-<
+NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&)
->
+NtupleWriter::beginLuminosityBlock(edm::LuminosityBlock const& lumi, edm::EventSetup const&)
-+
+  if(doLumiInfo){
-+
+    edm::Handle<LumiSummary> l;
-+
+    lumi.getByLabel("lumiProducer", l);
-+
-+
+    //add lumi of lumi blocks without any event to next lumiblock
-+
+    if(previouslumiblockwasfilled){
-+
+      intgRecLumi=0;
-+
+      intgDelLumi=0;
-+
+    }
-+
+    previouslumiblockwasfilled=false;
-+
-+
+    if (l.isValid()){;
-+
+      intgRecLumi+=l->intgRecLumi()*6.37;
-+
+      intgDelLumi+=l->intgDelLumi()*6.37;
-+
+      totalRecLumi+=l->intgRecLumi()*6.37;
-+
+      totalDelLumi+=l->intgDelLumi()*6.37;
-+
+    }
-+
+    //std::cout << "this lb: " <<l->intgRecLumi()*6.37 <<"   " << l->intgDelLumi()*6.37<<std::endl;
-+
+    //std::cout << "summed:  "<< intgRecLumi << "   " << intgDelLumi << std::endl;
-+
+  }
+// ------------ method called when ending the processing of a luminosity block  ------------

Diff Legend

-–
+Removed lines
-+
+Added lines
-<
+Changed lines
->
+Changed lines

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents): Revision 1.1 by peiffer, Mon Apr 2 12:40:12 2012 UTC vs. Revision 1.16 by peiffer, Mon Apr 30 09:01:16 2012 UTC

Diff Legend

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents):
Revision 1.1 by peiffer, Mon Apr 2 12:40:12 2012 UTC vs.
Revision 1.16 by peiffer, Mon Apr 30 09:01:16 2012 UTC