ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc
(Generate patch)

Comparing UserCode/UHHAnalysis/NtupleWriter/src/NtupleWriter.cc (file contents):
Revision 1.21 by peiffer, Tue Jun 26 08:13:28 2012 UTC vs.
Revision 1.26 by peiffer, Mon Jun 10 13:33:05 2013 UTC

# Line 17 | Line 17
17   //
18   //
19  
20 < #include "UHHAnalysis//NtupleWriter/interface/NtupleWriter.h"
21 <
22 <
20 > #include "UHHAnalysis/NtupleWriter/interface/NtupleWriter.h"
21 > #include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputer.h"
22 > #include "RecoBTau/JetTagComputer/interface/GenericMVAJetTagComputerWrapper.h"
23 > #include "RecoBTau/JetTagComputer/interface/JetTagComputer.h"
24 > #include "RecoBTau/JetTagComputer/interface/JetTagComputerRecord.h"
25 > #include "RecoBTag/SecondaryVertex/interface/CombinedSVComputer.h"
26  
27   //
28   // constants, enums and typedefs
# Line 50 | Line 53 | NtupleWriter::NtupleWriter(const edm::Pa
53    doMuons = iConfig.getParameter<bool>("doMuons");
54    doTaus = iConfig.getParameter<bool>("doTaus");
55    doJets = iConfig.getParameter<bool>("doJets");
56 +  doGenJets = iConfig.getParameter<bool>("doGenJets");
57    doJECUncertainty = iConfig.getParameter<bool>("doJECUncertainty");
58    doGenTopJets = iConfig.getParameter<bool>("doGenTopJets");  
59    doPhotons = iConfig.getParameter<bool>("doPhotons");
# Line 60 | Line 64 | NtupleWriter::NtupleWriter(const edm::Pa
64    doPV = iConfig.getParameter<bool>("doPV");
65    doTopJets = iConfig.getParameter<bool>("doTopJets");
66    doTrigger = iConfig.getParameter<bool>("doTrigger");
67 <
67 >  SVComputer_  = iConfig.getUntrackedParameter<edm::InputTag>("svComputer",edm::InputTag("combinedSecondaryVertex"));
68 >  doTagInfos = iConfig.getUntrackedParameter<bool>("doTagInfos",false);
69    // initialization of tree variables
70  
71    tr->Branch("run",&run);
# Line 108 | Line 113 | NtupleWriter::NtupleWriter(const edm::Pa
113        tr->Branch( jet_sources[j].c_str(), "std::vector<Jet>", &jets[j]);
114      }
115    }
116 +  if(doGenJets){
117 +    genjet_sources = iConfig.getParameter<std::vector<std::string> >("genjet_sources");
118 +    genjet_ptmin = iConfig.getParameter<double> ("genjet_ptmin");
119 +    genjet_etamax = iConfig.getParameter<double> ("genjet_etamax");
120 +    for(size_t j=0; j< genjet_sources.size(); ++j){  
121 +      tr->Branch( genjet_sources[j].c_str(), "std::vector<Particle>", &genjets[j]);
122 +    }
123 +  }
124    if(doTopJets){
125      topjet_sources = iConfig.getParameter<std::vector<std::string> >("topjet_sources");
126      topjet_ptmin = iConfig.getParameter<double> ("topjet_ptmin");
# Line 143 | Line 156 | NtupleWriter::NtupleWriter(const edm::Pa
156      }
157    }
158    if(doGenInfo){
159 +    genparticle_source= iConfig.getParameter<edm::InputTag>("genparticle_source");
160      tr->Branch("genInfo","GenInfo",&genInfo);
161      tr->Branch("GenParticles","std::vector<GenParticle>", &genps);
162    }
# Line 292 | Line 306 | NtupleWriter::analyze(const edm::Event&
306       }
307  
308       edm::Handle<reco::GenParticleCollection> genPartColl;
309 <     iEvent.getByLabel(edm::InputTag("genParticles"), genPartColl);
309 >     iEvent.getByLabel(genparticle_source, genPartColl);
310       int index=-1;
311       for(reco::GenParticleCollection::const_iterator iter = genPartColl->begin(); iter != genPartColl->end(); ++ iter){
312         index++;
313        
314 <       //write out only top quarks and status 3 particles (works fine only for MadGraph)
315 <       if(abs(iter->pdgId())==6 || iter->status()==3 || doAllGenParticles){
314 >       //write out only top quarks,final state leptons and status 3 particles (works fine only for MadGraph)
315 >       bool islepton = iter->status()==1 && abs(iter->pdgId())>=11 && abs(iter->pdgId())<=16 ;
316 >       if(abs(iter->pdgId())==6 || iter->status()==3 || islepton ||  doAllGenParticles){
317           GenParticle genp;
318           genp.set_charge(iter->charge());
319           genp.set_pt(iter->p4().pt());
# Line 383 | Line 398 | NtupleWriter::analyze(const edm::Event&
398           ele.set_EcalEnergy(pat_ele.ecalEnergy());
399           ele.set_mvaTrigV0(pat_ele.electronID("mvaTrigV0"));
400           ele.set_mvaNonTrigV0(pat_ele.electronID("mvaNonTrigV0"));
401 +         float AEff03 = 0.00;
402 +         if(isRealData){
403 +           AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAData2011);
404 +         }else{
405 +           AEff03 = ElectronEffectiveArea::GetElectronEffectiveArea(ElectronEffectiveArea::kEleGammaAndNeutralHadronIso03, pat_ele.superCluster()->eta(), ElectronEffectiveArea::kEleEAFall11MC);
406 +         }
407 +         ele.set_AEff(AEff03);
408  
409           eles[j].push_back(ele);
410         }
# Line 506 | Line 528 | NtupleWriter::analyze(const edm::Event&
528           tau.set_phi( pat_tau.phi());
529           tau.set_energy( pat_tau.energy());
530           tau.set_decayModeFinding ( pat_tau.tauID("decayModeFinding")>0.5);
531 <         tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
531 >         //tau.set_byVLooseCombinedIsolationDeltaBetaCorr  ( pat_tau.tauID("byVLooseCombinedIsolationDeltaBetaCorr")>0.5);
532           tau.set_byLooseCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr")>0.5);
533           tau.set_byMediumCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr")>0.5);
534           tau.set_byTightCombinedIsolationDeltaBetaCorr ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr")>0.5);
535 <         tau.set_againstElectronLoose  ( pat_tau.tauID("againstElectronLoose")>0.5);
536 <         tau.set_againstElectronMedium ( pat_tau.tauID("againstElectronMedium")>0.5);
537 <         tau.set_againstElectronTight ( pat_tau.tauID("againstElectronTight")>0.5);
538 <         tau.set_againstElectronMVA  ( pat_tau.tauID("againstElectronMVA")>0.5);
539 <         tau.set_againstMuonLoose ( pat_tau.tauID("againstMuonLoose")>0.5);
540 <         tau.set_againstMuonMedium ( pat_tau.tauID("againstMuonMedium")>0.5);
541 <         tau.set_againstMuonTight ( pat_tau.tauID("againstMuonTight")>0.5);
535 >         tau.set_byLooseIsolationMVA( pat_tau.tauID("byLooseIsolationMVA")>0.5);
536 >         tau.set_byMediumIsolationMVA( pat_tau.tauID("byMediumIsolationMVA")>0.5);
537 >         tau.set_byTightIsolationMVA( pat_tau.tauID("byTightIsolationMVA")>0.5);
538 >         tau.set_byLooseIsolationMVA2( pat_tau.tauID("byLooseIsolationMVA2")>0.5);
539 >         tau.set_byMediumIsolationMVA2( pat_tau.tauID("byMediumIsolationMVA2")>0.5);
540 >         tau.set_byTightIsolationMVA2( pat_tau.tauID("byTightIsolationMVA2")>0.5);
541 >         tau.set_byLooseCombinedIsolationDeltaBetaCorr3Hits(  pat_tau.tauID("byLooseCombinedIsolationDeltaBetaCorr3Hits")>0.5);
542 >         tau.set_byMediumCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byMediumCombinedIsolationDeltaBetaCorr3Hits")>0.5);
543 >         tau.set_byTightCombinedIsolationDeltaBetaCorr3Hits ( pat_tau.tauID("byTightCombinedIsolationDeltaBetaCorr3Hits")>0.5);
544 >         tau.set_againstElectronLooseMVA3  ( pat_tau.tauID("againstElectronLooseMVA3")>0.5);
545 >         tau.set_againstElectronMediumMVA3 ( pat_tau.tauID("againstElectronMediumMVA3")>0.5);
546 >         tau.set_againstElectronTightMVA3 ( pat_tau.tauID("againstElectronTightMVA3")>0.5);
547 >         tau.set_againstElectronVTightMVA3 ( pat_tau.tauID("againstElectronVTightMVA3")>0.5);
548 >         tau.set_againstMuonLoose2 ( pat_tau.tauID("againstMuonLoose2")>0.5);
549 >         tau.set_againstMuonMedium2 ( pat_tau.tauID("againstMuonMedium2")>0.5);
550 >         tau.set_againstMuonTight2 ( pat_tau.tauID("againstMuonTight2")>0.5);
551 >         tau.set_byIsolationMVAraw(  pat_tau.tauID("byIsolationMVAraw"));
552 >         tau.set_byIsolationMVA2raw(  pat_tau.tauID("byIsolationMVA2raw"));
553 >         tau.set_decayMode( pat_tau.decayMode() );
554 >         tau.set_byCombinedIsolationDeltaBetaCorrRaw( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw"));
555 >         tau.set_byCombinedIsolationDeltaBetaCorrRaw3Hits( pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits"));
556  
557 + //       std::cout << pat_tau.tauID("byCombinedIsolationDeltaBetaCorrRaw3Hits") << std::endl;
558 +        
559   //       reco::PFCandidateRef leadPFCand = pat_tau.leadPFCand();
560   //       if(!leadPFCand.isNull()){
561   //         tau.set_leadPFCand_px ( leadPFCand->px());
# Line 534 | Line 572 | NtupleWriter::analyze(const edm::Event&
572       }
573     }
574  
575 +   //-------------- gen jets -------------
576 +
577 +   if(doGenJets){
578 +     for(size_t j=0; j< genjet_sources.size(); ++j){
579 +      
580 +       genjets[j].clear();
581 +
582 +       edm::Handle< std::vector<reco::GenJet> > genjet_handle;
583 +       iEvent.getByLabel(genjet_sources[j], genjet_handle);
584 +       const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
585 +  
586 +       for (unsigned int i = 0; i < gen_jets.size(); ++i) {
587 +         const reco::GenJet* gen_jet = &gen_jets[i];
588 +         if(gen_jet->pt() < genjet_ptmin) continue;
589 +         if(fabs(gen_jet->eta()) > genjet_etamax) continue;
590 +
591 +         Particle jet;
592 +         jet.set_charge(gen_jet->charge());
593 +         jet.set_pt(gen_jet->pt());
594 +         jet.set_eta(gen_jet->eta());
595 +         jet.set_phi(gen_jet->phi());
596 +         jet.set_energy(gen_jet->energy());
597 +
598 +         genjets[j].push_back(jet);
599 +
600 +       }
601 +     }
602 +   }
603 +
604     // ------------- jets -------------
605     if(doJets){
606       for(size_t j=0; j< jet_sources.size(); ++j){
# Line 592 | Line 659 | NtupleWriter::analyze(const edm::Event&
659           jet.set_btag_jetBProbability(pat_jet.bDiscriminator("jetBProbabilityBJetTags"));
660           jet.set_btag_jetProbability(pat_jet.bDiscriminator("jetProbabilityBJetTags"));
661  
662 +        
663           const reco::GenJet *genj = pat_jet.genJet();
664           if(genj){
665 <           jet.set_genjet_pt(genj->pt());
666 <           jet.set_genjet_eta(genj->eta());
667 <           jet.set_genjet_phi(genj->phi());
668 <           jet.set_genjet_energy(genj->energy());
665 >
666 >           for(unsigned int k=0; k<genjets->size(); ++k){
667 >             if(genj->pt()==genjets->at(k).pt() && genj->eta()==genjets->at(k).eta())
668 >               jet.set_genjet_index(k);
669 >           }
670 > //         if( jet.genjet_index()<0){
671 > //           std::cout<< "genjet not found for " << genj->pt() << "  " << genj->eta() << std::endl;
672 > //         }
673 >
674             if(doAllGenParticles){
675               std::vector<const reco::GenParticle * > jetgenps = genj->getGenConstituents();
676               for(unsigned int l = 0; l<jetgenps.size(); ++l){
# Line 611 | Line 684 | NtupleWriter::analyze(const edm::Event&
684               if(jet.genparticles_indices().size()!= jetgenps.size())
685                 std::cout << "WARNING: Found only " << jet.genparticles_indices().size() << " from " << jetgenps.size() << " gen particles of this jet"<<std::endl;
686             }
687 +          
688           }
689 +        
690           jets[j].push_back(jet);
691         }
692       }
693     }
694  
695 +
696 +
697     // ------------- top jets -------------
698     if(doTopJets){
699       for(size_t j=0; j< topjet_sources.size(); ++j){
# Line 694 | Line 771 | NtupleWriter::analyze(const edm::Event&
771  
772           for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
773             Particle subjet_v4;
774 <           subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
775 <           subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
776 <           subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
777 <           subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
778 <           topjet.add_subjet(subjet_v4);
774 >
775 >           reco::Candidate const * subjetd =  pat_topjet.daughter(k);
776 >           pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
777 >           if(patsubjetd)
778 >           {
779 >              subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
780 >              subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
781 >              subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
782 >              subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
783 >              topjet.add_subjet(subjet_v4);
784 >              //btag info
785 >              topjet.add_subFlavour(patsubjetd->partonFlavour());
786 >              topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
787 >              if (doTagInfos)
788 >                {
789 >                  //ip tag info
790 >                  reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
791 >                  topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
792 >                  topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
793 >                  topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
794 >                  topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
795 >                  topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
796 >                  topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
797 >                  topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
798 >                  topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
799 >                  topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
800 >                  topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
801 >                  topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
802 >                  topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));    
803 >                  topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
804 >                  topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
805 >                  topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
806 >                  topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
807 >                  topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
808 >                  topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
809 >                  topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
810 >                  topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
811 >                  topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
812 >                  //sv tag info
813 >                  reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
814 >                  topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
815 >                  topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
816 >                  topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
817 >                  topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
818 >                  topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
819 >                  topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
820 >                  topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
821 >                  std::vector<TLorentzVector> vp4; vp4.clear();
822 >                  std::vector<float> vchi2; vchi2.clear();
823 >                  std::vector<float> vndof; vndof.clear();
824 >                  std::vector<float> vchi2ndof; vchi2ndof.clear();
825 >                  std::vector<float> vtrsize; vtrsize.clear();
826 >                  for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
827 >                    {
828 >                      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
829 >                      vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
830 >                      vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());  
831 >                      vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());  
832 >                      vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());  
833 >                      vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());  
834 >                    }
835 >                  topjet.add_subSecondaryVertex(vp4);
836 >                  topjet.add_subVertexChi2(vchi2);
837 >                  topjet.add_subVertexNdof(vndof);
838 >                  topjet.add_subVertexNormalizedChi2(vchi2ndof);
839 >                  topjet.add_subVertexTracksSize(vtrsize);
840 >                  //try computer
841 >                  edm::ESHandle<JetTagComputer> computerHandle;
842 >                  iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
843 >                  const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
844 >                  if(computer)
845 >                    {
846 >                      computer->passEventSetup(iSetup);
847 >                      std::vector<const reco::BaseTagInfo*>  baseTagInfos;
848 >                      baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
849 >                      baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );      
850 >                      JetTagComputer::TagInfoHelper helper(baseTagInfos);
851 >                      reco::TaggingVariableList vars = computer->taggingVariables(helper);
852 >                      topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
853 >                      topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
854 >                      topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
855 >                      topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
856 >                    }
857 >                }
858 >           }
859 >           else
860 >             {
861 >               //filling only standard information in case the subjet has not been pat-tified during the pattuples production
862 >               subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
863 >               subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
864 >               subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
865 >               subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
866 >               topjet.add_subjet(subjet_v4);
867 >             }
868 >          
869 >          
870           }
871           topjets[j].push_back(topjet);
872         }
873       }
874     }
875 <
876 <
875 >  
876 >  
877     // ------------- generator top jets -------------
878     if(doGenTopJets){
879       for(size_t j=0; j< gentopjet_sources.size(); ++j){
880        
881         gentopjets[j].clear();
882 <
882 >      
883         edm::Handle<reco::BasicJetCollection> reco_gentopjets;
884         //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
885         iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines