ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc
(Generate patch)

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.22 by yygao, Mon Mar 29 15:27:24 2010 UTC vs.
Revision 1.23 by yygao, Mon May 31 14:19:45 2010 UTC

# Line 86 | Line 86 | PVStudy::PVStudy(const edm::ParameterSet
86    nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
87    nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
88    zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
89 –  fHighPuritycut_            = iConfig.getUntrackedParameter<double>("fHighPuritycut");  
89    useHWTrk_                  = iConfig.getUntrackedParameter<bool>("useHWTrk",false);
90 +  avgTrkPtInPVMin_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMin");
91 +  avgTrkPtInPVMax_                  = iConfig.getUntrackedParameter<double>("avgTrkPtInPVMax");
92    bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
93    
94    //now do what ever initialization is needed
# Line 131 | Line 132 | PVStudy::PVStudy(const edm::ParameterSet
132      //Fill the variables in the twovtx pair (recvtx1, recvtx2)  
133      // Information related to the analyzing the two-vertex method
134  
134 –
135 –
135      pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
136      pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
137      pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
# Line 221 | Line 220 | PVStudy::PVStudy(const edm::ParameterSet
220        pvtxtree_->Branch("ndofPV_spl2_mct",&ndofPV_spl2_mct_,"ndofPV_spl2_mct[nrecPV_spl2_mct]/D");
221        pvtxtree_->Branch("normchi2PV_spl2_mct",&normchi2PV_spl2_mct_,"normchi2PV_spl2_mct[nrecPV_spl2_mct]/D");
222        pvtxtree_->Branch("avgPtPV_spl2_mct",&avgPtPV_spl2_mct_,"avgPtPV_spl2_mct[nrecPV_spl2_mct]/D");
224 –
223      }
224    }
225    
# Line 645 | Line 643 | PVStudy::analyze(const edm::Event& iEven
643    // Step 1:  Apply event cleaning for data and MC
644    //          WARNING: event selection cut are hard coded!!
645    //========================================================================
646 <  // =====Loop over the trackColl to get the fraction of HighPurity tracks
649 <  int nTracks = 0;
650 <  int nHighPurityTracks=0;
651 <  double fHighPurity=0;
652 <
653 <  for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
654 <    TrackRef tkref(trackCollectionHandle,i);
655 <    if(tkref->quality(reco::TrackBase::highPurity)) ++nHighPurityTracks;
656 <  }
657 <  if(nTracks>0)
658 <    fHighPurity =  double(nHighPurityTracks)/double(nTracks);
659 <  if(verbose_)
660 <    cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;  
661 <    
662 <  if( (fHighPurity<fHighPuritycut_ && nTracks>10) || vertexColl->begin()->isFake())  return;
663 <  
646 >
647    glob_runno_ = iEvent.id().run();
648    glob_evtno_ = iEvent.id().event();
649    glob_ls_   = iEvent.luminosityBlock();
# Line 723 | Line 706 | PVStudy::analyze(const edm::Event& iEven
706      for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
707          v!=vertexColl->end(); ++v) {  
708        if(v->isValid() && ! v->isFake()) {
726 –        // xsep
709          if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
710            min_xsep = fabs(v->x()- vertexColl->begin()->x());
729 –        // ysep
711          if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
712            min_ysep = fabs(v->y()- vertexColl->begin()->y());
732 –        // zsep
713          if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
714            min_zsep = fabs(v->z()- vertexColl->begin()->z());
735 –        // xsep_sign
715          double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
716          if(xsep_sign < min_xsep_sign)
717            min_xsep_sign = xsep_sign;
739 –        // ysep_sign
718          double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
719          if(ysep_sign < min_ysep_sign)
720            min_ysep_sign = ysep_sign;
743 –        // zsep_sign
721          double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
722          if(zsep_sign < min_zsep_sign)
723            min_zsep_sign = zsep_sign;
747 –        // ntrksep
724          if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
725            min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
750 –        // sumpt2sep
726          if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
727            min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
728        }
# Line 856 | Line 831 | PVStudy::analyze(const edm::Event& iEven
831        simy_[isimpv] = vsim->y;
832        simz_[isimpv] = vsim->y;
833        simptsq_[isimpv] = vsim->ptsq;
834 <      //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
860 <      fillMCHisto(vsim, isimpv, vertexColl, 3);
834 >      fillMCHisto(vsim, isimpv, vertexColl, 3, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
835      }
836    }
837    
# Line 872 | Line 846 | PVStudy::analyze(const edm::Event& iEven
846        int isimpv = 0;
847        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
848            vsim!=simpv.end(); vsim++, isimpv++){
849 <        //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
850 <        fillMCHisto(vsim, isimpv, matchedVertexColl1, 1);
877 <        fillMCHisto(vsim, isimpv, matchedVertexColl2, 2);
849 >        fillMCHisto(vsim, isimpv, matchedVertexColl1, 1, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
850 >        fillMCHisto(vsim, isimpv, matchedVertexColl2, 2, avgTrkPtInPVMin_, avgTrkPtInPVMax_);
851        }
852      }
853  
# Line 889 | Line 862 | PVStudy::analyze(const edm::Event& iEven
862      nrecPV_twovtx_ = 0;
863      for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
864          v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
865 <        ++v1, ++v2) {
865 >        ++v1, ++v2) {
866 >      
867 >      // ==================================================================
868 >      // With the option to fill the histograms at a given average pT range
869 >      // ==================================================================
870 >      if (avgPtRecVtx(*v1) < avgTrkPtInPVMin_ || avgPtRecVtx(*v1) > avgTrkPtInPVMax_ ) continue;
871 >      if (avgPtRecVtx(*v2) < avgTrkPtInPVMin_ || avgPtRecVtx(*v2) > avgTrkPtInPVMax_ ) continue;
872 >
873        h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
874        h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
875        fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
# Line 910 | Line 890 | PVStudy::analyze(const edm::Event& iEven
890        }
891        
892        fntrk_ = (nTrkPV1 + nTrkPV2)/2;
893 <
893 >      
894        h_summary->Fill1d("deltax", fres_[0] );
895        h_summary->Fill1d("deltay", fres_[1] );
896        h_summary->Fill1d("deltaz", fres_[2] );
# Line 924 | Line 904 | PVStudy::analyze(const edm::Event& iEven
904                                         error(ferror_[0], ferror_[1], ferror_[2]),
905                                         fntrk_));
906        // Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
907 +
908        fillHisto(res(fres_[0], fres_[1], fres_[2]),
909                  error(ferror_[0], ferror_[1], ferror_[2]),
910                  fntrk_,0);  
# Line 1144 | Line 1125 | void PVStudy::printSimTrks(const Handle<
1125    }
1126   }
1127  
1128 < void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
1128 > void PVStudy::fillHisto(res r, error er, int ntk, int datamode) {
1129    stringstream ss;
1130    ss << ntk;
1131    string suffix;
# Line 1334 | Line 1315 | void PVStudy::fillTrackHistoInPV(const r
1315    }
1316   }
1317  
1318 < void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
1318 > void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl,
1319 >                          int datatype, double avgTrkPtInPVMin_, double avgTrkPtInPVMax_)
1320   {
1321    std::string suffix;
1322    // Set the vertexColl and histogram suffix according to datamode
# Line 1375 | Line 1357 | void PVStudy::fillMCHisto(std::vector<si
1357    // === If match found fill the residual and pulls
1358    if(vsim->recVtx) {
1359      edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
1360 +    double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1361 +    if(avgPtPV>avgTrkPtInPVMax_ || avgPtPV<avgTrkPtInPVMin_) return;
1362 +
1363      // if study the resolution/pull for all tracks in the PVtx
1364      int nrectrk;
1380 –    
1365      if(useHWTrk_)
1366        nrectrk=nHWTrkRecVtx(*(vsim->recVtx));
1367      else
# Line 1388 | Line 1372 | void PVStudy::fillMCHisto(std::vector<si
1372      
1373      double ndofPV = vsim->recVtx->ndof();
1374      double normchi2PV =  vsim->recVtx->normalizedChi2();
1375 <    double avgPtPV = avgPtRecVtx(*(vsim->recVtx));
1392 <    
1375 >
1376      fres_mct[0] = vsim->recVtx->x()-vsim->x;
1377      fres_mct[1] = vsim->recVtx->y()-vsim->y;
1378      fres_mct[2] = vsim->recVtx->z()-vsim->z;

Diff Legend

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