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.22 by peiffer, Wed Jul 25 09:56:57 2012 UTC vs.
Revision 1.24 by peiffer, Wed Jun 5 14:13:58 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 628 | Line 641 | NtupleWriter::analyze(const edm::Event&
641       }
642     }
643  
644 +   //-------------- gen jets -------------
645 +
646 +   if(doGenJets){
647 +     for(size_t j=0; j< genjet_sources.size(); ++j){
648 +      
649 +       genjets[j].clear();
650 +
651 +       edm::Handle< std::vector<reco::GenJet> > genjet_handle;
652 +       iEvent.getByLabel(genjet_sources[j], genjet_handle);
653 +       const std::vector<reco::GenJet>& gen_jets = *(genjet_handle.product());
654 +  
655 +       for (unsigned int i = 0; i < gen_jets.size(); ++i) {
656 +         pat::Jet gen_jet = gen_jets[i];
657 +         if(gen_jet.pt() < genjet_ptmin) continue;
658 +         if(fabs(gen_jet.eta()) > genjet_etamax) continue;
659 +
660 +         Particle jet;
661 +         jet.set_charge(gen_jet.charge());
662 +         jet.set_pt(gen_jet.pt());
663 +         jet.set_eta(gen_jet.eta());
664 +         jet.set_phi(gen_jet.phi());
665 +         jet.set_energy(gen_jet.energy());
666 +
667 +         genjets[j].push_back(jet);
668 +
669 +       }
670 +     }
671 +   }
672 +
673     // ------------- top jets -------------
674     if(doTopJets){
675       for(size_t j=0; j< topjet_sources.size(); ++j){
# Line 705 | Line 747 | NtupleWriter::analyze(const edm::Event&
747  
748           for (unsigned int k = 0; k < pat_topjet.numberOfDaughters(); k++) {
749             Particle subjet_v4;
750 <           subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
751 <           subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
752 <           subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
753 <           subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
754 <           topjet.add_subjet(subjet_v4);
750 >
751 >           reco::Candidate const * subjetd =  pat_topjet.daughter(k);
752 >           pat::Jet const * patsubjetd = dynamic_cast<pat::Jet const *>(subjetd);
753 >           if(patsubjetd)
754 >           {
755 >              subjet_v4.set_pt(patsubjetd->correctedP4(0).pt());
756 >              subjet_v4.set_eta(patsubjetd->correctedP4(0).eta());
757 >              subjet_v4.set_phi(patsubjetd->correctedP4(0).phi());
758 >              subjet_v4.set_energy(patsubjetd->correctedP4(0).E());
759 >              topjet.add_subjet(subjet_v4);
760 >              //btag info
761 >              topjet.add_subFlavour(patsubjetd->partonFlavour());
762 >              topjet.add_subCSV(patsubjetd->bDiscriminator("combinedSecondaryVertexBJetTags"));
763 >              if (doTagInfos)
764 >                {
765 >                  //ip tag info
766 >                  reco::TaggingVariableList tvlIP=patsubjetd->tagInfoTrackIP("impactParameter")->taggingVariables();
767 >                  topjet.add_subTrackMomentum(tvlIP.getList(reco::btau::trackMomentum,false));
768 >                  topjet.add_subTrackEta(tvlIP.getList(reco::btau::trackEta,false));
769 >                  topjet.add_subTrackEtaRel(tvlIP.getList(reco::btau::trackEtaRel,false));
770 >                  topjet.add_subTrackDeltaR(tvlIP.getList(reco::btau::trackDeltaR,false));
771 >                  topjet.add_subTrackSip3dVal(tvlIP.getList(reco::btau::trackSip3dVal,false));
772 >                  topjet.add_subTrackSip3dSig(tvlIP.getList(reco::btau::trackSip3dSig,false));
773 >                  topjet.add_subTrackSip2dVal(tvlIP.getList(reco::btau::trackSip2dVal,false));
774 >                  topjet.add_subTrackSip2dSig(tvlIP.getList(reco::btau::trackSip2dSig,false));
775 >                  topjet.add_subTrackDecayLenVal(tvlIP.getList(reco::btau::trackDecayLenVal,false));
776 >                  topjet.add_subTrackChi2(tvlIP.getList(reco::btau::trackChi2,false));
777 >                  topjet.add_subTrackNTotalHits(tvlIP.getList(reco::btau::trackNTotalHits,false));
778 >                  topjet.add_subTrackNPixelHits(tvlIP.getList(reco::btau::trackNPixelHits,false));    
779 >                  topjet.add_subTrackPtRel(tvlIP.getList(reco::btau::trackPtRel,false));
780 >                  topjet.add_subTrackPPar(tvlIP.getList(reco::btau::trackPPar,false));
781 >                  topjet.add_subTrackPtRatio(tvlIP.getList(reco::btau::trackPtRatio,false));
782 >                  topjet.add_subTrackPParRatio(tvlIP.getList(reco::btau::trackPParRatio,false));
783 >                  topjet.add_subTrackJetDistVal(tvlIP.getList(reco::btau::trackJetDistVal,false));
784 >                  topjet.add_subTrackJetDistSig(tvlIP.getList(reco::btau::trackJetDistSig,false));
785 >                  topjet.add_subTrackGhostTrackDistVal(tvlIP.getList(reco::btau::trackGhostTrackDistVal,false));
786 >                  topjet.add_subTrackGhostTrackDistSig(tvlIP.getList(reco::btau::trackGhostTrackDistSig,false));
787 >                  topjet.add_subTrackGhostTrackWeight(tvlIP.getList(reco::btau::trackGhostTrackWeight,false));
788 >                  //sv tag info
789 >                  reco::TaggingVariableList tvlSV=patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->taggingVariables();
790 >                  topjet.add_subFlightDistance2dVal(tvlSV.getList(reco::btau::flightDistance2dVal,false));
791 >                  topjet.add_subFlightDistance2dSig(tvlSV.getList(reco::btau::flightDistance2dSig,false));
792 >                  topjet.add_subFlightDistance3dVal(tvlSV.getList(reco::btau::flightDistance3dVal,false));
793 >                  topjet.add_subFlightDistance3dSig(tvlSV.getList(reco::btau::flightDistance3dSig,false));
794 >                  topjet.add_subVertexJetDeltaR(tvlSV.getList(reco::btau::vertexJetDeltaR,false));
795 >                  topjet.add_subJetNSecondaryVertices(tvlSV.get(reco::btau::jetNSecondaryVertices,-9999));
796 >                  topjet.add_subVertexNTracks(tvlSV.get(reco::btau::vertexNTracks,-9999));
797 >                  std::vector<TLorentzVector> vp4; vp4.clear();
798 >                  std::vector<float> vchi2; vchi2.clear();
799 >                  std::vector<float> vndof; vndof.clear();
800 >                  std::vector<float> vchi2ndof; vchi2ndof.clear();
801 >                  std::vector<float> vtrsize; vtrsize.clear();
802 >                  for(unsigned int i=0; i<patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->nVertices(); i++)
803 >                    {
804 >                      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > p4 = patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).p4();
805 >                      vp4.push_back(TLorentzVector(p4.px(),p4.py(),p4.pz(),p4.e()));
806 >                      vchi2.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).chi2());  
807 >                      vndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).ndof());  
808 >                      vchi2ndof.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).normalizedChi2());  
809 >                      vtrsize.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex")->secondaryVertex(i).tracksSize());  
810 >                    }
811 >                  topjet.add_subSecondaryVertex(vp4);
812 >                  topjet.add_subVertexChi2(vchi2);
813 >                  topjet.add_subVertexNdof(vndof);
814 >                  topjet.add_subVertexNormalizedChi2(vchi2ndof);
815 >                  topjet.add_subVertexTracksSize(vtrsize);
816 >                  //try computer
817 >                  edm::ESHandle<JetTagComputer> computerHandle;
818 >                  iSetup.get<JetTagComputerRecord>().get( SVComputer_.label(), computerHandle );
819 >                  const GenericMVAJetTagComputer *computer = dynamic_cast<const GenericMVAJetTagComputer*>( computerHandle.product() );
820 >                  if(computer)
821 >                    {
822 >                      computer->passEventSetup(iSetup);
823 >                      std::vector<const reco::BaseTagInfo*>  baseTagInfos;
824 >                      baseTagInfos.push_back(patsubjetd->tagInfoTrackIP("impactParameter") );
825 >                      baseTagInfos.push_back(patsubjetd->tagInfoSecondaryVertex("secondaryVertex") );      
826 >                      JetTagComputer::TagInfoHelper helper(baseTagInfos);
827 >                      reco::TaggingVariableList vars = computer->taggingVariables(helper);
828 >                      topjet.add_subVertexMassJTC(vars.get(reco::btau::vertexMass,-9999));
829 >                      topjet.add_subVertexCategoryJTC(vars.get(reco::btau::vertexCategory,-9999));
830 >                      topjet.add_subVertexEnergyRatioJTC(vars.get(reco::btau::vertexEnergyRatio,-9999));
831 >                      topjet.add_subTrackSip3dSigAboveCharmJTC(vars.get(reco::btau::trackSip3dSigAboveCharm,-9999));
832 >                    }
833 >                }
834 >           }
835 >           else
836 >             {
837 >               //filling only standard information in case the subjet has not been pat-tified during the pattuples production
838 >               subjet_v4.set_pt(pat_topjet.daughter(k)->p4().pt());
839 >               subjet_v4.set_eta(pat_topjet.daughter(k)->p4().eta());
840 >               subjet_v4.set_phi(pat_topjet.daughter(k)->p4().phi());
841 >               subjet_v4.set_energy(pat_topjet.daughter(k)->p4().E());
842 >               topjet.add_subjet(subjet_v4);
843 >             }
844 >          
845 >          
846           }
847           topjets[j].push_back(topjet);
848         }
849       }
850     }
851 <
852 <
851 >  
852 >  
853     // ------------- generator top jets -------------
854     if(doGenTopJets){
855       for(size_t j=0; j< gentopjet_sources.size(); ++j){
856        
857         gentopjets[j].clear();
858 <
858 >      
859         edm::Handle<reco::BasicJetCollection> reco_gentopjets;
860         //edm::Handle<std::vector<reco::Jet> > reco_gentopjets;
861         iEvent.getByLabel(gentopjet_sources[j], reco_gentopjets);

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines