[ViewVC] Diff of: cvsroot/UserCode/Jeng/PVStudy/plugins/PVStudy.cc

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.3 by jengbou, Tue Sep 8 07:15:38 2009 UTC vs.
Revision 1.17 by yygao, Tue Jan 26 12:06:49 2010 UTC

+*/
+//
+// Original Author:  "Geng-yuan Jeng/UC Riverside"
-+
+//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
+//         Created:  Thu Aug 20 11:55:40 CDT 2009
+// $Id$
+//
+#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
+#include <SimDataFormats/Track/interface/SimTrack.h>
+#include <SimDataFormats/Track/interface/SimTrackContainer.h>
-+
+// BeamSpot
-+
+#include "DataFormats/BeamSpot/interface/BeamSpot.h"
-+
+// Trigger
-+
+#include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
-+
+#include "CondFormats/L1TObjects/interface/L1GtPrescaleFactors.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsAlgoTrigRcd.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtPrescaleFactorsTechTrigRcd.h"
-+
+#include "CondFormats/L1TObjects/interface/L1GtTriggerMask.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtTriggerMaskAlgoTrigRcd.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtTriggerMaskTechTrigRcd.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoAlgoTrigRcd.h"
-+
+#include "CondFormats/DataRecord/interface/L1GtTriggerMaskVetoTechTrigRcd.h"
-+
+#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
-+
+#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetup.h"
-+
+#include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
-–
+//root
-–
+#include "TROOT.h"
-–
+#include "TF1.h"
-–
+#include "TString.h"
-–
+#include "TStyle.h"
-+
+//root
-+
+#include <TROOT.h>
-+
+#include <TF1.h>
-+
+#include <TString.h>
-+
+#include <TStyle.h>
-+
+#include <TPaveStats.h>
-+
+#include <TPad.h>
+using namespace std;
-+
+typedef math::XYZTLorentzVectorF LorentzVector;
-+
+typedef math::XYZPoint Point;
-<
+PVStudy::PVStudy(const edm::ParameterSet& iConfig):
-<
+  nTrkMin_(10),nTrkMax_(100)
->
+PVStudy::PVStudy(const edm::ParameterSet& iConfig)
-<
+  simG4_           = iConfig.getParameter<edm::InputTag>( "simG4" );
-<
+  tracksLabel_     = iConfig.getUntrackedParameter<edm::InputTag>("tracks");
-<
+  vtxSample_       = iConfig.getUntrackedParameter<edm::InputTag>("vtxSample");
-<
+  verbose_         = iConfig.getUntrackedParameter<bool>("verbose", false);
-<
+  realData_        = iConfig.getUntrackedParameter<bool>("realData",false);
-<
+  analyze_         = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
-<
+  outputfilename_  = iConfig.getUntrackedParameter<string>("OutputFileName");
-<
->
+  //=======================================================================
->
+  // Get configuration for TrackTupleMaker
->
+  //=======================================================================
->
+  simG4_                     = iConfig.getParameter<edm::InputTag>( "simG4" );
->
+  trackCollectionTag_        = iConfig.getUntrackedParameter<edm::InputTag>("trackCollection");
->
+  splitTrackCollection1Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection1");
->
+  splitTrackCollection2Tag_  = iConfig.getUntrackedParameter<edm::InputTag>("splitTrackCollection2");
->
+  vertexCollectionTag_       = iConfig.getUntrackedParameter<edm::InputTag>("vertexCollection");
->
+  splitVertexCollection1Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection1");
->
+  splitVertexCollection2Tag_ = iConfig.getUntrackedParameter<edm::InputTag>("splitVertexCollection2");
->
+  verbose_                   = iConfig.getUntrackedParameter<bool>("verbose",false);
->
+  realData_                  = iConfig.getUntrackedParameter<bool>("realData",false);
->
+  analyze_                   = iConfig.getUntrackedParameter<bool>("analyzeOnTheFly",false);
->
+  saventuple_                = iConfig.getUntrackedParameter<bool>("saventuple",false);
->
+  outputfilename_            = iConfig.getUntrackedParameter<string>("OutputFileName");
->
+  ntrkdiffcut_               = iConfig.getUntrackedParameter<double>("ntrkdiffcut");
->
+  nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
->
+  nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
->
+  zsigncut_                  = iConfig.getUntrackedParameter<int>("zsigncut");
->
+  bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
->
+  applyEvtClean_            = iConfig.getUntrackedParameter<bool>("applyEvtClean",false);
->
+  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
->
+  //now do what ever initialization is needed
+  pvinfo.clear();
-–
+  file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
-–
+  ftree_ = new TTree("mytree","mytree");
-–
+  ftree_->AutoSave();
-–
+  ftree_->Branch("residual",&fres,"fres[3]/D");
-–
+  ftree_->Branch("error",&ferror,"ferror[3]/D");
-–
+  ftree_->Branch("nTrk",&fntrk,"fntrk/I");
-–
-–
+  edm::Service<TFileService> fs;
-–
+//   TFileDirectory subDir = fs->mkdir( "Summary" );
-–
+//   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );
-–
-–
+  h_nTrk           = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
-–
+  h_nTrkPV         = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
-–
+  h_trkPt          = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
-–
+  h_trkPtPV        = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);
-+
+  // Specify the data mode vector
-+
+  if(realData_) datamode.push_back(0);
-+
+  else {
-+
+    datamode.push_back(0);
-+
+    datamode.push_back(1);
-+
+    datamode.push_back(2);
-+
+    datamode.push_back(3);
-+
+  }
-+
+ // Create ntuple files if needed
-+
+  if(saventuple_) {
-+
+    file_ = TFile::Open(outputfilename_.c_str(),"RECREATE");
-+
+    ftree_ = new TTree("mytree","mytree");
-+
+    ftree_->AutoSave();
-+
+    ftree_->Branch("residual",&fres_,"fres_[3]/D");
-+
+    ftree_->Branch("error",&ferror_,"ferror_[3]/D");
-+
+    ftree_->Branch("nTrkPV",&fntrk_,"fntrk_/I");
-+
-+
+    // pvtxtree_ analyzing the pvtxs ootb
-+
+    pvtxtree_ = new TTree("pvtxtree","pvtxtree");
-+
+    //pvtx using all tracks
-+
+    pvtxtree_->Branch("nrecPV",&nrecPV_,"nrecPV/I");
-+
+    pvtxtree_->Branch("nTrkPV",&nTrkPV_,"nTrkPV[nrecPV]/I");
-+
+    pvtxtree_->Branch("sumptsq",&sumptsq_,"sumptsq[nrecPV]/D");
-+
+    pvtxtree_->Branch("isValid",&isValid_,"isValid/I");
-+
+    pvtxtree_->Branch("isFake",&isFake_,"isFake/I");
-+
+    pvtxtree_->Branch("recx",&recx_,"recx[nrecPV]/D");
-+
+    pvtxtree_->Branch("recy",&recy_,"recy[nrecPV]/D");
-+
+    pvtxtree_->Branch("recz",&recz_,"recz[nrecPV]/D");
-+
+    pvtxtree_->Branch("recx_err",&recx_err_,"recx_err[nrecPV]/D");
-+
+    pvtxtree_->Branch("recy_err",&recy_err_,"recy_err[nrecPV]/D");
-+
+    pvtxtree_->Branch("recz_err",&recz_err_,"recz_err[nrecPV]/D");
-+
-+
+    pvtxtree_->Branch("min_zsep",&min_zsep_,"min_zsep/D");
-+
+    pvtxtree_->Branch("min_ntrksep",&min_ntrksep_,"min_ntrksep/D");
-+
+    pvtxtree_->Branch("min_sumpt2sep",&min_sumpt2sep_,"min_sumpt2sep/D");
-+
-+
+    //pvtx using splittracks1
-+
+    pvtxtree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl/I");
-+
+    pvtxtree_->Branch("nTrkPV1_spl",&nTrkPV1_spl_,"nTrkPV1_spl[nrecPV1_spl]/I");
-+
+    pvtxtree_->Branch("sumptsq1_spl",&sumptsq1_spl_,"sumptsq1_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("isValid1_spl",&isValid1_spl_,"isValid1_spl/I");
-+
+    pvtxtree_->Branch("isFake1_spl",&isFake1_spl_,"isFake1_spl/I");
-+
+    pvtxtree_->Branch("recx1_spl",&recx1_spl_,"recx1_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("recy1_spl",&recy1_spl_,"recy1_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("recz1_spl",&recz1_spl_,"recz1_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("recx1_err_spl",&recx1_err_spl_,"recx1_err_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("recy1_err_spl",&recy1_err_spl_,"recy1_err_spl[nrecPV1_spl]/D");
-+
+    pvtxtree_->Branch("recz1_err_spl",&recz1_err_spl_,"recz1_err_spl[nrecPV1_spl]/D");
-+
-+
+    //pvtx using splittracks2
-+
+    pvtxtree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl/I");
-+
+    pvtxtree_->Branch("nTrkPV2_spl",&nTrkPV2_spl_,"nTrkPV2_spl[nrecPV2_spl]/I");
-+
+    pvtxtree_->Branch("sumptsq2_spl",&sumptsq2_spl_,"sumptsq2_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("isValid2_spl",&isValid2_spl_,"isValid2_spl/I");
-+
+    pvtxtree_->Branch("isFake2_spl",&isFake2_spl_,"isFake2_spl/I");
-+
+    pvtxtree_->Branch("recx2_spl",&recx2_spl_,"recx2_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("recy2_spl",&recy2_spl_,"recy2_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("recz2_spl",&recz2_spl_,"recz2_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("recx2_err_spl",&recx2_err_spl_,"recx2_err_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("recy2_err_spl",&recy2_err_spl_,"recy2_err_spl[nrecPV2_spl]/D");
-+
+    pvtxtree_->Branch("recz2_err_spl",&recz2_err_spl_,"recz2_err_spl[nrecPV2_spl]/D");
-+
-+
+    //pixeVertices
-+
+    pvtxtree_->Branch("nrecPV_pxlpvtx",&nrecPV_pxlpvtx_,"nrecPV_pxlpvtx/I");
-+
+    pvtxtree_->Branch("nTrkPV_pxlpvtx",&nTrkPV_pxlpvtx_,"nTrkPV_pxlpvtx[nrecPV_pxlpvtx]/I");
-+
+    pvtxtree_->Branch("sumptsq_pxlpvtx",&sumptsq_pxlpvtx_,"sumptsq_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("isValid_pxlpvtx",&isValid_pxlpvtx_,"isValid_pxlpvtx/I");
-+
+    pvtxtree_->Branch("isFake_pxlpvtx",&isFake_pxlpvtx_,"isFake_pxlpvtx/I");
-+
+    pvtxtree_->Branch("recx_pxlpvtx",&recx_pxlpvtx_,"recx_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("recy_pxlpvtx",&recy_pxlpvtx_,"recy_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("recz_pxlpvtx",&recz_pxlpvtx_,"recz_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("recx_err_pxlpvtx",&recx_err_pxlpvtx_,"recx_err_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("recy_err_pxlpvtx",&recy_err_pxlpvtx_,"recy_err_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
+    pvtxtree_->Branch("recz_err_pxlpvtx",&recz_err_pxlpvtx_,"recz_err_pxlpvtx[nrecPV_pxlpvtx]/D");
-+
-+
+    //Fill the variables in the twovtx pair (recvtx1, recvtx2)
-+
+    pvtxtree_->Branch("nrecPV_twovtx",&nrecPV_twovtx_,"nrecPV_twovtx/I");
-+
+    pvtxtree_->Branch("nTrkPV1_spl_twovtx",&nTrkPV1_spl_twovtx_,"nTrkPV1_spl_twovtx[nrecPV_twovtx]/I");
-+
+    pvtxtree_->Branch("nTrkPV2_spl_twovtx",&nTrkPV2_spl_twovtx_,"nTrkPV2_spl_twovtx[nrecPV_twovtx]/I");
-+
+    pvtxtree_->Branch("sumptsq1_spl_twovtx",&sumptsq1_spl_twovtx_,"sumptsq1_spl_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("sumptsq2_spl_twovtx",&sumptsq2_spl_twovtx_,"sumptsq2_spl_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("nTrkPV_twovtx",&nTrkPV_twovtx_,"nTrkPV_twovtx[nrecPV_twovtx]/I");
-+
+    pvtxtree_->Branch("deltax_twovtx",&deltax_twovtx_,"deltax_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("deltay_twovtx",&deltay_twovtx_,"deltay_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("deltaz_twovtx",&deltaz_twovtx_,"deltaz_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("errx_twovtx",&errx_twovtx_,"errx_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("erry_twovtx",&erry_twovtx_,"erry_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("errz_twovtx",&errz_twovtx_,"errz_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("pullx_twovtx",&pullx_twovtx_,"pullx_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("pully_twovtx",&pully_twovtx_,"pully_twovtx[nrecPV_twovtx]/D");
-+
+    pvtxtree_->Branch("pullz_twovtx",&pullz_twovtx_,"pullz_twovtx[nrecPV_twovtx]/D");
-+
-+
+    // MC variables
-+
+    if(!realData_) {
-+
+      pvtxtree_->Branch("nsimPV",&nsimPV_,"nsimPV/I");
-+
+      pvtxtree_->Branch("nsimTrkPV",&nsimTrkPV_,"nsimTrkPV[nsimPV]/I");
-+
+      pvtxtree_->Branch("simx",&simx_,"simx[nsimPV]/D");
-+
+      pvtxtree_->Branch("simy",&simy_,"simy[nsimPV]/D");
-+
+      pvtxtree_->Branch("simz",&simz_,"simz[nsimPV]/D");
-+
+      pvtxtree_->Branch("simptsq",&simptsq_,"simptsq[nsimPV]/D");
-+
-+
+      // For pvtxs with all the tracks
-+
+      pvtxtree_->Branch("nrecPV_mct",&nrecPV_mct_,"nrecPV_mct/I");
-+
+      pvtxtree_->Branch("deltax_mct",&deltax_mct_,"deltax_mct[nrecPV_mct]/D");
-+
+      pvtxtree_->Branch("deltay_mct",&deltay_mct_,"deltay_mct[nrecPV_mct]/D");
-+
+      pvtxtree_->Branch("deltaz_mct",&deltaz_mct_,"deltaz_mct[nrecPV_mct]/D");
-+
+      pvtxtree_->Branch("pullx_mct",&pullx_mct_,"pullx_mct[nrecPV_mct]/D");
-+
+      pvtxtree_->Branch("pully_mct",&pully_mct_,"pully_mct[nrecPV_mct]/D");
-+
+      pvtxtree_->Branch("pullz_mct",&pullz_mct_,"pullz_mct[nrecPV_mct]/D");
-+
+      // For pvtxs with splittracks1
-+
+      pvtxtree_->Branch("nrecPV_spl1_mct",&nrecPV_spl1_mct_,"nrecPV_spl1_mct/I");
-+
+      pvtxtree_->Branch("deltax_spl1_mct",&deltax_spl1_mct_,"deltax_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      pvtxtree_->Branch("deltay_spl1_mct",&deltay_spl1_mct_,"deltay_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      pvtxtree_->Branch("deltaz_spl1_mct",&deltaz_spl1_mct_,"deltaz_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      pvtxtree_->Branch("pullx_spl1_mct",&pullx_spl1_mct_,"pullx_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      pvtxtree_->Branch("pully_spl1_mct",&pully_spl1_mct_,"pully_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      pvtxtree_->Branch("pullz_spl1_mct",&pullz_spl1_mct_,"pullz_spl1_mct[nrecPV_spl1_mct]/D");
-+
+      // For pvtxs with splittracks1
-+
+      pvtxtree_->Branch("nrecPV_spl2_mct",&nrecPV_spl2_mct_,"nrecPV_spl2_mct/I");
-+
+      pvtxtree_->Branch("deltax_spl2_mct",&deltax_spl2_mct_,"deltax_spl2_mct[nrecPV_spl2_mct]/D");
-+
+      pvtxtree_->Branch("deltay_spl2_mct",&deltay_spl2_mct_,"deltay_spl2_mct[nrecPV_spl2_mct]/D");
-+
+      pvtxtree_->Branch("deltaz_spl2_mct",&deltaz_spl2_mct_,"deltaz_spl2_mct[nrecPV_spl2_mct]/D");
-+
+      pvtxtree_->Branch("pullx_spl2_mct",&pullx_spl2_mct_,"pullx_spl2_mct[nrecPV_spl2_mct]/D");
-+
+      pvtxtree_->Branch("pully_spl2_mct",&pully_spl2_mct_,"pully_spl2_mct[nrecPV_spl2_mct]/D");
-+
+      pvtxtree_->Branch("pullz_spl2_mct",&pullz_spl2_mct_,"pullz_spl2_mct[nrecPV_spl2_mct]/D");
-+
+    }
-+
+  }
-+
-+
+  setRootStyle();
-+
-+
+  //========================================
-+
+  // Booking histograms
-+
+  //========================================
-+
-+
+  // Create a root file for histograms
-+
+  theFile = new TFile(histoFileName_.c_str(), "RECREATE");
-+
+  // make diretories
-+
+  theFile->mkdir("Summary");
-+
+  theFile->cd();
-+
+  theFile->mkdir("Others");
-+
+  theFile->cd();
-+
+  if (analyze_) {
-+
+    theFile->mkdir("Results");
-+
+    theFile->cd();
-+
+  }
-+
-+
+  // Book Histograms:
-+
+  h_pvtrk   = new PVHistograms();
-+
+  h_pixvtx  = new PVHistograms();
-+
+  h_misc    = new PVHistograms();
-+
+  h_summary = new PVHistograms();
-+
+  h_others  = new PVHistograms();
-+
-+
+  for(int i=0;i<3;i++) {
-+
+    if(i == 0) h_pvtrk->Init("pvTrk");
-+
+    else {
-+
+      stringstream ss;
-+
+      ss << i;
-+
+      h_pvtrk->Init("pvTrk", ss.str(),"spl");
-+
+    }
-+
+  }
-+
+  h_pixvtx->Init("pixVtx");
-+
+  h_misc->Init("misc");
-+
-+
+  // Book MC only plots
+  if (!realData_) {
-<
+    h_deltax       = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
-<
+    h_deltay       = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
-<
+    h_deltaz       = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
-<
+    h_pullx        = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
-<
+    h_pully        = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
-<
+    h_pullz        = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
-<
-<
+    // Summary of analysis:
-<
+    if (analyze_) {
-<
+      h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
-<
+      h_resx_nTrk->SetMarkerStyle(21);
-<
+      h_resx_nTrk->SetMarkerColor(4);
-<
+      h_resx_nTrk->GetXaxis()->SetTitle("Num of tracks");
-<
+      h_resx_nTrk->GetYaxis()->SetTitle("#mum");
-<
+      h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
-<
+      h_resy_nTrk->SetMarkerStyle(21);
-<
+      h_resy_nTrk->SetMarkerColor(4);
-<
+      h_resy_nTrk->GetXaxis()->SetTitle("Num of tracks");
-<
+      h_resy_nTrk->GetYaxis()->SetTitle("#mum");
-<
+      h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
-<
+      h_resz_nTrk->SetMarkerStyle(21);
-<
+      h_resz_nTrk->SetMarkerColor(4);
-<
+      h_resz_nTrk->GetXaxis()->SetTitle("Num of tracks");
-<
+      h_resz_nTrk->GetYaxis()->SetTitle("#mum");
-<
+      h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
-<
+      h_pullx_nTrk->SetMarkerStyle(21);
-<
+      h_pullx_nTrk->SetMarkerColor(4);
-<
+      h_pullx_nTrk->SetBit(TH1::kCanRebin);
-<
+      h_pullx_nTrk->GetXaxis()->SetTitle("Num of tracks");
-<
+      h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
-<
+      h_pully_nTrk->SetMarkerStyle(21);
-<
+      h_pully_nTrk->SetMarkerColor(4);
-<
+      h_pully_nTrk->SetBit(TH1::kCanRebin);
-<
+      h_pully_nTrk->GetXaxis()->SetTitle("Num of tracks");
-<
+      h_pullz_nTrk   = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
-<
+      h_pullz_nTrk->SetMarkerStyle(21);
-<
+      h_pullz_nTrk->SetMarkerColor(4);
-<
+      h_pullz_nTrk->SetBit(TH1::kCanRebin);
-<
+      h_pullz_nTrk->GetXaxis()->SetTitle("Num of tracks");
->
+    h_gen = new PVHistograms();
->
+    h_gen->Init("generator");
->
+  }
->
+  if (analyze_) h_ana = new PVHistograms();
->
->
+  //Book histograms sensitive to data/mc
->
+  for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
->
+    string suffix;
->
+    edm::LogInfo("Debug")<<"datamode = "<< *it<<endl;
->
+    switch(*it) {
->
+    case 0: suffix = "";
->
+      break;
->
+    case 1: suffix = "_spl1_mct";
->
+      break;
->
+    case 2: suffix = "_spl2_mct";
->
+      break;
->
+    case 3: suffix = "_mct";
->
+      break;
-+
+    h_summary->Init("summary",suffix);
-+
-+
+    // Book residual, error and pull histograms for each nTrk bin
+    for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
+      stringstream ss;
+      ss << ntrk;
-<
+      string suffix = ss.str();
-<
+      h["resx_"+suffix]  = fs->make<TH1D> (TString("deltax_"+suffix),"x-distance (Rec - Sim)",100,-0.02,0.02);
-<
+      h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
-<
+      h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",100,-0.02,0.02);
-<
+      h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
-<
+      h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",100,-0.02,0.02);
-<
+      h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
-<
+      h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",100,-10.,10.);
-<
+      h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",100,-10.,10.);
-<
+      h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",100,-10.,10.);
-<
+      suffix.clear();
->
+      h_others->Init("others",suffix,ss.str());
-–
+  }
-–
+}
-+
+    // Book residual and pull histograms only when analyzeOntheFly is enabled
-+
+    if(analyze_) h_ana->InitA("analysis",suffix,"nTrk",nTrkMin_,nTrkMax_);
-+
+    suffix.clear();
-+
+  } // End of Book histograms sensitive to data/mc
-+
-+
-+
+}
+PVStudy::~PVStudy()
+  // do anything here that needs to be done at desctruction time
+  // (e.g. close files, deallocate resources etc.)
-+
+  theFile->cd();
-+
+  theFile->cd("Summary");
-+
+  h_pvtrk->Save();
-+
+  h_pixvtx->Save();
-+
+  h_misc->Save();
-+
+  h_summary->Save();
-+
+  if (!realData_)
-+
+    h_gen->Save();
-+
+  if (analyze_) {
-+
+    theFile->cd();
-+
+    theFile->cd("Results");
-+
+    h_ana->Save();
-+
+  }
-+
+  theFile->cd();
-+
+  theFile->cd("Others");
-+
+  h_others->Save();
-+
+  theFile->Close();
-<
->
+void PVStudy::setRootStyle() {
->
+  //
->
+  gROOT->SetStyle("Plain");
->
+  gStyle->SetFillColor(1);
->
+  gStyle->SetOptDate();
->
+  //   gStyle->SetOptStat(1111110);
->
+  //   gStyle->SetOptFit(0111);
->
+  //   gStyle->SetPadTickX(1);
->
+  //   gStyle->SetPadTickY(1);
->
+  gStyle->SetMarkerSize(0.5);
->
+  gStyle->SetMarkerStyle(8);
->
+  gStyle->SetGridStyle(3);
->
+  //gStyle->SetPadGridX(1);
->
+  //gStyle->SetPadGridY(1);
->
+  gStyle->SetPaperSize(TStyle::kA4);
->
+  gStyle->SetStatW(0.25);             // width of statistics box; default is 0.19
->
+  //   gStyle->SetStatH(0.10);             // height of statistics box; default is 0.1
->
+  gStyle->SetStatFormat("6.4g");      // leave default format for now
->
+  gStyle->SetTitleSize(0.055, "");    // size for pad title; default is 0.02
->
+  gStyle->SetLabelSize(0.03, "XYZ");  // size for axis labels; default is 0.04
->
+  gStyle->SetStatFontSize(0.08);      // size for stat. box
->
+  gStyle->SetTitleFont(32, "XYZ");    // times-bold-italic font (p. 153) for axes
->
+  gStyle->SetTitleFont(32, "");       // same for pad title
->
+  gStyle->SetLabelFont(32, "XYZ");    // same for axis labels
->
+  gStyle->SetStatFont(32);            // same for stat. box
->
+  gStyle->SetLabelOffset(0.006, "Y"); // default is 0.005
->
+  //
->
+  return;
->
+}
+//
+// member functions
+//
+std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
-<
+  std::vector<PVStudy::simPrimaryVertex> simpv;
->
+ std::vector<PVStudy::simPrimaryVertex> simpv;
+  const HepMC::GenEvent* evt=evtMC->GetEvent();
+  if (evt) {
-<
+    if(verbose_) {
-<
+      cout << "process id " << evt->signal_process_id()<<endl;
-<
+      cout << "signal process vertex " << ( evt->signal_process_vertex() ?
-<
+                                            evt->signal_process_vertex()->barcode() : 0 )   <<endl;
-<
+      cout << "number of vertices " << evt->vertices_size() << endl;
-<
+    }
->
+    edm::LogInfo("SimPVs") << "[getSimPVs] process id " << evt->signal_process_id()<<endl;
->
+    edm::LogInfo("SimPVs") << "[getSimPVs] signal process vertex " << ( evt->signal_process_vertex() ?
->
+                                                                        evt->signal_process_vertex()->barcode() : 0 ) <<endl;
->
+    edm::LogInfo("SimPVs") << "[getSimPVs] number of vertices " << evt->vertices_size() << endl;
-<
+    int idx=0;
->
+    int idx=0; int npv=0;
+    for(HepMC::GenEvent::vertex_const_iterator vitr= evt->vertices_begin();
-<
+        vitr != evt->vertices_end(); ++vitr )
-<
+      { // loop for vertex ...
-<
+        HepMC::FourVector pos = (*vitr)->position();
-<
+        //HepLorentzVector pos = (*vitr)->position();
-<
+        if (pos.t()>0) { continue;} // ?
-<
-<
+        bool hasMotherVertex=false;
-<
+        //      if(verbose_) cout << "mothers" << endl;
-<
+        for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
-<
+              mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
-<
+          HepMC::GenVertex * mv=(*mother)->production_vertex();
-<
+          if (mv) {hasMotherVertex=true;}
-<
+          //      if(verbose_) {
-<
+          //        cout << "\t";
-<
+          //        (*mother)->print();
-<
+          //      }
-<
+        }
-<
-<
+        if(hasMotherVertex) {continue;}
-<
-<
+        // could be a new vertex, check  all primaries found so far to avoid multiple entries
-<
+        const double mm2cm=0.1;
-<
+        simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
-<
+        simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
-<
+        for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
-<
+            v0!=simpv.end(); v0++){
-<
+          if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
-<
+            vp=&(*v0);
-<
+            break;
->
+        vitr != evt->vertices_end(); ++vitr ) { // loop for vertex ...
->
+      HepMC::FourVector pos = (*vitr)->position();
->
+      //HepLorentzVector pos = (*vitr)->position();
->
->
+      // t component of PV:
->
+      for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
->
+            mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
->
+        //        edm::LogInfo("SimPVs") << "Status = " << (*mother)->status() << endl;
->
+        HepMC::GenVertex * mv=(*mother)->production_vertex();
->
+        if( ((*mother)->status() == 3) && (!mv)) {
->
+          //        edm::LogInfo("SimPVs") << "npv= " << npv << endl;
->
+          if (npv == 0) {
->
+            h_gen->Fill1d("genPart_cT", pos.t()); // mm
->
+            h_gen->Fill1d("genPart_T", pos.t()/299.792458); // ns
-+
+          npv++;
-+
+      }
-+
+      //       if (pos.t()>0) { continue;} // for 22X when t of PV was not smeared
-+
-+
+      bool hasMotherVertex=false;
-+
+      if (verbose_) cout << "[getSimPVs] mothers of vertex[" << ++idx << "]: " << endl;
-+
+      for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
-+
+            mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
-+
+        HepMC::GenVertex * mv=(*mother)->production_vertex();
-+
+        // if (verbose_) cout << "Status = " << (*mother)->status() << endl;
-+
+        if (mv) {
-+
+          hasMotherVertex=true;
-+
+          if(!verbose_) break; //if verbose_, print all particles of gen vertices
-+
+        }
-+
+        if(verbose_) {
-+
+          cout << "\t";
-+
+          (*mother)->print();
-+
+        }
-+
+      }
-<
+        if(!vp){
-<
+          // this is a new vertex
-<
+          if(verbose_) cout << "this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
-<
+          simpv.push_back(sv);
-<
+          vp=&simpv.back();
-<
+        }else{
-<
+          if(verbose_) cout << "this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
-<
+        }
-<
+        vp->genVertex.push_back((*vitr)->barcode());
-<
+        // collect final state descendants
-<
+        for ( HepMC::GenVertex::particle_iterator
-<
+                daughter  = (*vitr)->particles_begin(HepMC::descendants);
-<
+              daughter != (*vitr)->particles_end(HepMC::descendants);
-<
+              ++daughter ) {
-<
+          if (isFinalstateParticle(*daughter)){
-<
+            if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
-<
+                 == vp->finalstateParticles.end()){
-<
+              vp->finalstateParticles.push_back((*daughter)->barcode());
-<
+              HepMC::FourVector m=(*daughter)->momentum();
-<
+              // the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
-<
+              // but adding FourVectors seems not to be foreseen
-<
+              vp->ptot.setPx(vp->ptot.px()+m.px());
-<
+              vp->ptot.setPy(vp->ptot.py()+m.py());
-<
+              vp->ptot.setPz(vp->ptot.pz()+m.pz());
-<
+              vp->ptot.setE(vp->ptot.e()+m.e());
-<
+              vp->ptsq+=(m.perp())*(m.perp());
-<
+              if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
-<
+                vp->nGenTrk++;
-<
+              }
->
+      if(hasMotherVertex) continue;
->
->
+      // could be a new vertex, check  all primaries found so far to avoid multiple entries
->
+      const double mm2cm=0.1;
->
+      simPrimaryVertex sv(pos.x()*mm2cm,pos.y()*mm2cm,pos.z()*mm2cm); // sim unit mm, rec unit cm
->
+      simPrimaryVertex *vp=NULL;  // will become non-NULL if a vertex is found and then point to it
->
+      for(vector<simPrimaryVertex>::iterator v0=simpv.begin();
->
+          v0!=simpv.end(); v0++){
->
+        if( (fabs(sv.x-v0->x)<1e-5) && (fabs(sv.y-v0->y)<1e-5) && (fabs(sv.z-v0->z)<1e-5)){
->
+          vp=&(*v0);
->
+          break;
->
+        }
->
+      }
->
->
+      if(!vp){
->
+        // this is a new vertex
->
+        edm::LogInfo("SimPVs") << "[getSimPVs] this is a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
->
+        simpv.push_back(sv);
->
+        vp=&simpv.back();
->
+      }else{
->
+        edm::LogInfo("SimPVs") << "[getSimPVs] this is not a new vertex " << sv.x << " " << sv.y << " " << sv.z << endl;
->
+      }
->
+      vp->genVertex.push_back((*vitr)->barcode());
->
+      // collect final state descendants
->
+      for ( HepMC::GenVertex::particle_iterator daughter  = (*vitr)->particles_begin(HepMC::descendants);
->
+            daughter != (*vitr)->particles_end(HepMC::descendants);
->
+            ++daughter ) {
->
+        if (isFinalstateParticle(*daughter)){
->
+          if ( find(vp->finalstateParticles.begin(), vp->finalstateParticles.end(),(*daughter)->barcode())
->
+               == vp->finalstateParticles.end()){
->
+            vp->finalstateParticles.push_back((*daughter)->barcode());
->
+            HepMC::FourVector m=(*daughter)->momentum();
->
+            // the next four lines used to be "vp->ptot+=m;" in the days of CLHEP::HepLorentzVector
->
+            // but adding FourVectors seems not to be foreseen
->
+            vp->ptot.setPx(vp->ptot.px()+m.px());
->
+            vp->ptot.setPy(vp->ptot.py()+m.py());
->
+            vp->ptot.setPz(vp->ptot.pz()+m.pz());
->
+            vp->ptot.setE(vp->ptot.e()+m.e());
->
+            vp->ptsq+=(m.perp())*(m.perp());
->
+            if ( (m.perp()>0.8) && (fabs(m.pseudoRapidity())<2.5) && isCharged( *daughter ) ){
->
+              vp->nGenTrk++;
-<
+        idx++;
-<
+      }
->
+      }//loop MC vertices daughters
->
+    }//loop MC vertices
+  return simpv;
+              ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
-<
+          primary =  ( gv->position().t()==0);
->
+          //primary =  ( gv->position().t()==0 );
->
+          primary = true;
->
+          h_gen->Fill1d("genPart_cT", gv->position().t()); // mm
->
+          h_gen->Fill1d("genPart_T", gv->position().t()/299.792458); // ns
->
+          //resonance= ( gp->mother() && isResonance(gp->mother()));  // in CLHEP/HepMC days
+          // no more mother pointer in the improved HepMC GenParticle
+          resonance= ( isResonance(*(gp->production_vertex()->particles_in_const_begin())));
+PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
+  using namespace edm;
-<
+  using reco::TrackCollection;
-<
-<
+  ftree_->SetBranchAddress("residual",&fres);
-<
+  ftree_->SetBranchAddress("error",&ferror);
-<
+  ftree_->SetBranchAddress("nTrk",&fntrk);
-<
-<
+  Handle<TrackCollection> tracks;
-<
+  iEvent.getByLabel(tracksLabel_,tracks);
->
+  using namespace reco;
->
->
+  edm::LogInfo("Debug")<<"[PVStudy]"<<endl;
->
+  //=======================================================
->
+  // Initialize Root-tuple variables if needed
->
+  //=======================================================
->
+  //if(saventuple_)
->
+    SetVarToZero();
->
->
+  //=======================================================
->
+  // Track accessors
->
+  //=======================================================
->
+  //trackCollection
->
+  static const reco::TrackCollection s_empty_trackColl;
->
+  const reco::TrackCollection *trackColl = &s_empty_trackColl;
->
+  edm::Handle<reco::TrackCollection> trackCollectionHandle;
->
+  iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle);
->
+  if( iEvent.getByLabel(trackCollectionTag_, trackCollectionHandle)) {
->
+    trackColl = trackCollectionHandle.product();
->
+  } else {
->
+    edm::LogInfo("Debug") << "[PVStudy] trackCollection cannot be found -> using empty collection of same type." <<endl;
->
+  }
->
+  //splitTrackCollection1
->
+  static const reco::TrackCollection s_empty_splitTrackColl1;
->
+  const reco::TrackCollection *splitTrackColl1 = &s_empty_splitTrackColl1;
->
+  edm::Handle<reco::TrackCollection> splitTrackCollection1Handle;
->
+  iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle);
->
+  if( iEvent.getByLabel(splitTrackCollection1Tag_, splitTrackCollection1Handle)) {
->
+    splitTrackColl1 = splitTrackCollection1Handle.product();
->
+  } else {
->
+    edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection1 cannot be found -> using empty collection of same type." <<endl;
->
+  }
->
+  //splitTrackCollection2
->
+  static const reco::TrackCollection s_empty_splitTrackColl2;
->
+  const reco::TrackCollection *splitTrackColl2 = &s_empty_splitTrackColl2;
->
+  edm::Handle<reco::TrackCollection> splitTrackCollection2Handle;
->
+  iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle);
->
+  if( iEvent.getByLabel(splitTrackCollection2Tag_, splitTrackCollection2Handle)) {
->
+    splitTrackColl2 = splitTrackCollection2Handle.product();
->
+  } else {
->
+   edm::LogInfo("Debug") << "[PVStudy] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
->
+  }
->
->
+  //=======================================================
->
+  // PVTX accessors
->
+  //=======================================================
->
+  //vertexCollection
->
+  static const reco::VertexCollection s_empty_vertexColl;
->
+  const reco::VertexCollection *vertexColl = &s_empty_vertexColl;
->
+  edm::Handle<reco::VertexCollection> vertexCollectionHandle;
->
+  iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle);
->
+  if( iEvent.getByLabel(vertexCollectionTag_, vertexCollectionHandle)) {
->
+    vertexColl = vertexCollectionHandle.product();
->
+  } else {
->
+   edm::LogInfo("Debug") << "[PVStudy] vertexCollection cannot be found -> using empty collection of same type." <<endl;
->
+  }
->
+  //splitVertexCollection1
->
+  static const reco::VertexCollection s_empty_splitVertexColl1;
->
+  const reco::VertexCollection *splitVertexColl1 = &s_empty_splitVertexColl1;
->
+  edm::Handle<reco::VertexCollection> splitVertexCollection1Handle;
->
+  iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle);
->
+  if( iEvent.getByLabel(splitVertexCollection1Tag_, splitVertexCollection1Handle)) {
->
+    splitVertexColl1 = splitVertexCollection1Handle.product();
->
+  } else {
->
+    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection1 cannot be found -> using empty collection of same type." <<endl;
->
+  }
->
+  //splitVertexCollection2
->
+  static const reco::VertexCollection s_empty_splitVertexColl2;
->
+  const reco::VertexCollection *splitVertexColl2 = &s_empty_splitVertexColl2;
->
+  edm::Handle<reco::VertexCollection> splitVertexCollection2Handle;
->
+  iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle);
->
+  if( iEvent.getByLabel(splitVertexCollection2Tag_, splitVertexCollection2Handle)) {
->
+    splitVertexColl2 = splitVertexCollection2Handle.product();
->
+  } else {
->
+    edm::LogInfo("Debug") << "[PVStudy] splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
->
+  }
-<
+  Handle<reco::VertexCollection> recVtxs;
-<
+  iEvent.getByLabel(vtxSample_, recVtxs);
->
->
+  //=======================================================
->
+  // GET pixelVertices
->
+  //=======================================================
->
+  static const reco::VertexCollection s_empty_pixelVertexColl;
->
+  const reco::VertexCollection *pixelVertexColl = &s_empty_pixelVertexColl;
->
+  edm::Handle<reco::VertexCollection>  pixelVertexCollectionHandle;
->
+  iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle);
->
+  if( iEvent.getByLabel(pixelVertexCollectionTag_, pixelVertexCollectionHandle)) {
->
+    pixelVertexColl = pixelVertexCollectionHandle.product();
->
+  } else {
->
+    edm::LogInfo("Debug") << "[PVStudy] pixelVertexCollection cannot be found. -> using empty collection of same type." <<endl;
->
+  }
-+
+  //=======================================================
-+
+  // BeamSpot accessors
-+
+  //=======================================================
-+
-+
+  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
-+
+  iEvent.getByLabel(bsSrc,recoBeamSpotHandle);
-+
+  reco::BeamSpot bs = *recoBeamSpotHandle;
-+
+  const Point beamSpot = recoBeamSpotHandle.isValid() ? Point(recoBeamSpotHandle->x0(), recoBeamSpotHandle->y0(), recoBeamSpotHandle->z0()) : Point(0, 0, 0);
-+
-+
+  edm::LogInfo("Debug")<<"[PVStudy] End accessing the track, beamSpot, primary vertex and pixelvertices collections"<<endl;
-+
-+
+  //=======================================================
-+
+  // MC simvtx accessor
-+
+  //=======================================================
+  if (!realData_) {
-<
+    Handle<SimVertexContainer> simVtxs;
->
+    edm::Handle<SimVertexContainer> simVtxs;
+    iEvent.getByLabel( simG4_, simVtxs);
-<
+    Handle<SimTrackContainer> simTrks;
->
+    edm::Handle<SimTrackContainer> simTrks;
+    iEvent.getByLabel( simG4_, simTrks);
-<
+  if(verbose_) printRecVtxs(recVtxs);
-<
->
+  //=======================================================
->
+  // GET PDT
->
+  //=======================================================
+  try{
+    iSetup.getData(pdt);
+  }catch(const Exception&){
-<
+    cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
->
+    edm::LogInfo("Debug") << "[PVStudy] Some problem occurred with the particle data table. This may not work !." <<endl;
-<
+  if (!realData_) {
-<
+    bool MC=false;
-<
+    Handle<HepMCProduct> evtMC;
-<
+    iEvent.getByLabel("generator",evtMC);
-<
+    if (!evtMC.isValid()) {
-<
+      MC=false;
-<
+      if(verbose_){
-<
+        cout << "no HepMCProduct found"<< endl;
-<
+      }
-<
+    } else {
-<
+      if(verbose_){
-<
+        cout << "generator HepMCProduct found"<< endl;
-<
+      }
-<
+      MC=true;
-<
+    }
-<
-<
+    if(MC){
-<
+      // make a list of primary vertices:
-<
+      std::vector<simPrimaryVertex> simpv;
-<
+      simpv=getSimPVs(evtMC,"");
-<
+      //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
-<
-<
+      int nsimtrk=0;
-<
+      int nrectrk=0;
-<
+      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
-<
+          vsim!=simpv.end(); vsim++){
-<
+        nsimtrk+=vsim->nGenTrk;
-<
+        // look for a matching reconstructed vertex
-<
+        vsim->recVtx=NULL;
-<
+        for(reco::VertexCollection::const_iterator vrec=recVtxs->begin();
-<
+            vrec!=recVtxs->end(); ++vrec){
-<
+          nrectrk=vrec->tracksSize();
-<
+          if(verbose_){
-<
+            cout << "sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
-<
+            cout << "Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
-<
+          }
-<
+          if ( matchVertex(*vsim,*vrec) ){
-<
+            // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
-<
+            if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
-<
+                   || (!vsim->recVtx) )
-<
+              {
-<
+                vsim->recVtx=&(*vrec);
-<
+              }
-<
+          }
-<
+        }
-<
+        if (vsim->recVtx){
-<
+          if(verbose_) {
-<
+            cout <<"primary matched " << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
-<
+          }
-<
+          fres[0] = vsim->recVtx->x()-vsim->x;
-<
+          fres[1] = vsim->recVtx->y()-vsim->y;
-<
+          fres[2] = vsim->recVtx->z()-vsim->z;
-<
+          ferror[0] = vsim->recVtx->xError();
-<
+          ferror[1] = vsim->recVtx->yError();
-<
+          ferror[2] = vsim->recVtx->zError();
-<
+          fntrk = nrectrk;
-<
+          ftree_->Fill();
-<
-<
+          h_deltax->Fill( vsim->recVtx->x()-vsim->x );
-<
+          h_deltay->Fill( vsim->recVtx->y()-vsim->y );
-<
+          h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
-<
+          h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
-<
+          h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
-<
+          h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
-<
+          pvinfo.push_back(PVStudy::PVInfo(res(vsim->recVtx->x()-vsim->x,
-<
+                                               vsim->recVtx->y()-vsim->y,
-<
+                                               vsim->recVtx->z()-vsim->z),
-<
+                                           error(vsim->recVtx->xError(),
-<
+                                                 vsim->recVtx->yError(),
-<
+                                                 vsim->recVtx->zError()),
-<
+                                           nrectrk));
-<
+          // Fill histo according to its track multiplicity
-<
+          fillHisto(res(vsim->recVtx->x()-vsim->x,
-<
+                        vsim->recVtx->y()-vsim->y,
-<
+                        vsim->recVtx->z()-vsim->z),
-<
+                    pull((vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
-<
+                         (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
-<
+                         (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError()),
-<
+                    nrectrk);
-<
+        }
-<
+        else{  // no rec vertex found for this simvertex
-<
+          if(verbose_) {
-<
+            cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
-<
+          }
->
+  //=======================================================
->
+  // Apply event cleaning for firstcoll data and MC
->
+  //=======================================================
->
+  if(applyEvtClean_) {
->
+    // =====Select on the trigger bits
->
+    edm::ESHandle<L1GtTriggerMenu> menuRcd;
->
+    iSetup.get<L1GtTriggerMenuRcd>().get(menuRcd) ;
->
+    const L1GtTriggerMenu* menu = menuRcd.product();
->
+    edm::Handle< L1GlobalTriggerReadoutRecord > gtRecord;
->
+    iEvent.getByLabel( edm::InputTag("gtDigis"), gtRecord);
->
->
+    bool isTechBit0 = false;
->
+    bool isTechBit40 = false;
->
+    bool isBeamHalo = false;
->
->
+    TechnicalTriggerWord tw = gtRecord->technicalTriggerWord();
->
+    if ( ! tw.empty() ) {
->
+      // loop over dec. bit to get total rate (no overlap)
->
+      for ( int i = 0; i < 64; ++i ) {
->
+        if ( tw[i] ) {
->
+          //cout<<"technical number "<<i<<"  "<<endl;
->
+          if (i == 0) isTechBit0  = true;
->
+          if (i < 40 && i > 35) isBeamHalo = true;   // The beamHalo bits are 36-39
->
+          if (i == 40 || i == 41) isTechBit40 = true;
-+
+    } // =====End select on the trigger bits
-+
+    // =====Loop over the trackColl to tet the fraction of HighPurity tracks
-+
+    int nTracks = 0;
-+
+    int nHighPurityTracks=0;
-+
+    double fHighPurity=0;
-+
-+
+    for (unsigned int i=0; i<trackCollectionHandle->size(); i++, nTracks++) {
-+
+      TrackRef tkref(trackCollectionHandle,i);
-+
+      if( (tkref->qualityMask() & 4 ) == 4) ++nHighPurityTracks;
-<
+  }// With MC
-<
-<
+  // Loop RecVtxs and find matched SimVtxs
-<
+  for(reco::VertexCollection::const_iterator v=recVtxs->begin();
-<
+      v!=recVtxs->end(); ++v){
-<
+    try {
-<
+      for(reco::Vertex::trackRef_iterator t = v->tracks_begin();
-<
+          t!=v->tracks_end(); t++) {
-<
+        // illegal charge
-<
+        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
-<
+          //      h_trkPtPV->Fill(0.);
-<
+        }
-<
+        else {
-<
+          h_trkPtPV->Fill((**t).pt());
-<
+        }
->
+    if(nTracks>0)
->
+      fHighPurity =  double(nHighPurityTracks)/double(nTracks);
->
+    if(verbose_)
->
+      cout<<"fraction of HighPurity track is "<<fHighPurity<<endl;
->
->
+    // Apply the event selection for MC events
->
+    if(realData_) {
->
+      if(!isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) {
->
+        glob_runno_ = iEvent.id().run();
->
+        glob_evtno_ = iEvent.id().event();
->
+        glob_ls_   = iEvent.luminosityBlock();
->
+        glob_bx_  = iEvent.bunchCrossing();
->
+        return;
->
+      }
->
+      else {
->
+        // TechBit 0, beamHalo are not simulted in MC
->
+        if ( !isTechBit0 || !isTechBit40 || isBeamHalo || (fHighPurity<0.2&&nTracks>10) || vertexColl->begin()->isFake()) return;
-<
+    catch (...) {
-<
+      // exception thrown when trying to use linked track
-<
+      //          h_trkPtPV->Fill(0.);
->
+  } // End of Apply event cleaning cuts for firstcoll data
->
->
->
->
+  //=======================================================
->
+  // Fill trackparameters of the input tracks to pvtx fitter
->
+  //=======================================================
->
+  edm::LogInfo("Debug")<<"[PVStudy] Start filling track parameters of the input tracks to pvtx fitter."<<endl;
->
+  //fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs)
->
+  // datatype: unsplittracks (0); splittracks1 (1);  splittracks2 (2);
->
+  fillTrackHisto(trackColl, 0, beamSpot);
->
+  fillTrackHisto(splitTrackColl1, 1, beamSpot);
->
+  fillTrackHisto(splitTrackColl2, 2, beamSpot);
->
+  edm::LogInfo("Debug")<<"[PVStudy] End filling track parameters of the input tracks to pvtx fitter."<<endl;
->
->
->
+  //=======================================================
->
+  // Fill pixelVertices related histograms
->
+  //=======================================================
->
+  nrecPV_pxlpvtx_ = int (pixelVertexColl->size());
->
+  if(pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake())) {
->
+    //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
->
+    fillTrackHistoInPV(pixelVertexColl, 4, false, true, beamSpot);
->
+    h_pixvtx->Fill1d("dzErr_pxlpvtx", pixelVertexColl->begin()->zError());
->
+    // Get the dZ error of the tracks in the leading pixelVertexColl
->
+    for(reco::Vertex::trackRef_iterator t = (pixelVertexColl->begin())->tracks_begin();
->
+        t!= (pixelVertexColl->begin())->tracks_end(); t++) {
->
->
+      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
->
+        //        h_pvtrk->Fill1d("trkPtPV", 0.);
->
+      }
->
+      else
->
+        h_pixvtx->Fill1d("trkdzErr_pxlpvtx", (**t).dzError());
-<
-<
+    h_nTrkPV->Fill(v->tracksSize());
-<
->
+    // Fill the track->dz() in the PV difference from first pixelpvtx
->
+    for(reco::Vertex::trackRef_iterator t = (vertexColl->begin())->tracks_begin();
->
+        t!= (vertexColl->begin())->tracks_end(); t++) {
->
+      // illegal charge
->
+      if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
->
+        //        h_pvtrk->Fill1d("trkPtPV", 0.);
->
+      }
->
+      else {
->
+        if(pixelVertexColl->size()>0) {
->
+          h_pixvtx->Fill1d("trkdz_pxlpvtxdz", (**t).dz() -  pixelVertexColl->begin()->z());
->
+          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_pxlpvtxdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/pixelVertexColl->begin()->zError());
->
+          h_pixvtx->Fill1d("trkdz_pxlpvtxdz_trkdzerr", fabs((**t).dz() -  pixelVertexColl->begin()->z())/(**t).dzError());
->
+          h_pixvtx->Fill1d("trkdzErr_pvtx", (**t).dzError());
->
+        }
->
+      }
->
+    }
->
+  }
->
->
+  //=======================================================
->
+  // Fill number of reconstructed vertices
->
+  //=======================================================
->
->
+  edm::LogInfo("Debug")<<"[PVStudy] Printing vertexCollection: "<<endl;
->
+  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection1: "<<endl;
->
+  edm::LogInfo("Debug")<<"[PVStudy] Printing splitVertexCollection2: "<<endl;
->
+  edm::LogInfo("Debug")<<"[PVStudy] Start filling pvtx parameters."<<endl;
->
+  if (verbose_) {
->
+    printRecVtxs(vertexCollectionHandle);
->
+    printRecVtxs(splitVertexCollection1Handle);
->
+    printRecVtxs(splitVertexCollection2Handle);
-+
+  nrecPV_ = int (vertexColl->size());
-+
+  nrecPV1_spl_ = int (splitVertexColl1->size());
-+
+  nrecPV2_spl_ = int (splitVertexColl2->size());
-+
-+
+  h_pvtrk->Fill1d("nrecPV", nrecPV_);
-+
+  h_pvtrk->Fill1d("nrecPV1_spl", nrecPV1_spl_);
-+
+  h_pvtrk->Fill1d("nrecPV2_spl", nrecPV2_spl_);
-+
+  h_misc->Fill1d("nrecPVDiff", double(nrecPV1_spl_)-double(nrecPV2_spl_));
-–
+  h_nTrk->Fill(tracks->size());
-<
+  for(TrackCollection::const_iterator itTrack = tracks->begin();
-<
+      itTrack != tracks->end();
-<
+      ++itTrack) {
-<
+    int charge = 0;
-<
+    charge = itTrack->charge();
-<
+    h_trkPt->Fill(itTrack->pt());
->
+ //=================================================================
->
+  // Fill track parameter ntuple/hist for tracks used in recoVertices
->
+  //=================================================================
->
->
+  //fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
->
+  if(vertexColl->size() > 0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()))
->
+    fillTrackHistoInPV(vertexColl, 0, true, true, beamSpot);
->
->
+  if(splitVertexColl1->size() > 0 && splitVertexColl1->begin()->isValid() && !(splitVertexColl1->begin()->isFake()))
->
+    fillTrackHistoInPV(splitVertexColl1, 1, false, true, beamSpot);
->
->
+  if(splitVertexColl2->size() > 0 && splitVertexColl2->begin()->isValid() && !(splitVertexColl2->begin()->isFake()))
->
+    fillTrackHistoInPV(splitVertexColl2, 2, false, true, beamSpot);
->
->
+  //=======================================================
->
+  // Compare offlinePrimaryVertices with pixelVertices
->
+  //=======================================================
->
+  if( (pixelVertexColl->size()>0 && pixelVertexColl->begin()->isValid() && !(pixelVertexColl->begin()->isFake()))
->
+      && (vertexColl->size()>0 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake())) ) {
->
+    h_pixvtx->Fill1d("nrecPV_minus_nrecPxlPV", double (nrecPV_ - nrecPV_pxlpvtx_));
->
+    // difference in reconstructed position of the leading pvtx
->
+    //edm::LogInfo("Debug")<<"recx_[0] = "<< recx_[0] << "recx_pxlpvtx_[0] = "<< recx_pxlpvtx_[0]<<endl;
->
+    //edm::LogInfo("Debug")<<"recy_[0] = "<< recy_[0] << "recy_pxlpvtx_[0] = "<< recy_pxlpvtx_[0]<<endl;
->
+    h_pixvtx->Fill1d("recxPV_minus_recxPxlPV", recx_[0] - recx_pxlpvtx_[0]);
->
+    h_pixvtx->Fill1d("recyPV_minus_recyPxlPV", recy_[0] - recy_pxlpvtx_[0]);
->
+    h_pixvtx->Fill1d("reczPV_minus_reczPxlPV", recz_[0] - recz_pxlpvtx_[0]);
-+
-+
+ //==========================================================
-+
+  // Fill secondary/primary min separations for multi-vertices
-+
+  //==========================================================
-+
-+
+  if(vertexColl->size()>1 && vertexColl->begin()->isValid() && !(vertexColl->begin()->isFake()) ) {
-+
-+
+    double min_xsep = 9999.0;
-+
+    double min_ysep = 9999.0;
-+
+    double min_zsep = 9999.0;
-+
+    double min_xsep_sign = 9999.0;
-+
+    double min_ysep_sign = 9999.0;
-+
+    double min_zsep_sign = 9999.0;
-+
+    double min_ntrksep = 9999.0;
-+
+    double min_sumpt2sep = 9999999.0;
-+
-+
+    edm::LogInfo("Debug")<<"[PVStudy] leading pvtx z = "<< vertexColl->begin()->z() <<endl;
-+
-+
+    // Looping through the secondary vertices to find the mininum separation to the primary
-+
+    for(reco::VertexCollection::const_iterator v=vertexColl->begin() + 1 ;
-+
+        v!=vertexColl->end(); ++v) {
-+
+      if(v->isValid() && ! v->isFake()) {
-+
+        // xsep
-+
+        if(fabs(v->x()- vertexColl->begin()->x())<min_xsep)
-+
+          min_xsep = fabs(v->x()- vertexColl->begin()->x());
-+
+        // ysep
-+
+        if(fabs(v->y()- vertexColl->begin()->y())<min_ysep)
-+
+          min_ysep = fabs(v->y()- vertexColl->begin()->y());
-+
+        // zsep
-+
+        if(fabs(v->z()- vertexColl->begin()->z())<min_zsep)
-+
+          min_zsep = fabs(v->z()- vertexColl->begin()->z());
-+
+        // xsep_sign
-+
+        double xsep_sign = fabs(v->x()- vertexColl->begin()->x())/max(v->xError(), vertexColl->begin()->xError());
-+
+        if(xsep_sign < min_xsep_sign)
-+
+          min_xsep_sign = xsep_sign;
-+
+        // ysep_sign
-+
+        double ysep_sign = fabs(v->y()- vertexColl->begin()->y())/max(v->yError(), vertexColl->begin()->yError());
-+
+        if(ysep_sign < min_ysep_sign)
-+
+          min_ysep_sign = ysep_sign;
-+
+        // zsep_sign
-+
+        double zsep_sign = fabs(v->z()- vertexColl->begin()->z())/max(v->zError(), vertexColl->begin()->zError());
-+
+        if(zsep_sign < min_zsep_sign)
-+
+          min_zsep_sign = zsep_sign;
-+
+        // ntrksep
-+
+        if( (double(vertexColl->begin()->tracksSize()) - double(v->tracksSize())) < min_ntrksep)
-+
+          min_ntrksep = double(vertexColl->begin()->tracksSize()) - double(v->tracksSize());
-+
+        // sumpt2sep
-+
+        if(fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin()))) < min_sumpt2sep)
-+
+          min_sumpt2sep = fabs(sumPtSquared(*v) - sumPtSquared(*(vertexColl->begin())));
-+
+      }
-+
+    }
-+
+    h_misc->Fill1d("min_xsep", min_xsep);
-+
+    h_misc->Fill1d("min_ysep", min_ysep);
-+
+    h_misc->Fill1d("min_zsep", min_zsep);
-+
+    h_misc->Fill1d("min_xsep_sign", min_xsep_sign);
-+
+    h_misc->Fill1d("min_ysep_sign", min_ysep_sign);
-+
+    h_misc->Fill1d("min_zsep_sign", min_zsep_sign);
-+
+    h_misc->Fill1d("min_ntrksep", min_ntrksep);
-+
+    h_misc->Fill1d("min_sumpt2sep", min_sumpt2sep);
-+
-+
+    min_zsep_ = min_zsep;
-+
+    min_ntrksep_ = min_ntrksep;
-+
+    min_sumpt2sep_ = min_sumpt2sep;
-+
-+
-+
+  } // End of  if(vertexColl->size()>1) {
-+
-+
+  edm::LogInfo("Debug")<<"[PVStudy] End filling pvtx parameters."<<endl;
-+
-+
-+
-+
+  //=======================================================
-+
+  //           PrimaryVertex Matching
-+
+  // 1. In z  |z1-z2|< zsigncut * max(sigmaz1, sigmaz2)
-+
+  // 2. |(nTrkPV1 - nTrkPV2)/(nTrkPV1+nTrkPV2)|<ntrkdiffcut_
-+
+  // Assume the first match is the primary vertex,
-+
+  // since vertexColl are sorted in decreasing order of sum(pT)
-+
+  //=======================================================
-+
-+
+  edm::LogInfo("Debug")<<"[PVStudy] matching pvtxs "<<endl;
-+
+  reco::VertexCollection s_empty_matchedVertexColl1;
-+
+  reco::VertexCollection *matchedVertexColl1 = &s_empty_matchedVertexColl1;
-+
+  matchedVertexColl1->reserve(splitVertexColl1->size());
-+
+  reco::VertexCollection s_empty_matchedVertexColl2;
-+
+  reco::VertexCollection *matchedVertexColl2 = &s_empty_matchedVertexColl2;
-+
+  matchedVertexColl2->reserve(splitVertexColl2->size());
-+
-+
+  bool stopmatching = false;
-+
+  for (size_t ivtx = 0; ivtx<splitVertexCollection1Handle->size() && !stopmatching; ++ivtx) {
-+
+    reco::VertexRef recvtx1(splitVertexCollection1Handle, ivtx);
-+
+    if( !(recvtx1->isValid()) || recvtx1->isFake()) break;
-+
+    for (size_t jvtx = ivtx; jvtx < splitVertexCollection2Handle->size(); ++jvtx) {
-+
+      //------------------------------------------------------------------------
-+
+      //== If only considering matching the first vertex of the splitVertexColl
-+
+      //------------------------------------------------------------------------
-+
+      if (ivtx!=0 || jvtx!=0) { stopmatching = true; break; }
-+
+      reco::VertexRef recvtx2(splitVertexCollection2Handle, jvtx);
-+
+      if( !(recvtx2->isValid()) || recvtx2->isFake()) break;
-+
+      edm::LogInfo("Debug")<<"[PVStudy] Matching splitVertexColl1: #  "<< ivtx<< " and splitVertexColl1: # "<<jvtx<<endl;
-+
+      edm::LogInfo("Debug")<<"[PVStudy] recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
-+
+      edm::LogInfo("Debug")<<"[PVStudy] recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
-+
-+
+      double twovtxsig = (recvtx1->z()-recvtx2->z())/max(recvtx1->zError(), recvtx2->zError());
-+
+      h_misc->Fill1d("twovtxzsign", twovtxsig);
-+
+      if(matchVertex(recvtx1, recvtx2, zsigncut_)) {
-+
+        edm::LogInfo("Debug")<<"[PVStudy] The two splitted vertices match in Z. "<<endl;
-+
-+
+        int nTrkPV1 = recvtx1->tracksSize();
-+
+        int nTrkPV2 = recvtx2->tracksSize();
-+
+        double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
-+
+        h_misc->Fill1d("nTrkPVDiff", nTrkPV1-nTrkPV2);
-+
+        h_misc->Fill1d("nTrkPVRelDiff", ntrkreldiff);
-+
+        if(fabs(ntrkreldiff)<ntrkdiffcut_) {
-+
+          edm::LogInfo("Debug")<<"[PVStudy] (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<endl;
-+
-+
+          matchedVertexColl1->push_back(*recvtx1);
-+
+          matchedVertexColl2->push_back(*recvtx2);
-+
+          stopmatching = true; // stop the matching when the first match is found!
-+
+          break;
-+
+        }
-+
+        else
-+
+          edm::LogInfo("Debug")<<"WARNING: (nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2) = " << ntrkreldiff<<", exceeding cut "<<ntrkdiffcut_<<endl;
-+
+      }
-+
+      else {
-+
+        edm::LogInfo("Debug")<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
-+
+        edm::LogInfo("Debug")<<"recvtx1->z() = " << recvtx1->z() << "+/- "<< recvtx1->zError() << "." << endl;
-+
+        edm::LogInfo("Debug")<<"recvtx2->z() = " << recvtx2->z() << "+/- "<< recvtx2->zError() << "." << endl;
-+
+      }
-+
+    }
-+
+  }
-+
+  edm::LogInfo("Debug")<<"[PVStudy] End matching pvtxs"<<endl;
-+
-+
+  //=======================================================
-+
+  //  Analyze matchedVertexColls
-+
+  //=======================================================
-+
+  edm::LogInfo("Debug")<<"[PVStudy] Begin analyzing the matchedVertexColl with size = "<< matchedVertexColl1->size()<< endl;
-+
-+
+  // Compare the reconstructed vertex position and calculate resolution/pulls
-+
+  if(matchedVertexColl1->size() != 0 && matchedVertexColl2->size() != 0) {
-+
+    //fillPvtxHisto(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple)
-+
+    fillTrackHistoInPV(matchedVertexColl1, 1, true, false, beamSpot);
-+
+    fillTrackHistoInPV(matchedVertexColl2, 2, true, false, beamSpot);
-+
-+
+    reco::VertexCollection::const_iterator v1;
-+
+    reco::VertexCollection::const_iterator v2;
-+
+    nrecPV_twovtx_ = 0;
-+
+    for(v1 = matchedVertexColl1->begin(), v2 = matchedVertexColl2->begin();
-+
+        v1!=matchedVertexColl1->end(), v2 != matchedVertexColl2->end();
-+
+        ++v1, ++v2) {
-+
+      h_misc->Fill1d("ndofPVDiff", v1->ndof() - v2->ndof());
-+
+      h_misc->Fill1d("ndofPVRelDiff", (v1->ndof()-v2->ndof())/(v1->ndof()+v2->ndof()));
-+
+      fres_[0] = (v1->x()-v2->x())/sqrt(2.0);
-+
+      fres_[1] = (v1->y()-v2->y())/sqrt(2.0);
-+
+      fres_[2] = (v1->z()-v2->z())/sqrt(2.0);
-+
+      ferror_[0] = sqrt(pow(v1->xError(),2)+pow(v2->xError(),2))/sqrt(2);
-+
+      ferror_[1] = sqrt(pow(v1->yError(),2)+pow(v2->yError(),2))/sqrt(2);
-+
+      ferror_[2] = sqrt(pow(v1->zError(),2)+pow(v2->zError(),2))/sqrt(2);
-+
+      fntrk_ =  (v1->tracksSize()+v2->tracksSize())/2;
-+
-+
+      h_summary->Fill1d("deltax", fres_[0] );
-+
+      h_summary->Fill1d("deltay", fres_[1] );
-+
+      h_summary->Fill1d("deltaz", fres_[2] );
-+
+      h_summary->Fill1d("pullx", fres_[0]/ferror_[0]);
-+
+      h_summary->Fill1d("pully", fres_[1]/ferror_[1]);
-+
+      h_summary->Fill1d("pullz", fres_[2]/ferror_[2]);
-+
+      h_summary->Fill1d("errPVx", ferror_[0] );
-+
+      h_summary->Fill1d("errPVy", ferror_[1] );
-+
+      h_summary->Fill1d("errPVz", ferror_[2] );
-+
+      pvinfo.push_back(PVStudy::PVInfo(res(fres_[0], fres_[1], fres_[2]),
-+
+                                       error(ferror_[0], ferror_[1], ferror_[2]),
-+
+                                       fntrk_));
-+
+      // Fill histo according to its track multiplicity, set datamode = 0 for pvtx using all tracks
-+
+      fillHisto(res(fres_[0], fres_[1], fres_[2]),
-+
+                error(ferror_[0], ferror_[1], ferror_[2]),
-+
+                fntrk_,0);
-+
+      // Fill the ntuple variables
-+
+      nTrkPV1_spl_twovtx_[nrecPV_twovtx_] = v1->tracksSize();
-+
+      nTrkPV2_spl_twovtx_[nrecPV_twovtx_] = v2->tracksSize();
-+
+      sumptsq1_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v1);
-+
+      sumptsq2_spl_twovtx_[nrecPV_twovtx_] = sumPtSquared(*v2);
-+
+      nTrkPV_twovtx_[nrecPV_twovtx_] = fntrk_;
-+
+      deltax_twovtx_[nrecPV_twovtx_] = fres_[0];
-+
+      deltay_twovtx_[nrecPV_twovtx_] = fres_[1];
-+
+      deltaz_twovtx_[nrecPV_twovtx_] = fres_[2];
-+
+      errx_twovtx_[nrecPV_twovtx_] = ferror_[0];
-+
+      erry_twovtx_[nrecPV_twovtx_] = ferror_[1];
-+
+      errz_twovtx_[nrecPV_twovtx_] = ferror_[2];
-+
+      pullx_twovtx_[nrecPV_twovtx_] = fres_[0]/ferror_[0];
-+
+      pully_twovtx_[nrecPV_twovtx_] = fres_[1]/ferror_[1];
-+
+      pullz_twovtx_[nrecPV_twovtx_] = fres_[2]/ferror_[2];
-+
+      nrecPV_twovtx_++;
-+
+    } // End of analyzing res/pull
-+
-+
+    //=======================================================
-+
+    // Fill simulated vertices
-+
+    //=======================================================
-+
+    if (!realData_) {
-+
+      bool MC=false;
-+
+      Handle<HepMCProduct> evtMC;
-+
+      iEvent.getByLabel("generator",evtMC);
-+
+      if (!evtMC.isValid()) {
-+
+        MC=false;
-+
+        edm::LogInfo("Debug") << "[PVStudy] no HepMCProduct found"<< endl;
-+
+      } else {
-+
+        edm::LogInfo("Debug") << "[PVStudy] generator HepMCProduct found"<< endl;
-+
+        MC=true;
-+
+      }
-+
-+
+      if(MC){
-+
+        // make a list of primary vertices:
-+
+        std::vector<simPrimaryVertex> simpv;
-+
+        simpv=getSimPVs(evtMC,"");
-+
+        //     simpv=getSimPVs(evtMC, simVtxs, simTrks);
-+
+        h_gen->Fill1d("nsimPV", simpv.size());
-+
+        nsimPV_ = simpv.size();
-+
-+
+        int isimpv = 0;
-+
+        for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
-+
+            vsim!=simpv.end(); vsim++, isimpv++){
-+
+          nsimTrkPV_[isimpv]  =vsim->nGenTrk;
-+
+          simx_[isimpv] = vsim->x;
-+
+          simy_[isimpv] = vsim->y;
-+
+          simz_[isimpv] = vsim->y;
-+
+          simptsq_[isimpv] = vsim->ptsq;
-+
-+
+          //fillMCHisto((std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype) {
-+
+          fillMCHisto(vsim, isimpv, splitVertexColl1, 1);
-+
+          fillMCHisto(vsim, isimpv, splitVertexColl2, 2);
-+
+          fillMCHisto(vsim, isimpv, vertexColl, 3);
-+
+        }
-+
+      } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
-+
+    } // End of   if (!realData_) {
-+
-+
+  } // End of  if(matchedVertexColl1->size() == 1 && matchedVertexColl2->size() == 1 ) {
-+
+  else
-+
+    edm::LogInfo("Debug")<<"[PVStudy] WARNING: Cannot find matching pvtxs"<<endl;
-+
-+
+  edm::LogInfo("Debug")<<"[PVStudy] End analyzing the matchedVertexColl"<<endl;
-+
-+
-+
-+
+  // Finally fill the ftree_
-+
+  if(saventuple_) {
-+
+    ftree_->Fill();
-+
+    pvtxtree_->Fill();
-+
+  }
-+
-+
+// ------------ method called once each job just before starting event loop  ------------
+void
+PVStudy::beginJob()
+// ------------ method called once each job just after ending the event loop  ------------
+void
+PVStudy::endJob() {
-<
+  double sqt2 = sqrt(2.);
-<
+  if (!realData_) {
-<
+    if (verbose_) cout << "Analyzing PV info" << endl;
-<
+    for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
-<
+      if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
-<
+      PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk);
-<
+      if (analyze_) {
-<
+//      if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
-<
+//      if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
-<
+//      if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
-<
+        if ( ipv.res_.x() > 0 ) h_resx_nTrk->Fill(ntrk,ipv.res_.x(), 1.);
-<
+        if ( ipv.res_.y() > 0 ) h_resy_nTrk->Fill(ntrk,ipv.res_.y(), 1.);
-<
+        if ( ipv.res_.z() > 0 ) h_resz_nTrk->Fill(ntrk,ipv.res_.z(), 1.);
-<
+        if ( ipv.pull_.x() > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pull_.x(), 1.);
-<
+        if ( ipv.pull_.y() > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pull_.y(), 1.);
-<
+        if ( ipv.pull_.z() > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pull_.z(), 1.);
->
+  edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
->
+  // Looping through the datamodes available
->
+  for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
->
+    string suffix;
->
+    edm::LogInfo("Analysis")<<"[endJob] datamode = "<< *it<<endl;
->
+    switch(*it) {
->
+    case 0: suffix = "";
->
+      break;
->
+    case 1: suffix = "_spl1_mct";
->
+      break;
->
+    case 2: suffix = "_spl2_mct";
->
+      break;
->
+    case 3: suffix = "_mct";
->
+      break;
->
+    }
->
+    suffix += "_nTrk";
->
+    if(analyze_) {
->
+      for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
->
+        edm::LogInfo("Analysis") << "[endJob] ntrk = " << ntrk << endl;
->
+        PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
->
+        if ( ipv.res_.x() > 0 ) h_ana->Fill2d("resx"+suffix, ntrk,ipv.res_.x());
->
+        if ( ipv.res_.y() > 0 ) h_ana->Fill2d("resy"+suffix, ntrk,ipv.res_.y());
->
+        if ( ipv.res_.z() > 0 ) h_ana->Fill2d("resz"+suffix, ntrk,ipv.res_.z());
->
+        if ( ipv.pull_.x() > 0 ) h_ana->Fill2d("pullx"+suffix, ntrk,ipv.pull_.x());
->
+        if ( ipv.pull_.y() > 0 ) h_ana->Fill2d("pully"+suffix, ntrk,ipv.pull_.y());
->
+        if ( ipv.pull_.z() > 0 ) h_ana->Fill2d("pullz"+suffix, ntrk,ipv.pull_.z());
-+
+    suffix.clear();
-+
+  }
-+
+  if(saventuple_) {
-+
+    file_->cd();
-+
+    ftree_->Write();
-+
+    pvtxtree_->Write();
-+
+    file_->Close();
-<
+  file_->cd();
-<
+  ftree_->Write();
-<
+  file_->Close();
->
->
-+
+// Match a recovertex with a simvertex
-+
+// ? Should the cut be parameter dependent?
-+
+// cut can be within n sigma of the vertex.zerror()
+bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim,
+                          const reco::Vertex & vrec)
-<
+  return (fabs(vsim.z-vrec.z())<0.0500); // =500um
->
+  if(vrec.isFake() || !vrec.isValid()) return false;
->
+  else
->
+    return (fabs(vsim.z-vrec.z())<0.0500); // =500um
-+
+// Match two reconstructed vertices
-+
+bool PVStudy::matchVertex( const reco::VertexRef & vrec1,
-+
+                           const reco::VertexRef & vrec2,
-+
+                           int zsigncut)
-+
+{
-+
+  double cut = zsigncut*max(vrec1->zError(), vrec2->zError());
-+
+  edm::LogInfo("Debug")<<"[matchVertex] The matching criteria in z is "<<cut<<endl;
-+
+  return (fabs(vrec1->z()-vrec2->z())<cut);
-+
+}
-+
-+
+bool PVStudy::isResonance(const HepMC::GenParticle * p){
+  double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
-<
+  if(verbose_) cout << "isResonance   " << p->pdg_id() << " " << ctau << endl;
->
+  edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
+  return  ctau >0 && ctau <1e-6;
-<
+void PVStudy::fillHisto(res r, pull p, int ntk){
->
+void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
+  stringstream ss;
+  ss << ntk;
-<
+  string suffix = ss.str();
->
+  string suffix;
->
+  if(datamode == 0 )  suffix = "_" + ss.str();
->
+  if(datamode == 1 )  suffix = "_spl1_mct_" + ss.str();
->
+  if(datamode == 2 )  suffix = "_spl2_mct_" + ss.str();
->
+  if(datamode == 3 )  suffix = "_mct_" + ss.str();
->
+  if (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
+  // Fill histograms of res and pull of ntk
-<
+  h["resx_"+suffix]->Fill(r.x());
-<
+  h["resy_"+suffix]->Fill(r.y());
-<
+  h["resz_"+suffix]->Fill(r.z());
-<
+  h["pullx_"+suffix]->Fill(p.x());
-<
+  h["pully_"+suffix]->Fill(p.y());
-<
+  h["pullz_"+suffix]->Fill(p.z());
->
+  h_others->Fill1d("deltax"+suffix, r.x());
->
+  h_others->Fill1d("deltay"+suffix, r.y());
->
+  h_others->Fill1d("deltaz"+suffix, r.z());
->
+  h_others->Fill1d("pullx"+suffix, r.x()/er.x());
->
+  h_others->Fill1d("pully"+suffix, r.y()/er.y());
->
+  h_others->Fill1d("pullz"+suffix, r.z()/er.z());
->
+  h_others->Fill1d("errPVx"+suffix, er.x());
->
+  h_others->Fill1d("errPVy"+suffix, er.y());
->
+  h_others->Fill1d("errPVz"+suffix, er.z());
-<
+PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
->
->
+PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
+  map<int,double> data;
+  data.clear();
-<
+  double fpar[2];
->
+  double fpar[2] = {-999,-999};
+  double cm2um = 10000;
+  stringstream ss;
+  ss << ntk;
-<
+  string suffix = ss.str();
-<
->
+  string suffix = ss.str();
->
+  if(datamode == 0 )  suffix = "_" + suffix;
->
+  if(datamode == 1 )  suffix = "_spl1_mct_" + suffix;
->
+  if(datamode == 2 )  suffix = "_spl2_mct_" + suffix;
->
+  if(datamode == 3 )  suffix = "_mct_" + suffix;
->
+  if(analyze_) {
+    for ( int i = 0; i < 6; ++i) {
-–
+      cout << "Here" << endl;
+      switch (i) {
+      case 0:
-<
+        fit(h["resx_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("deltax"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("deltax"+suffix), fpar);
+        data[i] = fpar[0]*cm2um;
+        break;
+      case 1:
-<
+        fit(h["resy_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("deltay"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("deltay"+suffix), fpar);
+        data[i] = fpar[0]*cm2um;
+        break;
+      case 2:
-<
+        fit(h["resz_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("deltaz"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("deltaz"+suffix), fpar);
+        data[i] = fpar[0]*cm2um;
+        break;
+      case 3:
-<
+        fit(h["pullx_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("pullx"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("pullx"+suffix), fpar);
+        data[i] = fpar[0];
+        break;
+      case 4:
-<
+        fit(h["pully_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("pully"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("pully"+suffix), fpar);
+        data[i] = fpar[0];
+        break;
+      case 5:
-<
+        fit(h["pullz_"+suffix], fpar);
->
+        if (!h_others->ReadHisto1D("pullz"+suffix)) break;
->
+        fit(h_others->ReadHisto1D("pullz"+suffix), fpar);
+        data[i] = fpar[0];
+        break;
-<
+      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
->
+      if (data[i] > 0) edm::LogInfo("Analysis") << "[Analysis] data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
+  else
+                            ntk);
-<
+void PVStudy::fit(TH1 *&hdist, double data[]){
->
->
+void PVStudy::fit(TH1 *hdist, double fitPars[]){
+  TAxis *axis0 = hdist->GetXaxis();
+  int nbin = axis0->GetLast();
+  double nOF = hdist->GetBinContent(nbin+1);
+  double nUF = hdist->GetBinContent(0);
-<
+  if ( verbose_ ){
-<
+    cout << "Last bin = " << nbin;
-<
+    cout << "; hdist: Overflow Entries = " << nOF;
-<
+    cout << "; Underflow Entries = " << nUF;
-<
+    cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
-<
+  }
-<
+  if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
-<
+    hdist->Fit("gaus","Q");
-<
+    TF1 *fgaus = hdist->GetFunction("gaus");
-<
+    if (verbose_) cout << "got function" << endl;
-<
+    data[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
->
+  //   double fitRange = axis0->GetBinUpEdge(nbin);
->
+  // double fitRange = axis0->GetXmax();
->
+  double fitRange = 2*hdist->GetRMS();
->
+  double sigMax[2] = {0.,0.};
->
->
+  edm::LogInfo("Analysis") << "[Fit] Last bin = " << nbin
->
+                           << "; hdist: Overflow Entries = " << nOF
->
+                           << "; Underflow Entries = " << nUF
->
+                           << "; hdist->GetEntries() = " << hdist->GetEntries()
->
+                           << "; fitting range = " << -fitRange << " to " << fitRange << endl;
->
->
+  if (hdist->GetEntries() - nOF - nUF >= 10) { // FIXME: for reasonable Gaussian fit
->
+    for (int bi = 0; bi < nbin; bi++) {
->
+      if ( (axis0->GetBinLowEdge(bi) < 0) && (hdist->GetBinContent(bi) > 0) ) {
->
+        sigMax[0] = axis0->GetBinLowEdge(bi);
->
+        if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
->
+      }
->
+      if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) {
->
+        sigMax[0] = axis0->GetBinUpEdge(bi);
->
+        if ( abs(sigMax[0]) > abs(sigMax[1]) ) sigMax[1] = abs(sigMax[0]);
->
+      }
->
+    }
->
+    //edm::LogInfo("Analysis") << "[Fit] Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;
->
->
+    TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
->
+    fgaus->SetParameter(1, 0.);
->
+    fgaus->SetParLimits(1, -fitRange/10., fitRange/10.);
->
+    fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
->
+    fgaus->SetLineColor(4);
->
+    hdist->Fit(fgaus,"QLRM");
->
+    gPad->Update();
->
+    TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
->
+    s->SetOptStat(1111111);
->
+    s->SetOptFit(0111);
->
+    gPad->Update();
->
+    //edm::LogInfo("Analysis") << "[Fit] got function" << endl;
->
+    fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
->
+  }
->
+  else
->
+    fitPars[0] = -999;
->
+}
->
->
->
+void PVStudy::fillTrackHisto(const reco::TrackCollection *trackColl, int datatype, const Point & bs) {
->
+  string suffix;
->
+  suffix = ""; // for unsplit tracks
->
+  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
->
+  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
->
->
+  h_pvtrk->Fill1d("nTrk"+suffix, trackColl->size());
->
+  for(reco::TrackCollection::const_iterator itTrack = trackColl->begin();
->
+      itTrack != trackColl->end();
->
+      ++itTrack) {
->
+    h_pvtrk->Fill1d("trkPt"+suffix, itTrack->pt());
->
+    h_pvtrk->Fill1d("trkDxy"+suffix, itTrack->dxy());
->
+    h_pvtrk->Fill1d("trkDxyCorr"+suffix, itTrack->dxy(bs));
->
+    h_pvtrk->Fill1d("trkDz"+suffix, itTrack->dz());
->
+    h_pvtrk->Fill1d("trkEta"+suffix, itTrack->eta());
->
+    h_pvtrk->Fill1d("trkPhi"+suffix, itTrack->phi());
->
+  }
->
+}
->
->
+void PVStudy::fillTrackHistoInPV(const reco::VertexCollection *vertexColl, int datatype, bool fillHisto, bool fillNtuple, const Point & bs) {
->
+  string suffix;
->
+  suffix = ""; // for unsplit tracks
->
+  if(datatype == 1) suffix = "1_spl"; // for splittracks 1
->
+  if(datatype == 2) suffix = "2_spl"; // for splittracks 2
->
+  int ivtx = 0;
->
+  for(reco::VertexCollection::const_iterator v=vertexColl->begin();
->
+      v!=vertexColl->end(); ++v, ++ivtx) {
->
+    if(v->isFake()) continue;
->
+    if(fillNtuple) {
->
+      if(datatype == 0) {
->
+        nTrkPV_[ivtx] = v->tracksSize();
->
+        sumptsq_[ivtx] = sumPtSquared(*v);
->
+        isValid_[ivtx] = v->isValid();
->
+        isFake_[ivtx] = v->isFake();
->
+        recx_[ivtx] = v->x();
->
+        recy_[ivtx] = v->y();
->
+        recz_[ivtx] = v->z();
->
+        recx_err_[ivtx] = v->xError();
->
+        recy_err_[ivtx] = v->yError();
->
+        recz_err_[ivtx] = v->zError();
->
+      }
->
+      if(datatype == 1) {
->
+        nTrkPV1_spl_[ivtx] = v->tracksSize();
->
+        sumptsq1_spl_[ivtx] = sumPtSquared(*v);
->
+        isValid1_spl_[ivtx] = v->isValid();
->
+        isFake1_spl_[ivtx] = v->isFake();
->
+        recx1_spl_[ivtx] = v->x();
->
+        recy1_spl_[ivtx] = v->y();
->
+        recz1_spl_[ivtx] = v->z();
->
+        recx1_err_spl_[ivtx] = v->xError();
->
+        recy1_err_spl_[ivtx] = v->yError();
->
+        recz1_err_spl_[ivtx] = v->zError();
->
->
+      }
->
+      if(datatype == 2) {
->
+        nTrkPV2_spl_[ivtx] = v->tracksSize();
->
+        sumptsq2_spl_[ivtx] = sumPtSquared(*v);
->
+        isValid2_spl_[ivtx] = v->isValid();
->
+        isFake2_spl_[ivtx] = v->isFake();
->
+        recx2_spl_[ivtx] = v->x();
->
+        recy2_spl_[ivtx] = v->y();
->
+        recz2_spl_[ivtx] = v->z();
->
+        recx2_err_spl_[ivtx] = v->xError();
->
+        recy2_err_spl_[ivtx] = v->yError();
->
+        recz2_err_spl_[ivtx] = v->zError();
->
+      }
->
+      if(datatype == 4) { // for pixelVertices
->
+        nTrkPV_pxlpvtx_[ivtx] = v->tracksSize();
->
+        sumptsq_pxlpvtx_[ivtx] = sumPtSquared(*v);
->
+        isValid_pxlpvtx_[ivtx] = v->isValid();
->
+        isFake_pxlpvtx_[ivtx] = v->isFake();
->
+        recx_pxlpvtx_[ivtx] = v->x();
->
+        recy_pxlpvtx_[ivtx] = v->y();
->
+        recz_pxlpvtx_[ivtx] = v->z();
->
+        recx_err_pxlpvtx_[ivtx] = v->xError();
->
+        recy_err_pxlpvtx_[ivtx] = v->yError();
->
+        recz_err_pxlpvtx_[ivtx] = v->zError();
->
+      }
->
+    }
->
+  }
->
+  // For histogram only fill the leading pvtx parameters
->
+  if(fillHisto) {
->
+    h_pvtrk->Fill1d("nTrkPV"+suffix, vertexColl->begin()->tracksSize());
->
+    h_pvtrk->Fill1d("ndofPV"+suffix, vertexColl->begin()->ndof());
->
+    int nHWTrkPV_ = 0;
->
->
+    try {
->
+      for(reco::Vertex::trackRef_iterator t = vertexColl->begin()->tracks_begin();
->
+          t!=vertexColl->begin()->tracks_end(); t++) {
->
+        // illegal charge
->
+        if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
->
+          //      h_pvtrk->Fill1d("trkPtPV", 0.);
->
+        }
->
+        else {
->
+          h_pvtrk->Fill1d("trkPtPV"+suffix, (**t).pt());
->
+          h_pvtrk->Fill1d("trkDxyPV"+suffix, (**t).dxy());
->
+          h_pvtrk->Fill1d("trkDxyCorrPV"+suffix, (**t).dxy(bs));
->
+          h_pvtrk->Fill1d("trkDzPV"+suffix, (**t).dz());
->
+          h_pvtrk->Fill1d("trkEtaPV"+suffix, (**t).eta());
->
+          h_pvtrk->Fill1d("trkPhiPV"+suffix, (**t).phi());
->
->
+          if(vertexColl->begin()->trackWeight(*t) > 0.5)
->
+            nHWTrkPV_++;
->
+        }
->
+      }
->
+    }
->
+    catch (...) {
->
+      // exception thrown when trying to use linked track
->
+      //          h_pvtrk->Fill1d("trkPtPV", 0.);
->
+    }
->
+    h_pvtrk->Fill1d("nHWTrkPV"+suffix, nHWTrkPV_);
->
+  }
->
->
->
+}
->
->
+void PVStudy::fillMCHisto(std::vector<simPrimaryVertex>::iterator vsim, int isimpv, const reco::VertexCollection *vtxColl, int datatype)
->
+{
->
+  std::string suffix;
->
+  // Set the vertexColl and histogram suffix according to datamode
->
+  if (datatype == 1) {
->
+    suffix = "_spl1_mct";
->
+    nrecPV_spl1_mct_ = 0;
->
+  }
->
+  if (datatype == 2) {
->
+    suffix = "_spl2_mct";
->
+    nrecPV_spl2_mct_ = 0;
->
+  }
->
+  if (datatype == 3) {
->
+    suffix = "_mct";
->
+    nrecPV_mct_ = 0;
-+
-+
+  //========================================================
-+
+  //  look for a matching reconstructed vertex in vertexColl
-+
+  //========================================================
-+
-+
+  for(reco::VertexCollection::const_iterator vrec=vtxColl->begin();
-+
+      vrec!=vtxColl->end(); ++vrec){
-+
+    int nrectrk = vrec->tracksSize();
-+
+    vsim->recVtx=NULL;
-+
-+
+    edm::LogInfo("Debug") << "[fillMCHisto] sim primary vertex x = " << vsim->x << "; y = " << vsim->y << "; z = " << vsim->z <<  endl;
-+
+    edm::LogInfo("Debug") << "[fillMCHisto] Is matched? " << (matchVertex(*vsim,*vrec)?"Yes":"No") << endl;
-+
-+
+    if ( matchVertex(*vsim,*vrec)  ) {
-+
+      vsim->recVtx=&(*vrec);
-+
-+
+      // if the matching critera are fulfilled, accept the rec-vertex that is closest in z
-+
+      //if(    ((vsim->recVtx) && (fabs(vsim->recVtx->position().z()-vsim->z)>fabs(vrec->z()-vsim->z)))
-+
+      //|| (!vsim->recVtx) )
-+
+      //vsim->recVtx=&(*vrec);
-+
+      //}
-+
+      edm::LogInfo("Debug") <<"[fillMCHisto] primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
-+
-+
+      double fres_mct[3];
-+
+      double ferror_mct[3];
-+
-+
+      fres_mct[0] = vsim->recVtx->x()-vsim->x;
-+
+      fres_mct[1] = vsim->recVtx->y()-vsim->y;
-+
+      fres_mct[2] = vsim->recVtx->z()-vsim->z;
-+
+      ferror_mct[0] = vsim->recVtx->xError();
-+
+      ferror_mct[1] = vsim->recVtx->yError();
-+
+      ferror_mct[2] = vsim->recVtx->zError();
-+
-+
+      h_summary->Fill1d("deltax"+suffix, fres_mct[0] );
-+
+      h_summary->Fill1d("deltay"+suffix, fres_mct[1] );
-+
+      h_summary->Fill1d("deltaz"+suffix, fres_mct[2] );
-+
+      h_summary->Fill1d("pullx"+suffix, fres_mct[0]/ferror_mct[0] );
-+
+      h_summary->Fill1d("pully"+suffix, fres_mct[1]/ferror_mct[1] );
-+
+      h_summary->Fill1d("pullz"+suffix, fres_mct[2]/ferror_mct[2] );
-+
+      h_summary->Fill1d("errPVx"+suffix, ferror_mct[0] );
-+
+      h_summary->Fill1d("errPVy"+suffix, ferror_mct[1] );
-+
+      h_summary->Fill1d("errPVz"+suffix, ferror_mct[2] );
-+
+      pvinfo.push_back(PVStudy::PVInfo(res(fres_mct[0], fres_mct[1], fres_mct[2]),
-+
+                                       error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
-+
+                                       nrectrk));
-+
+      // Fill histo according to its track multiplicity
-+
+      fillHisto(res(fres_mct[0], fres_mct[1], fres_mct[2]),
-+
+                error(ferror_mct[0], ferror_mct[1], ferror_mct[2]),
-+
+                nrectrk, datatype);
-+
-+
+      if(saventuple_) {
-+
+        //Fill the values for variables in ftree_
-+
+        if(datatype == 1) {
-+
+          nTrkPV_spl1_mct_[nrecPV_spl1_mct_] =   nrectrk;
-+
+          deltax_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
-+
+          deltay_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
-+
+          deltaz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[0];
-+
+          pullx_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[0]/ferror_mct[0];
-+
+          pully_spl1_mct_[nrecPV_spl1_mct_] = fres_mct[1]/ferror_mct[1];
-+
+          pullz_spl1_mct_[nrecPV_spl1_mct_] =  fres_mct[2]/ferror_mct[2];
-+
+          nrecPV_spl1_mct_++;
-+
+        }
-+
-+
+        if(datatype == 2) {
-+
+          nTrkPV_spl2_mct_[nrecPV_spl2_mct_] = nrectrk;
-+
+          deltax_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
-+
+          deltay_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
-+
+          deltaz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[0];
-+
+          pullx_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[0]/ferror_mct[0];
-+
+          pully_spl2_mct_[nrecPV_spl2_mct_] = fres_mct[1]/ferror_mct[1];
-+
+          pullz_spl2_mct_[nrecPV_spl2_mct_] =  fres_mct[2]/ferror_mct[2];
-+
+          nrecPV_spl2_mct_++;
-+
+        }
-+
+        if(datatype == 3) {
-+
+          nTrkPV_mct_[nrecPV_mct_] = nrectrk;
-+
+          deltax_mct_[nrecPV_mct_] =  fres_mct[0];
-+
+          deltay_mct_[nrecPV_mct_] =  fres_mct[0];
-+
+          deltaz_mct_[nrecPV_mct_] =  fres_mct[0];
-+
+          pullx_mct_[nrecPV_mct_] = fres_mct[0]/ferror_mct[0];
-+
+          pully_mct_[nrecPV_mct_] = fres_mct[1]/ferror_mct[1];
-+
+          pullz_mct_[nrecPV_mct_] =  fres_mct[2]/ferror_mct[2];
-+
+          nrecPV_mct_++;
-+
+        }
-+
+      } // End of  if(saventuple_) {
-+
+    } //  if ( matchVertex(*vsim,*vrec) ) {
-+
+    else {  // no rec vertex found for this simvertex
-+
+      edm::LogInfo("Debug") <<"[fillMCHisto] primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
-+
+    }
-+
+  }
-+
+}
-+
-+
+void PVStudy::SetVarToZero() {
-+
+  //Greg's variables
-+
+  fntrk_ = 0;
-+
+  //pvtx position (x,y,z) residual and error
-+
+  for(int i = 0; i<3;i++)    {
-+
+    fres_[i] = 0;
-+
+    ferror_[i] = 0;
-+
+  }
-+
-+
+  //Yanyan's variables
-+
+  // Number of reconstructed vertices
-+
+  nrecPV_ = 0;
-+
+  nrecPV1_spl_ = 0;
-+
+  nrecPV2_spl_ = 0;
-+
+  nrecPV_twovtx_ = 0;
-+
+  nsimPV_ = 0;
-+
-+
+  nrecPV_mct_ = 0;
-+
+  nrecPV_spl1_mct_ = 0;
-+
+  nrecPV_spl2_mct_ = 0;
-+
-+
+  // Mininum separation between the secondary pvtxes and leading pvtx
-+
+  min_zsep_ = 9999.0;
-+
+  min_ntrksep_ = 9999.0;
-+
+  min_sumpt2sep_ = -9999.0;
-+
-+
-+
+  for (int i = 0; i < nMaxPVs_; i++) {
-+
+    // recoVertices with all tracks
-+
+    nTrkPV_[i] = 0; // Number of tracks in the pvtx
-+
+    sumptsq_[i] = 0;
-+
+    isValid_[i] = -1;
-+
+    isFake_[i] = -1;
-+
+    recx_[i] = 0;
-+
+    recy_[i] = 0;
-+
+    recz_[i] = 0;
-+
+    recx_err_[i] = 0;
-+
+    recy_err_[i] = 0;
-+
+    recz_err_[i] = 0;
-+
-+
+    // recoVertices with splitTrack1
-+
+    nTrkPV1_spl_[i] = 0; // Number of tracks in the pvtx
-+
+    sumptsq1_spl_[i] = 0;
-+
+    isValid1_spl_[i] = -1;
-+
+    isFake1_spl_[i] = -1;
-+
+    recx1_spl_[i] = 0;
-+
+    recy1_spl_[i] = 0;
-+
+    recz1_spl_[i] = 0;
-+
+    recx1_err_spl_[i] = 0;
-+
+    recy1_err_spl_[i] = 0;
-+
+    recz1_err_spl_[i] = 0;
-+
-+
+    // recoVertices with splitTrack2
-+
+    nTrkPV2_spl_[i] = 0; // Number of tracks in the pvtx
-+
+    sumptsq2_spl_[i] = 0;
-+
+    isValid2_spl_[i] = -1;
-+
+    isFake2_spl_[i] = -1;
-+
+    recx2_spl_[i] = 0;
-+
+    recy2_spl_[i] = 0;
-+
+    recz2_spl_[i] = 0;
-+
+    recx2_err_spl_[i] = 0;
-+
+    recy2_err_spl_[i] = 0;
-+
+    recz2_err_spl_[i] = 0;
-+
-+
+    //pixelVertices
-+
+    nTrkPV_pxlpvtx_[i] = 0; // Number of tracks in the pvtx
-+
+    sumptsq_pxlpvtx_[i] = 0;
-+
+    isValid_pxlpvtx_[i] = -1;
-+
+    isFake_pxlpvtx_[i] = -1;
-+
+    recx_pxlpvtx_[i] = 0;
-+
+    recy_pxlpvtx_[i] = 0;
-+
+    recz_pxlpvtx_[i] = 0;
-+
+    recx_err_pxlpvtx_[i] = 0;
-+
+    recy_err_pxlpvtx_[i] = 0;
-+
+    recz_err_pxlpvtx_[i] = 0;
-+
-+
+    // matched two-vertices
-+
+    nTrkPV1_spl_twovtx_[i] = 0;
-+
+    nTrkPV2_spl_twovtx_[i] = 0;
-+
+    nTrkPV_twovtx_[i] = 0;
-+
+    sumptsq1_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
-+
+    sumptsq2_spl_twovtx_[i] = 0; // Sum of pT squred in the pvtx
-+
+    deltax_twovtx_[i] = 0;
-+
+    deltay_twovtx_[i] = 0;
-+
+    deltaz_twovtx_[i] = 0;
-+
+    errx_twovtx_[i] = 0;
-+
+    erry_twovtx_[i] = 0;
-+
+    errz_twovtx_[i] = 0;
-+
+    pullx_twovtx_[i] = 0;
-+
+    pully_twovtx_[i] = 0;
-+
+    pullz_twovtx_[i] = 0;
-+
-+
+    //simpv
-+
+    nsimTrkPV_[i] = 0;
-+
+    simx_[i] = 0;
-+
+    simy_[i] = 0;
-+
+    simz_[i] = 0;
-+
+    simptsq_[i] = 0;
-+
-+
+    // residual and pulls with mc truth required
-+
+    deltax_mct_[i] = 0;
-+
+    deltay_mct_[i] = 0;
-+
+    deltaz_mct_[i] = 0;
-+
+    pullx_mct_[i] = 0;
-+
+    pully_mct_[i] = 0;
-+
+    pullz_mct_[i] = 0;
-+
-+
+    deltax_spl1_mct_[i] = 0;
-+
+    deltay_spl1_mct_[i] = 0;
-+
+    deltaz_spl1_mct_[i] = 0;
-+
+    pullx_spl1_mct_[i] = 0;
-+
+    pully_spl1_mct_[i] = 0;
-+
+    pullz_spl1_mct_[i] = 0;
-+
-+
+    deltax_spl2_mct_[i] = 0;
-+
+    deltay_spl2_mct_[i] = 0;
-+
+    deltaz_spl2_mct_[i] = 0;
-+
+    pullx_spl2_mct_[i] = 0;
-+
+    pully_spl2_mct_[i] = 0;
-+
+    pullz_spl2_mct_[i] = 0;
-+
+  }
-+
+}
-+
-+
+double PVStudy::sumPtSquared(const reco::Vertex & v)  {
-+
+  double sum = 0.;
-+
+  double pT;
-+
+  for (reco::Vertex::trackRef_iterator it = v.tracks_begin(); it != v.tracks_end(); it++) {
-+
+    pT = (**it).pt();
-+
+    sum += pT*pT;
-+
+  }
-+
+  return sum;
-<
+//define this as a plug-in
-<
+DEFINE_FWK_MODULE(PVStudy);
->
->

Diff Legend

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

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents): Revision 1.3 by jengbou, Tue Sep 8 07:15:38 2009 UTC vs. Revision 1.17 by yygao, Tue Jan 26 12:06:49 2010 UTC

Diff Legend

Comparing UserCode/Jeng/PVStudy/plugins/PVStudy.cc (file contents):
Revision 1.3 by jengbou, Tue Sep 8 07:15:38 2009 UTC vs.
Revision 1.17 by yygao, Tue Jan 26 12:06:49 2010 UTC