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

Comparing UserCode/CmsHi/JetAnalysis/src/HiInclusiveJetAnalyzer.cc (file contents):
Revision 1.13 by yilmaz, Thu May 10 10:47:17 2012 UTC vs.
Revision 1.23 by yilmaz, Fri May 25 10:06:59 2012 UTC

# Line 43 | Line 43 | using namespace edm;
43   using namespace reco;
44  
45   HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) {
46 –  
46  
47    jetTag_ = iConfig.getParameter<InputTag>("jetTag");
48 +  matchTag_ = iConfig.getUntrackedParameter<InputTag>("matchTag",jetTag_);
49 +
50    vtxTag_ = iConfig.getUntrackedParameter<edm::InputTag>("vtxTag",edm::InputTag("hiSelectedVertex"));  
51 +  trackTag_ = iConfig.getParameter<InputTag>("trackTag");
52 +  useQuality_ = iConfig.getUntrackedParameter<bool>("useQuality",1);
53 +  trackQuality_ = iConfig.getUntrackedParameter<string>("trackQuality","highPurity");
54  
55    isMC_ = iConfig.getUntrackedParameter<bool>("isMC",false);
56 <  
56 >  fillGenJets_ = iConfig.getUntrackedParameter<bool>("fillGenJets",false);
57 >
58 >  doTrigger_ = iConfig.getUntrackedParameter<bool>("doTrigger",false);
59 >
60 >  rParam = iConfig.getParameter<double>("rParam");
61 >  hardPtMin_ = iConfig.getUntrackedParameter<double>("hardPtMin",4);  
62 >
63    if(isMC_){
64      genjetTag_ = iConfig.getParameter<InputTag>("genjetTag");
65      eventInfoTag_ = iConfig.getParameter<InputTag>("eventInfoTag");
# Line 67 | Line 77 | HiInclusiveJetAnalyzer::HiInclusiveJetAn
77  
78    pfCandidateLabel_ = iConfig.getUntrackedParameter<edm::InputTag>("pfCandidateLabel",edm::InputTag("particleFlowTmp"));
79  
80 <  if(!isMC_){
80 >  if(doTrigger_){
81      L1gtReadout_ = iConfig.getParameter<edm::InputTag>("L1gtReadout");
82      hltResName_ = iConfig.getUntrackedParameter<string>("hltTrgResults","TriggerResults::HLT");
83      
# Line 136 | Line 146 | HiInclusiveJetAnalyzer::beginJob() {
146    t->Branch("jty",jets_.jty,"jty[nref]/F");
147    t->Branch("jtphi",jets_.jtphi,"jtphi[nref]/F");
148    t->Branch("jtpu",jets_.jtpu,"jtpu[nref]/F");
149 +  t->Branch("jtm",jets_.jtm,"jtm[nref]/F");
150 +
151 +  // jet ID information, jet composition
152 +  t->Branch("discr_fr01", jets_.discr_fr01,"discr_fr01[nref]/F");
153 +
154 +  t->Branch("trackMax", jets_.trackMax,"trackMax[nref]/F");
155 +  t->Branch("trackSum", jets_.trackSum,"trackSum[nref]/F");
156 +  t->Branch("trackN", jets_.trackN,"trackN[nref]/I");
157 +  t->Branch("trackHardSum", jets_.trackHardSum,"trackHardSum[nref]/F");
158 +  t->Branch("trackHardN", jets_.trackHardN,"trackHardN[nref]/I");
159 +
160 +  t->Branch("chargedMax", jets_.chargedMax,"chargedMax[nref]/F");
161 +  t->Branch("chargedSum", jets_.chargedSum,"chargedSum[nref]/F");
162 +  t->Branch("chargedN", jets_.chargedN,"chargedN[nref]/I");
163 +  t->Branch("chargedHardSum", jets_.chargedHardSum,"chargedHardSum[nref]/F");
164 +  t->Branch("chargedHardN", jets_.chargedHardN,"chargedHardN[nref]/I");
165 +
166 +  t->Branch("photonMax", jets_.photonMax,"photonMax[nref]/F");
167 +  t->Branch("photonSum", jets_.photonSum,"photonSum[nref]/F");
168 +  t->Branch("photonN", jets_.photonN,"photonN[nref]/I");
169 +  t->Branch("photonHardSum", jets_.photonHardSum,"photonHardSum[nref]/F");
170 +  t->Branch("photonHardN", jets_.photonHardN,"photonHardN[nref]/I");
171 +
172 +  t->Branch("neutralMax", jets_.neutralMax,"neutralMax[nref]/F");
173 +  t->Branch("neutralSum", jets_.neutralSum,"neutralSum[nref]/F");
174 +  t->Branch("neutralN", jets_.neutralN,"neutralN[nref]/I");
175 +
176 +  t->Branch("eMax", jets_.eMax,"eMax[nref]/F");
177 +  t->Branch("eSum", jets_.eSum,"eSum[nref]/F");
178 +  t->Branch("eN", jets_.eN,"eN[nref]/I");
179 +
180 +  t->Branch("muMax", jets_.muMax,"muMax[nref]/F");
181 +  t->Branch("muSum", jets_.muSum,"muSum[nref]/F");
182 +  t->Branch("muN", jets_.muN,"muN[nref]/I");
183 +
184 +  t->Branch("matchedPt", jets_.matchedPt,"matchedPt[nref]/F");
185 +  t->Branch("matchedR", jets_.matchedR,"matchedR[nref]/F");
186  
187    // b-jet discriminators
188    if (doLifeTimeTagging_) {
# Line 183 | Line 230 | HiInclusiveJetAnalyzer::beginJob() {
230      t->Branch("muptrel", jets_.muptrel, "muptrel[nref]/F");
231      t->Branch("muchg",   jets_.muchg,   "muchg[nref]/I");
232    }
233 +
234    
235    if(isMC_){
236      t->Branch("beamId1",&jets_.beamId1,"beamId1/I");    
# Line 202 | Line 250 | HiInclusiveJetAnalyzer::beginJob() {
250      t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I");
251      t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I");
252  
253 <    // For all gen jets, matched or unmatched
254 <    t->Branch("ngen",&jets_.ngen,"ngen/I");
255 <    t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
256 <    t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
257 <    t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
258 <    t->Branch("geny",jets_.geny,"geny[ngen]/F");
259 <    t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
260 <    t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
261 <    t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
262 <
263 <    if(doSubEvent_){
264 <       t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
253 >    if(fillGenJets_){
254 >       // For all gen jets, matched or unmatched
255 >       t->Branch("ngen",&jets_.ngen,"ngen/I");
256 >       t->Branch("genmatchindex",jets_.genmatchindex,"genmatchindex[ngen]/I");
257 >       t->Branch("genpt",jets_.genpt,"genpt[ngen]/F");
258 >       t->Branch("geneta",jets_.geneta,"geneta[ngen]/F");
259 >       t->Branch("geny",jets_.geny,"geny[ngen]/F");
260 >       t->Branch("genphi",jets_.genphi,"genphi[ngen]/F");
261 >       t->Branch("gendphijt",jets_.gendphijt,"gendphijt[ngen]/F");
262 >       t->Branch("gendrjt",jets_.gendrjt,"gendrjt[ngen]/F");
263 >      
264 >       if(doSubEvent_){
265 >          t->Branch("gensubid",jets_.gensubid,"gensubid[ngen]/I");
266 >       }
267      }
268 <
268 >    
269      if(saveBfragments_  ) {
270        t->Branch("bMult",&jets_.bMult,"bMult/I");
271        t->Branch("bJetIndex",jets_.bJetIndex,"bJetIndex[bMult]/I");
# Line 303 | Line 353 | HiInclusiveJetAnalyzer::analyze(const Ev
353    
354     edm::Handle<pat::JetCollection> patjets;
355     if(usePat_)iEvent.getByLabel(jetTag_, patjets);
356 +
357 +   edm::Handle<pat::JetCollection> matchedjets;
358 +   iEvent.getByLabel(matchTag_, matchedjets);
359    
360     edm::Handle<reco::JetView> jets;
361     iEvent.getByLabel(jetTag_, jets);
362  
363     edm::Handle<reco::PFCandidateCollection> pfCandidates;
364 <   if(doLifeTimeTagging_){
365 <     iEvent.getByLabel(pfCandidateLabel_,pfCandidates);  
366 <   }
364 >   iEvent.getByLabel(pfCandidateLabel_,pfCandidates);  
365 >
366 >   edm::Handle<reco::TrackCollection> tracks;
367 >   iEvent.getByLabel(trackTag_,tracks);
368 >
369     // FILL JRA TREE
370  
371     jets_.b = b;
372     jets_.nref = 0;
373    
374 <   if(!isMC_){
374 >   if(doTrigger_){
375       fillL1Bits(iEvent);
376       fillHLTBits(iEvent);
377     }
# Line 416 | Line 471 | HiInclusiveJetAnalyzer::analyze(const Ev
471           jets_.mudr[jets_.nref]    =  reco::deltaR(jet,muon);
472           jets_.muptrel[jets_.nref] =  getPtRel(muon, jet);
473           jets_.muchg[jets_.nref]   =  muon.charge();
474 +       }else{
475 +         jets_.mupt[jets_.nref]    =  0.0;
476 +         jets_.mueta[jets_.nref]   =  0.0;
477 +         jets_.muphi[jets_.nref]   =  0.0;
478 +         jets_.mue[jets_.nref]     =  0.0;
479 +         jets_.mudr[jets_.nref]    =  9.9;
480 +         jets_.muptrel[jets_.nref] =  0.0;
481 +         jets_.muchg[jets_.nref]   = 0;
482         }
483 +
484       }
485  
486 +
487 +
488 +     // Jet ID variables
489 +
490 +     jets_.muMax[jets_.nref] = 0;
491 +     jets_.muSum[jets_.nref] = 0;
492 +     jets_.muN[jets_.nref] = 0;
493 +
494 +     jets_.eMax[jets_.nref] = 0;
495 +     jets_.eSum[jets_.nref] = 0;
496 +     jets_.eN[jets_.nref] = 0;
497 +
498 +     jets_.neutralMax[jets_.nref] = 0;
499 +     jets_.neutralSum[jets_.nref] = 0;
500 +     jets_.neutralN[jets_.nref] = 0;
501 +
502 +     jets_.photonMax[jets_.nref] = 0;
503 +     jets_.photonSum[jets_.nref] = 0;
504 +     jets_.photonN[jets_.nref] = 0;
505 +     jets_.photonHardSum[jets_.nref] = 0;
506 +     jets_.photonHardN[jets_.nref] = 0;
507 +
508 +     jets_.chargedMax[jets_.nref] = 0;
509 +     jets_.chargedSum[jets_.nref] = 0;
510 +     jets_.chargedN[jets_.nref] = 0;
511 +     jets_.chargedHardSum[jets_.nref] = 0;
512 +     jets_.chargedHardN[jets_.nref] = 0;
513 +
514 +     jets_.trackMax[jets_.nref] = 0;
515 +     jets_.trackSum[jets_.nref] = 0;
516 +     jets_.trackN[jets_.nref] = 0;
517 +     jets_.trackHardSum[jets_.nref] = 0;
518 +     jets_.trackHardN[jets_.nref] = 0;
519 +
520 +
521 +     for(unsigned int icand = 0; icand < tracks->size(); ++icand){
522 +        const reco::Track& track = (*tracks)[icand];
523 +        if(useQuality_ ){
524 +           bool goodtrack = track.quality(reco::TrackBase::qualityByName(trackQuality_));
525 +           if(!goodtrack) continue;
526 +        }
527 +
528 +        double dr = deltaR(jet,track);
529 +        if(dr < rParam){
530 +           double ptcand = track.pt();
531 +           jets_.trackSum[jets_.nref] += ptcand;
532 +           jets_.trackN[jets_.nref] += 1;
533 +
534 +           if(ptcand > hardPtMin_){
535 +              jets_.trackHardSum[jets_.nref] += ptcand;
536 +              jets_.trackHardN[jets_.nref] += 1;
537 +
538 +           }
539 +           if(ptcand > jets_.trackMax[jets_.nref]) jets_.trackMax[jets_.nref] = ptcand;
540 +
541 +        }
542 +     }
543 +
544 +     for(unsigned int icand = 0; icand < pfCandidates->size(); ++icand){
545 +        const reco::PFCandidate& track = (*pfCandidates)[icand];
546 +        double dr = deltaR(jet,track);
547 +        if(dr < rParam){
548 +           double ptcand = track.pt();
549 +           int pfid = track.particleId();
550 +
551 +           switch(pfid){
552 +
553 +           case 1:
554 +              jets_.chargedSum[jets_.nref] += ptcand;
555 +              jets_.chargedN[jets_.nref] += 1;
556 +              if(ptcand > hardPtMin_){
557 +                 jets_.chargedHardSum[jets_.nref] += ptcand;
558 +                 jets_.chargedHardN[jets_.nref] += 1;
559 +              }
560 +              if(ptcand > jets_.chargedMax[jets_.nref]) jets_.chargedMax[jets_.nref] = ptcand;
561 +              break;
562 +
563 +           case 2:
564 +              jets_.eSum[jets_.nref] += ptcand;
565 +              jets_.eN[jets_.nref] += 1;
566 +              if(ptcand > jets_.eMax[jets_.nref]) jets_.eMax[jets_.nref] = ptcand;
567 +              break;
568 +
569 +           case 3:
570 +              jets_.muSum[jets_.nref] += ptcand;
571 +              jets_.muN[jets_.nref] += 1;
572 +              if(ptcand > jets_.muMax[jets_.nref]) jets_.muMax[jets_.nref] = ptcand;
573 +              break;
574 +
575 +           case 4:
576 +              jets_.photonSum[jets_.nref] += ptcand;
577 +              jets_.photonN[jets_.nref] += 1;
578 +              if(ptcand > hardPtMin_){
579 +                 jets_.photonHardSum[jets_.nref] += ptcand;
580 +                 jets_.photonHardN[jets_.nref] += 1;
581 +              }
582 +              if(ptcand > jets_.photonMax[jets_.nref]) jets_.photonMax[jets_.nref] = ptcand;
583 +              break;
584 +
585 +           case 5:
586 +              jets_.neutralSum[jets_.nref] += ptcand;
587 +              jets_.neutralN[jets_.nref] += 1;
588 +              if(ptcand > jets_.neutralMax[jets_.nref]) jets_.neutralMax[jets_.nref] = ptcand;
589 +              break;
590 +
591 +           default:
592 +             break;
593 +
594 +           }
595 +        }
596 +     }
597 +
598 +     double drMin = 100;
599 +     for(unsigned int imatch = 0 ; imatch < matchedjets->size(); ++imatch){
600 +        const reco::Jet& mjet = (*matchedjets)[imatch];
601 +
602 +        double dr = deltaR(jet,mjet);
603 +        if(dr < drMin){
604 +           jets_.matchedPt[jets_.nref] = mjet.pt();
605 +           jets_.matchedR[jets_.nref] = dr;
606 +           drMin = dr;
607 +        }
608 +     }
609 +  
610 +
611 +
612 +
613 +     //     if(etrk.quality(reco::TrackBase::qualityByName(qualityString_))) pev_.trkQual[pev_.nTrk]=1;
614 +
615 +        /////////////////////////////////////////////////////////////////
616 +        // Jet core pt^2 discriminant for fake jets
617 +        // Edited by Yue Shi Lai <ylai@mit.edu>
618 +
619 +        // Initial value is 0
620 +        jets_.discr_fr01[jets_.nref] = 0;
621 +        // Start with no directional adaption, i.e. the fake rejection
622 +        // axis is the jet axis
623 +        float pseudorapidity_adapt = jets_.jteta[jets_.nref];
624 +        float azimuth_adapt = jets_.jtphi[jets_.nref];
625 +
626 +        // Unadapted discriminant with adaption search
627 +        for (size_t iteration = 0; iteration < 2; iteration++) {
628 +                float pseudorapidity_adapt_new = pseudorapidity_adapt;
629 +                float azimuth_adapt_new = azimuth_adapt;
630 +                float max_weighted_perp = 0;
631 +                float perp_square_sum = 0;
632 +
633 +                for (size_t index_pf_candidate = 0;
634 +                         index_pf_candidate < pfCandidates->size();
635 +                         index_pf_candidate++) {
636 +                        const reco::PFCandidate &p =
637 +                                (*pfCandidates)[index_pf_candidate];
638 +
639 +                        switch (p.particleId()) {
640 +                          //case 1:     // Charged hadron
641 +                          //case 3:     // Muon
642 +                        case 4: // Photon
643 +                                {
644 +                                        const float dpseudorapidity =
645 +                                                p.eta() - pseudorapidity_adapt;
646 +                                        const float dazimuth =
647 +                                                reco::deltaPhi(p.phi(), azimuth_adapt);
648 +                                        // The Gaussian scale factor is 0.5 / (0.1 * 0.1)
649 +                                        // = 50
650 +                                        const float angular_weight =
651 +                                                exp(-50.0F * (dpseudorapidity * dpseudorapidity +
652 +                                                                          dazimuth * dazimuth));
653 +                                        const float weighted_perp =
654 +                                                angular_weight * p.pt() * p.pt();
655 +                                        const float weighted_perp_square =
656 +                                                weighted_perp * p.pt();
657 +
658 +                                        perp_square_sum += weighted_perp_square;
659 +                                        if (weighted_perp >= max_weighted_perp) {
660 +                                                pseudorapidity_adapt_new = p.eta();
661 +                                                azimuth_adapt_new = p.phi();
662 +                                                max_weighted_perp = weighted_perp;
663 +                                        }
664 +                                }
665 +                        default:                                                  
666 +                          break;
667 +                        }
668 +                }
669 +                // Update the fake rejection value
670 +                jets_.discr_fr01[jets_.nref] = std::max(
671 +                        jets_.discr_fr01[jets_.nref], perp_square_sum);
672 +                // Update the directional adaption
673 +                pseudorapidity_adapt = pseudorapidity_adapt_new;
674 +                azimuth_adapt = azimuth_adapt_new;
675 +        }
676 +
677 +
678       jets_.jtpt[jets_.nref] = jet.pt();                            
679       jets_.jteta[jets_.nref] = jet.eta();
680       jets_.jtphi[jets_.nref] = jet.phi();
681       jets_.jty[jets_.nref] = jet.eta();
682       jets_.jtpu[jets_.nref] = jet.pileup();
683 +     jets_.jtm[jets_.nref] = jet.mass();
684          
685       if(isMC_ && usePat_){
686  

Diff Legend

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