PVStudy/plugins/PVStudy.cc

// -*- C++ -*-
//
// Package:    PVStudy
// Class:      PVStudy
// 
/**\class PVStudy PVStudy.cc UserCode/PVStudy/plugins/PVStudy.cc

Description: <one line class summary>

Implementation:
<Notes on implementation>
*/
//
// Original Author:  "Geng-yuan Jeng/UC Riverside"
//                    "Yanyan Gao/Fermilab ygao@fnal.gov"
//         Created:  Thu Aug 20 11:55:40 CDT 2009
// $Id: PVStudy.cc,v 1.6 2009/09/09 23:05:38 jengbou Exp $
//
//


// system include files
#include <memory>
#include <string>
#include <vector>
#include <iostream>
#include <sstream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
#include "UserCode/PVStudy/interface/PVStudy.h"
//
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"

// simulated vertices,..., add <use name=SimDataFormats/Vertex> and <../Track>
#include <SimDataFormats/Vertex/interface/SimVertex.h>
#include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
#include <SimDataFormats/Track/interface/SimTrack.h>
#include <SimDataFormats/Track/interface/SimTrackContainer.h>

//root
#include <TROOT.h>
#include <TF1.h>
#include <TString.h>
#include <TStyle.h>
#include <TPaveStats.h>
#include <TPad.h>

using namespace std;

PVStudy::PVStudy(const edm::ParameterSet& iConfig)
{
  //=======================================================================
  // 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");

  //now do what ever initialization is needed
  pvinfo.clear();


  // 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");
    ftree_->Branch("nrecPV",&nrecPV_,"nrecPV_/I");
    ftree_->Branch("nrecPV1_spl",&nrecPV1_spl_,"nrecPV1_spl_/I");
    ftree_->Branch("nrecPV2_spl",&nrecPV2_spl_,"nrecPV2_spl_/I");
  
    // If it is MC, get the resA, resB
    if(!realData_) {
      ftree_->Branch("residual_spl1_mct",&fres_spl1_mct_,"fres_spl1_mct_[3]/D");
      ftree_->Branch("error_spl1_mct",&ferror_spl1_mct_,"ferror_spl1_mct_[3]/D");
      ftree_->Branch("nTrkPV_spl1_mct",&fntrk_spl1_mct_,"fntrk_spl1_mct_/I");
      ftree_->Branch("residual_spl2_mct",&fres_spl2_mct_,"fres_spl2_mct_[3]/D");
      ftree_->Branch("error_spl2_mct",&ferror_spl2_mct_,"ferror_spl2_mct_[3]/D");
      ftree_->Branch("nTrkPV_spl2_mct",&fntrk_spl2_mct_,"fntrk_spl2_mct_/I");
      ftree_->Branch("residual_mct",&fres_mct_,"fres_mct_[3]/D");
      ftree_->Branch("error_mct",&ferror_mct_,"ferror_mct_[3]/D");
      ftree_->Branch("nTrkPV_mct",&fntrk_mct_,"fntrk_mct_/I");
    }
  }

  edm::Service<TFileService> fs;
  TFileDirectory subDir = fs->mkdir( "Summary" );
  TFileDirectory subDir1 = fs->mkdir( "Others" );
  //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );

  setRootStyle();
  
  //Book histograms for track and pvtx parameters for each splitted pvtx
  //i=0 : x (as in unsplit track collection)
  //i=1 : x1_spl_ 
  //i=2 : x2_spl_ 
  
  for(int i=0;i<3;i++)
    {  
      string suffix;
      if(i == 0) suffix = "";
      else {
        stringstream ss;
        ss << i;
        suffix =ss.str()+"_spl";  
      }
      h["nTrk"+suffix]   = subDir.make<TH1D>(TString("nTrk"+suffix), TString("Num of rec tracks"+suffix),300,0,300);
      h["trkPt"+suffix]  = subDir.make<TH1D>(TString("trkPt"+suffix), TString("Pt of rec tracks "+suffix),100,0,100);
      h["trkEta"+suffix] = subDir.make<TH1D>(TString("trkEta"+suffix), TString("#Eta of rec tracks "+suffix),100,-3,3);
      h["trkPhi"+suffix] = subDir.make<TH1D>(TString("trkPhi"+suffix), TString("#Phi of rec tracks "+suffix),100,-3.2,3.2);
      h["trkDxy"+suffix] = subDir.make<TH1D>(TString("trkDxy"+suffix), TString("Dxy of rec tracks "+suffix),100,-5,5);  
      h["trkDz"+suffix] = subDir.make<TH1D>(TString("trkDz"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50);  
      
      h["nTrkPV"+suffix]   = subDir.make<TH1D>(TString("nTrkPV"+suffix), TString("Num of rec tracks in PV"+suffix),300,0,300);
      h["trkPtPV"+suffix]  = subDir.make<TH1D>(TString("trkPtPV"+suffix), TString("Pt of rec tracks in "+suffix),100,0,100);
      h["trkEtaPV"+suffix] = subDir.make<TH1D>(TString("trkEtaPV"+suffix), TString("#Eta of rec tracks in PV"+suffix),100,-3,3);
      h["trkPhiPV"+suffix] = subDir.make<TH1D>(TString("trkPhiPV"+suffix), TString("#Phi of rec tracks in PV"+suffix),100,-3.2,3.2);
      h["trkDxyPV"+suffix] = subDir.make<TH1D>(TString("trkDxyPV"+suffix), TString("Dxy of rec tracks in PV"+suffix),100,-5,5);  
      h["trkDzPV"+suffix] = subDir.make<TH1D>(TString("trkDzPV"+suffix), TString("Dz of rec tracks "+suffix),100,-50,50); 
      h["nrecPV"+suffix]   = subDir.make<TH1D>(TString("nrecPV"+suffix), TString("Num of rec pvtx"+suffix),50,0,50);    
      suffix.clear();
    }
  h["nrecPVDiff"]     = subDir.make<TH1D>("nrecPVDiff","nrecPV1-nRecPV2",21,-10.5,10.5);
  h["nTrkPVDiff"]     = subDir.make<TH1D>("nTrkPVDiff","nTrkPV1-nTrkPV2",41,-20.5,20.5);
  h["nTrkPVRelDiff"]  = subDir.make<TH1D>("nTrkPVRelDiff","(nTrkPV1-nTrkPV2)/(nTrkPV1+nTrkPV2)",100,-1,1);
  
  //Book histograms sensitive to data/mc
  for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {
    string suffix; 
    if(verbose_)  cout<<"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["deltax"+suffix] = subDir.make<TH1D>(TString("deltax"+suffix), TString("x-residual pvtx"+suffix),800,-0.04,0.04);
    h["deltay"+suffix] = subDir.make<TH1D>(TString("deltay"+suffix), TString("y-residual pvtx"+suffix),800,-0.04,0.04);
    h["deltaz"+suffix] = subDir.make<TH1D>(TString("deltaz"+suffix), TString("z-residual pvtx"+suffix),800,-0.04,0.04);
    h["pullx"+suffix]  = subDir.make<TH1D>(TString("pullx"+suffix), TString("x-pull pvtx"+suffix),800,-10,10);
    h["pully"+suffix]  = subDir.make<TH1D>(TString("pully"+suffix), TString("y-pull pvtx"+suffix),800,-10,10);
    h["pullz"+suffix]  = subDir.make<TH1D>(TString("pullz"+suffix), TString("z-pull pvtx"+suffix),800,-10,10);
    h["errPVx"+suffix] = subDir.make<TH1D>(TString("errPVx"+suffix), TString("X"+suffix+" vertex error"),100,0.,0.02);
    h["errPVy"+suffix] = subDir.make<TH1D>(TString("errPVy"+suffix), TString("Y"+suffix+" vertex error"),100,0.,0.02);
    h["errPVz"+suffix] = subDir.make<TH1D>(TString("errPVz"+suffix), TString("Z"+suffix+" vertex error"),100,0.,0.02);
    
    // Book residual, error and pull histograms for each nTrk bin
    for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
      stringstream ss;
      ss << ntrk;
      string suffix2 = suffix+"_"+ss.str();
      h["deltax"+suffix2]  = subDir1.make<TH1D>(TString("deltax"+suffix2), TString("x-residual of pvtx"+suffix),100,-0.02,0.02);
      h["deltax"+suffix2]->GetXaxis()->SetTitle("cm");
      h["deltay"+suffix2]  = subDir1.make<TH1D>(TString("deltay"+suffix2), TString("y-residual of pvtx"+suffix),100,-0.02,0.02);
      h["deltay"+suffix2]->GetXaxis()->SetTitle("cm");
      h["deltaz"+suffix2]  = subDir1.make<TH1D>(TString("deltaz"+suffix2), TString("z-residual of pvtx"+suffix),100,-0.02,0.02);
      h["deltaz"+suffix2]->GetXaxis()->SetTitle("cm");
      h["pullx"+suffix2] = subDir1.make<TH1D>(TString("pullx"+suffix2), TString("x-pull of pvtx"+suffix),100,-10.,10.);
      h["pully"+suffix2] = subDir1.make<TH1D>(TString("pully"+suffix2), TString("y-pull of pvtx"+suffix),100,-10.,10.);
      h["pullz"+suffix2] = subDir1.make<TH1D>(TString("pullz"+suffix2), TString("z-pull of pvtx"+suffix),100,-10.,10.);
      h["errPVx"+suffix2] = subDir1.make<TH1D>(TString("errPVx"+suffix2), TString("X"+suffix+" vertex error"),100,0.,0.02);
      h["errPVy"+suffix2] = subDir1.make<TH1D>(TString("errPVy"+suffix2), TString("Y"+suffix+" vertex error"),100,0.,0.02);
      h["errPVz"+suffix2] = subDir1.make<TH1D>(TString("errPVz"+suffix2), TString("Z"+suffix+" vertex error"),100,0.,0.02);
      suffix2.clear();
    } // End of   for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
    
    // Book residual and pull histograms only when analyzeOntheFly is enabled
    if(analyze_) {
      h2["resx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resx"+suffix+"_nTrk"), TString("x-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
      h2["resx"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["resx"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["resx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resx"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["resy"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resy"+suffix+"_nTrk"), TString("y-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
      h2["resy"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["resy"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["resy"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resy"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["resz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("resz"+suffix+"_nTrk"), TString("z-resolution vs number of tracks in pvtx"+suffix),150,0,150,400,0.,200);
      h2["resz"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["resz"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["resz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resz"+suffix+"_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["pullx"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullx"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
      h2["pullx"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["pullx"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["pullx"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pullx"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["pully"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pully"+suffix+"_nTrk"), TString("y-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
      h2["pully"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["pully"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["pully"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pully"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["pullz"+suffix+"_nTrk"]= subDir.make<TH2D>(TString("pullz"+suffix+"_nTrk"), TString("x-pull vs number of tracks"+suffix),150,0,150,100,0.,2.);
      h2["pullz"+suffix+"_nTrk"]->SetMarkerStyle(21);
      h2["pullz"+suffix+"_nTrk"]->SetMarkerColor(4);
      h2["pullz"+suffix+"_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pullz"+suffix+"_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
    } // End of   if(analyze_) {
    suffix.clear();
  } // End of Book histograms sensitive to data/mc  
  // Book MC only plots
  if (!realData_) {  
    h["genPart_T"]      = subDir.make<TH1D>("genPart_T","t component of gen particles",300,-0.5,0.5);
    h["genPart_T"]->GetXaxis()->SetTitle("t (nanosecond)");
    h["genPart_cT"]     = subDir.make<TH1D>("genPart_cT","ct component of gen particles",300,-150.,150.);
    h["genPart_cT"]->GetXaxis()->SetTitle("ct (mm)");    
    h["nsimPV"]         = subDir.make<TH1D>("nsimPV","Num of sim PV",51,-0.5,50.5);  
  }
}

PVStudy::~PVStudy()
{
 
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)

}

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;
  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;
    }

    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();

      // t component of PV:
      for ( HepMC::GenVertex::particle_iterator mother  = (*vitr)->particles_begin(HepMC::parents);
            mother != (*vitr)->particles_end(HepMC::parents); ++mother ) {
        //        std::cout << "Status = " << (*mother)->status() << std::endl;
        HepMC::GenVertex * mv=(*mother)->production_vertex();
        if( ((*mother)->status() == 3) && (!mv)) {
          //        std::cout << "npv= " << npv << std::endl;
          if (npv == 0) {
            h["genPart_cT"]->Fill(pos.t()); // mm
            h["genPart_T"]->Fill(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 << "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();
        //        std::cout << "Status = " << (*mother)->status() << std::endl;
        if (mv) {
          hasMotherVertex=true;
          if(!verbose_) break; //if verbose_, print all particles of gen vertices
        }
        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;
        }
      }

      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++;
            }
          }
        }
      }
    }
  }
  return simpv;
}

std::vector<PVStudy::simPrimaryVertex> PVStudy::getSimPVs(const Handle<HepMCProduct> evtMC, 
                                                          const Handle<SimVertexContainer> simVtxs, 
                                                          const Handle<SimTrackContainer> simTrks)
{
  // simvertices don't have enough information to decide, 
  // genVertices don't have the simulated coordinates  ( with VtxSmeared they might)
  // go through simtracks to get the link between simulated and generated vertices
  std::vector<PVStudy::simPrimaryVertex> simpv;
  int idx=0;
  for(SimTrackContainer::const_iterator t=simTrks->begin();
      t!=simTrks->end(); ++t){
    if ( !(t->noVertex()) && !(t->type()==-99) ){
      double ptsq=0;
      bool primary=false;   // something coming directly from the primary vertex
      bool resonance=false; // resonance
      bool track=false;     // undecayed, charged particle
      HepMC::GenParticle* gp=evtMC->GetEvent()->barcode_to_particle( (*t).genpartIndex() );
      if (gp) {
        HepMC::GenVertex * gv=gp->production_vertex();
        if (gv) {
          for ( HepMC::GenVertex::particle_iterator 
                  daughter =  gv->particles_begin(HepMC::descendants);
                daughter != gv->particles_end(HepMC::descendants);
                ++daughter ) {
            if (isFinalstateParticle(*daughter)){
              ptsq+=(*daughter)->momentum().perp()*(*daughter)->momentum().perp();
            }
          }
          primary =  ( gv->position().t()==0);
          //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())));
          if (gp->status()==1){
            //track=((pdt->particle(gp->pdg_id()))->charge() != 0);
            track=not isCharged(gp);
          }
        }
      }

      const HepMC::FourVector & v=(*simVtxs)[t->vertIndex()].position();
      //const HepLorentzVector & v=(*simVtxs)[t->vertIndex()].position();
      if(primary or resonance){
        {
          // check all primaries found so far to avoid multiple entries
          bool newVertex=true;
          for(std::vector<PVStudy::simPrimaryVertex>::iterator v0=simpv.begin();
              v0!=simpv.end(); v0++){
            if( (fabs(v0->x-v.x())<0.001) && (fabs(v0->y-v.y())<0.001) && (fabs(v0->z-v.z())<0.001) ){
              if (track) {
                v0->simTrackIndex.push_back(idx);
                if (ptsq>(*v0).ptsq){(*v0).ptsq=ptsq;}
              }
              newVertex=false;
            }
          }
          if(newVertex && !resonance){
            PVStudy::simPrimaryVertex anotherVertex(v.x(),v.y(),v.z());
            if (track) anotherVertex.simTrackIndex.push_back(idx);
            anotherVertex.ptsq=ptsq;
            simpv.push_back(anotherVertex);
          }
        }// 
      }
       
    }// simtrack has vertex and valid type
    idx++;
  }//simTrack loop
  return simpv;
}

// ------------ method called to for each event  ------------
void
PVStudy::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;
  using namespace reco;
  
  if (verbose_)
    cout<<"PVStudy::analyze():"<<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 {
    cout << "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 {
    cout << "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 {
    cout << "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 {
    cout << "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 {
    cout << "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 {
    cout << "splitVertexCollection2 cannot be found -> using empty collection of same type." <<endl;
  }

  if(verbose_) cout<<"Done accessing the track and vertex collections"<<endl;
  
  //=======================================================
  // MC simvtx accessor
  //=======================================================
  if (!realData_) {
    edm::Handle<SimVertexContainer> simVtxs;
    iEvent.getByLabel( simG4_, simVtxs);
    
    edm::Handle<SimTrackContainer> simTrks;
    iEvent.getByLabel( simG4_, simTrks);
  }

  //=======================================================
  // GET PDT
  //=======================================================
  try{
    iSetup.getData(pdt);
  }catch(const Exception&){
    cout << "Some problem occurred with the particle data table. This may not work !." <<endl;
  }

  //=======================================================
  // Fill trackparameters of the input tracks to pvtx fitter
  //=======================================================

  // GeneralTracks  
  h["nTrk"]->Fill(trackColl->size());
  for(TrackCollection::const_iterator itTrack = trackColl->begin();
      itTrack != trackColl->end();                      
      ++itTrack) {
    h["trkPt"]->Fill(itTrack->pt());
    h["trkDxy"]->Fill(itTrack->dxy());
    h["trkDz"]->Fill(itTrack->dz());
    h["trkEta"]->Fill(itTrack->eta());
    h["trkPhi"]->Fill(itTrack->phi());
  }
  
  //SplittedTracks1 
  h["nTrk1_spl"]->Fill(splitTrackColl1->size());
  for(TrackCollection::const_iterator itTrack = splitTrackColl1->begin();
      itTrack != splitTrackColl1->end();                      
      ++itTrack) {
    h["trkPt1_spl"]->Fill(itTrack->pt());
    h["trkDxy1_spl"]->Fill(itTrack->dxy());
    h["trkDz1_spl"]->Fill(itTrack->dz());
    h["trkEta1_spl"]->Fill(itTrack->eta());
    h["trkPhi1_spl"]->Fill(itTrack->phi());
  }
  //SplittedTracks2  
  h["nTrk2_spl"]->Fill(splitTrackColl2->size());
  for(TrackCollection::const_iterator itTrack = splitTrackColl2->begin();
      itTrack != splitTrackColl2->end();                      
      ++itTrack) {
    h["trkPt2_spl"]->Fill(itTrack->pt());
    h["trkDxy2_spl"]->Fill(itTrack->dxy());
    h["trkDz2_spl"]->Fill(itTrack->dz());
    h["trkEta2_spl"]->Fill(itTrack->eta());
    h["trkPhi2_spl"]->Fill(itTrack->phi());
  }
  
  if(verbose_) 
    cout<<"Done filling track parameters of the input tracks to pvtx fitter."<<endl;
 
  //=======================================================
  // Fill reconstructed vertices 
  //=======================================================
  if(verbose_)  {
    cout<<"PVStudy::analyze() Printing vertexCollection: "<<endl;
    printRecVtxs(vertexCollectionHandle);
    cout<<"PVStudy::analyze() Printing splitVertexCollection1: "<<endl;
    printRecVtxs(splitVertexCollection1Handle);
    cout<<"PVStudy::analyze() Printing splitVertexCollection2: "<<endl;
    printRecVtxs(splitVertexCollection2Handle);
  }
  
  h["nrecPV"]->Fill(vertexColl->size());
  h["nrecPV1_spl"]->Fill(splitVertexColl1->size());
  h["nrecPV2_spl"]->Fill(splitVertexColl2->size());
  h["nrecPVDiff"]->Fill(double(splitVertexColl1->size())-double(splitVertexColl2->size())); 
  //Fill the tree variables
  nrecPV_ = vertexColl->size();
  nrecPV1_spl_ = splitVertexColl1->size();
  nrecPV2_spl_ = splitVertexColl2->size();
  
  // Fill the track paramters for the vertexColl
  for(reco::VertexCollection::const_iterator v=vertexColl->begin(); 
      v!=vertexColl->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());
          h["trkDxyPV"]->Fill((**t).dxy());
          h["trkDzPV"]->Fill((**t).dz());
          h["trkEtaPV"]->Fill((**t).eta());
          h["trkPhiPV"]->Fill((**t).phi());
        }
      }
    }
    catch (...) {
      // exception thrown when trying to use linked track
      //          h["trkPtPV"]->Fill(0.);
    }
  }
  
  
  // == ignore multi-vertices (for now)
  //take single pvt collection and compare
  if(splitVertexColl1->size() == 1 && splitVertexColl2->size() == 1 ) {
    const reco::Vertex recvtx1 = *(splitVertexColl1->begin());
    const reco::Vertex recvtx2 = *(splitVertexColl2->begin());
    if (matchVertex(recvtx1, recvtx2))  {
      int nTrkPV1 = recvtx1.tracksSize();
      int nTrkPV2 = recvtx2.tracksSize();     
      h["nTrkPV1_spl"]->Fill(nTrkPV1);
      h["nTrkPV2_spl"]->Fill(nTrkPV2);
      double ntrkreldiff = double(nTrkPV1-nTrkPV2)/double(nTrkPV1+nTrkPV2);
      h["nTrkPVDiff"]->Fill(nTrkPV1-nTrkPV2);
      h["nTrkPVRelDiff"]->Fill(ntrkreldiff);
      
      // Fill the track paramters for the splitVertexColl1
      try {
        for(reco::Vertex::trackRef_iterator t = recvtx1.tracks_begin(); 
            t!=recvtx1.tracks_end(); t++) {
          // illegal charge
          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
            //    h["trkPtPV"]->Fill(0.);
          }
          else {
            h["trkPtPV1_spl"]->Fill((**t).pt());
            h["trkDxyPV1_spl"]->Fill((**t).dxy());
            h["trkDzPV1_spl"]->Fill((**t).dz());
            h["trkEtaPV1_spl"]->Fill((**t).eta());
            h["trkPhiPV1_spl"]->Fill((**t).phi());
          }
        }
      }
      catch (...) {
        // exception thrown when trying to use linked track
        //        h["trkPtPV"]->Fill(0.);
      }
      
      // Fill the track paramters for the splitVertexColl2
      try {
        for(reco::Vertex::trackRef_iterator t = recvtx2.tracks_begin(); 
            t!=recvtx2.tracks_end(); t++) {
          // illegal charge
          if ( (**t).charge() < -1 || (**t).charge() > 1 ) {
            //    h["trkPtPV"]->Fill(0.);
          }
          else {
            h["trkPtPV2_spl"]->Fill((**t).pt());
            h["trkDxyPV2_spl"]->Fill((**t).dxy());
            h["trkDzPV2_spl"]->Fill((**t).dz());
            h["trkEtaPV2_spl"]->Fill((**t).eta());
            h["trkPhiPV2_spl"]->Fill((**t).phi());
          }
        }
      }
      catch (...) {
        // exception thrown when trying to use linked track
        //        h["trkPtPV"]->Fill(0.);
      }
      
      if(verbose_) 
        cout<<"Done filling track parameters in reconstruced vertices."<<endl;

      if(fabs(ntrkreldiff)<ntrkdiffcut_) {
        fres_[0] = (recvtx1.x()-recvtx2.x())/sqrt(2.0);
        fres_[1] = (recvtx1.y()-recvtx2.y())/sqrt(2.0);
        fres_[2] = (recvtx1.z()-recvtx2.z())/sqrt(2.0);       
        ferror_[0] = sqrt(pow(recvtx1.xError(),2)+pow(recvtx2.xError(),2))/sqrt(2);
        ferror_[1] = sqrt(pow(recvtx1.yError(),2)+pow(recvtx2.yError(),2))/sqrt(2);
        ferror_[2] = sqrt(pow(recvtx1.zError(),2)+pow(recvtx2.zError(),2))/sqrt(2);
        fntrk_ =  (nTrkPV1+nTrkPV2)/2;
        
        h["deltax"]->Fill( fres_[0] );
        h["deltay"]->Fill( fres_[1] );
        h["deltaz"]->Fill( fres_[2] );
        h["pullx"]->Fill( fres_[0]/ferror_[0]);
        h["pully"]->Fill( fres_[1]/ferror_[1]);
        h["pullz"]->Fill( fres_[2]/ferror_[2]);
        h["errPVx"]->Fill( ferror_[0] );
        h["errPVy"]->Fill( ferror_[1] );
        h["errPVz"]->Fill( 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 
        fillHisto(res(fres_[0], fres_[1], fres_[2]),
                  error(ferror_[0], ferror_[1], ferror_[2]),
                  fntrk_,0);  

        ftree_->Fill();  
      } // End of  if(fabs(ntrkreldiff)<ntrkdiffcut_) {
    } // End of  if (matchVertex(recvtx1, recvtx2))  {
    else   {
      cout<<"WARNING: The two reconstructed vertices do not match in z: "<<endl;
      cout<<"recvtx1.z() = " << recvtx1.z() <<"; recvtx2.z() = " <<  recvtx2.z() <<endl;
    }
  } // End of fill histograms for two-vertices Method

  if(verbose_) 
    cout<<"Done filling res/pull with two-vertices method"<<endl;

  //=======================================================
  // Fill simulated vertices 
  //=======================================================
  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);
      h["nsimPV"]->Fill(simpv.size());
      int nsimtrk=0;
      int nrectrk=0;
      
      for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin();
          vsim!=simpv.end(); vsim++){
        nsimtrk+=vsim->nGenTrk;
        // Loop over the datamodes, for MC, fill only 1, 2, 3
        for (int i=1;i<=3;i++) {
          static const reco::VertexCollection s_empty_vtxColl;
          const reco::VertexCollection *vtxColl = &s_empty_vtxColl;
          std::string suffix;
          // Set the vertexColl and histogram suffix according to datamode
          if (i == 1) {
            vtxColl = splitVertexColl1; 
            suffix = "_spl1_mct";
          }
          if (i == 2) {
            vtxColl =  splitVertexColl2; 
            suffix = "_spl2_mct";
          }
          if (i == 3) {
            vtxColl =  vertexColl; 
            suffix = "_mct";
          }
          //========================================================
          //  look for a matching reconstructed vertex in vertexColl
          //========================================================    
          vsim->recVtx=NULL;
          for(reco::VertexCollection::const_iterator vrec=vtxColl->begin(); 
              vrec!=vtxColl->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);
            }// End of finding the matched reco_vertex
            
            if (vsim->recVtx){
              if(verbose_) {
                cout <<"primary matched in vertexColl" << vsim->x << " " << vsim->y << " " << vsim->z <<  endl;
              } 
              double fres_mct[3];
              double ferror_mct[3];
              int fntrk_mct; 
              
              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();
              fntrk_mct = nrectrk;
              
              //Fill the values for variables in ftree_
              if(i == 1) {
                for(int j = 0;j<3;j++)
                  fres_spl1_mct_[j] = fres_mct[j];
                fntrk_spl1_mct_ =  fntrk_mct;
              }
              if(i == 2) {
                for(int j = 0;j<3;j++)
                  fres_spl2_mct_[j] = fres_mct[j];
                fntrk_spl2_mct_ =  fntrk_mct;
              }
              if(i == 3) {
                for(int j=0;j<3;j++)
                  fres_mct_[j] = fres_mct[j];
                fntrk_mct_ =  fntrk_mct;
              }
              
              h["deltax"+suffix]->Fill( fres_mct[0] );
              h["deltay"+suffix]->Fill( fres_mct[1] );
              h["deltaz"+suffix]->Fill( fres_mct[2] );
              h["pullx"+suffix]->Fill( fres_mct[0]/ferror_mct[0] );
              h["pully"+suffix]->Fill( fres_mct[1]/ferror_mct[1] );
              h["pullz"+suffix]->Fill( fres_mct[2]/ferror_mct[2] );
              h["errPVx"+suffix]->Fill( ferror_mct[0] );
              h["errPVy"+suffix]->Fill( ferror_mct[1] );
              h["errPVz"+suffix]->Fill( 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]),
                                               fntrk_mct));
              // 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]),
                        fntrk_mct, i);
              ftree_->Fill();
              suffix.clear();
            }
            else {  // no rec vertex found for this simvertex
              if(verbose_) 
                cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
            }
          }
        } // === End of Looping through vertexColl
      } 
    } // End of  for(std::vector<simPrimaryVertex>::iterator vsim=simpv.begin()
  } // End of Filling MC variables
}

// ------------ 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() {
  if (verbose_) cout << "Analyzing PV info" << endl;
  // Looping through the datamodes available
  for (vector<int>::const_iterator it= datamode.begin(); it != datamode.end() ; ++it) {  
    string suffix; 
    if(verbose_)  cout<<"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;
    }
    for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
      if (verbose_) cout << "Fill point of ntrk = " << ntrk << endl;
      PVStudy::PVAnaInfo ipv = GetPVAnaInfo(ntrk,*it);
      if (analyze_) {
        if ( ipv.res_.x() > 0 ) h2["resx"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
        if ( ipv.res_.y() > 0 ) h2["resy"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
        if ( ipv.res_.z() > 0 ) h2["resz"+suffix+"_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
        if ( ipv.pull_.x() > 0 ) h2["pullx"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
        if ( ipv.pull_.y() > 0 ) h2["pully"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
        if ( ipv.pull_.z() > 0 ) h2["pullz"+suffix+"_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
      }
    } // End of  for (int ntrk=nTrkMin_;ntrk<=nTrkMax_;++ntrk) {
    suffix.clear();
  }
  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
}

// Match two reconstructed vertices
bool PVStudy::matchVertex( const reco::Vertex & vrec1, 
                           const reco::Vertex & vrec2)
{
  return (fabs(vrec1.z()-vrec2.z())<0.0500); // =500um
}


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;
  return  ctau >0 && ctau <1e-6;
}

bool PVStudy::isFinalstateParticle(const HepMC::GenParticle * p){
  return ( !p->end_vertex() && p->status()==1 );
}

bool PVStudy::isCharged(const HepMC::GenParticle * p){
  const ParticleData * part = pdt->particle( p->pdg_id() );
  if (part){
    return part->charge()!=0;
  }else{
    // the new/improved particle table doesn't know anti-particles
    return  pdt->particle( -p->pdg_id() )!=0;
  }
}

void PVStudy::printSimVtxs(const Handle<SimVertexContainer> simVtxs){
  int i=0;
  for(SimVertexContainer::const_iterator vsim=simVtxs->begin();
      vsim!=simVtxs->end(); ++vsim){
    cout << i++ << ")" << scientific
         << " evtid=" << vsim->eventId().event() 
         << " sim x=" << vsim->position().x()
         << " sim y=" << vsim->position().y()
         << " sim z=" << vsim->position().z()
         << " sim t=" << vsim->position().t()
         << " parent=" << vsim->parentIndex() 
         << endl;
  }
}

void PVStudy::printRecVtxs(const Handle<reco::VertexCollection> recVtxs){
  int ivtx=0;
  for(reco::VertexCollection::const_iterator v=recVtxs->begin(); 
      v!=recVtxs->end(); ++v){
    cout << "Recvtx "<< std::setw(3) << std::setfill(' ')<<ivtx++
         << "#trk " << std::setw(3) << v->tracksSize() 
         << " chi2 " << std::setw(4) << v->chi2() 
         << " ndof " << std::setw(3) << v->ndof() << endl 
         << " x "  << std::setw(8) <<std::fixed << std::setprecision(4) << v->x() 
         << " dx " << std::setw(8) << v->xError()<< endl
         << " y "  << std::setw(8) << v->y() 
         << " dy " << std::setw(8) << v->yError()<< endl
         << " z "  << std::setw(8) << v->z() 
         << " dz " << std::setw(8) << v->zError()
         << endl;
  }
}

void PVStudy::printSimTrks(const Handle<SimTrackContainer> simTrks){
  cout <<  " simTrks   type, (momentum), vertIndex, genpartIndex"  << endl;
  int i=1;
  for(SimTrackContainer::const_iterator t=simTrks->begin();
      t!=simTrks->end(); ++t){
    //HepMC::GenParticle* gp=evtMC->GetEvent()->particle( (*t).genpartIndex() );
    cout << i++ << ")" 
         << (*t)
         << " index="
         << (*t).genpartIndex();
    //if (gp) {
    //  HepMC::GenVertex *gv=gp->production_vertex();
    //  cout  <<  " genvertex =" << (*gv);
    //}
    cout << endl;
  }
}

void PVStudy::fillHisto(res r, error er, int ntk, int datamode){
  stringstream ss;
  ss << ntk;
  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 (ntk < nTrkMin_ || ntk > nTrkMax_ ) return;
  // Fill histograms of res and pull of ntk
  h["deltax"+suffix]->Fill(r.x());
  h["deltay"+suffix]->Fill(r.y());
  h["deltaz"+suffix]->Fill(r.z()); 
  h["pullx"+suffix]->Fill(r.x()/er.x());
  h["pully"+suffix]->Fill(r.y()/er.y());
  h["pullz"+suffix]->Fill(r.z()/er.z());  
  h["errPVx"+suffix]->Fill(er.x());
  h["errPVy"+suffix]->Fill(er.y());
  h["errPVz"+suffix]->Fill(er.z());
}

PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk, int datamode) {
  map<int,double> data;
  data.clear();
  double fpar[2] = {-999,-999};
  double cm2um = 10000;
  stringstream ss;
  ss << ntk;
  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) {
      switch (i) {
      case 0:
        fit(h["deltax"+suffix], fpar);
        data[i] = fpar[0]*cm2um;
        break;
      case 1:
        fit(h["deltay"+suffix], fpar);
        data[i] = fpar[0]*cm2um;
        break;
      case 2:
        fit(h["deltaz"+suffix], fpar);
        data[i] = fpar[0]*cm2um;
        break;
      case 3:
        fit(h["pullx"+suffix], fpar);
        data[i] = fpar[0];
        break;
      case 4:
        fit(h["pully"+suffix], fpar);
        data[i] = fpar[0];
        break;
      case 5:
        fit(h["pullz"+suffix], fpar);
        data[i] = fpar[0];
        break;
      }
      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" ") << endl;
    }
  }
  else
    for ( int i = 0; i < 6; ++i) {
      data[i] = -999;
    }

  return PVStudy::PVAnaInfo(res(data[0], data[1], data[2]),
                            error(0.,0.,0.),
                            pull(data[3], data[4], data[5]),
                            error(0.,0.,0.),
                            ntk);
}

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);
  //   double fitRange = axis0->GetBinUpEdge(nbin);
  double fitRange = axis0->GetXmax();
  double sigMax[2] = {0.,0.};

  if ( verbose_ ){
    cout << "Last bin = " << nbin;
    cout << "; hdist: Overflow Entries = " << nOF;
    cout << "; Underflow Entries = " << nUF;
    cout << "; hdist->GetEntries() = " << hdist->GetEntries();
    cout << "; 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]);
      }
    }
    if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[1] << endl;

    TF1 *fgaus = new TF1("fgaus","gaus",-fitRange, fitRange);
    fgaus->SetParameter(1, 0.);
    fgaus->SetParLimits(1, -fitRange/40., fitRange/40.);
    fgaus->SetParLimits(2, 0., sqrt(2.)*sigMax[1]);
    fgaus->SetLineColor(4);
    hdist->Fit(fgaus,"Q");
    gPad->Update();
    TPaveStats *s = (TPaveStats*) hdist->GetListOfFunctions()->FindObject("stats");
    s->SetOptStat(1111111);
    s->SetOptFit(0111);
    gPad->Update();
    if (verbose_) cout << "got function" << endl;
    fitPars[0] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
  }
  else 
    fitPars[0] = -999;
}

void PVStudy::SetVarToZero() {
  
  // number of reconstructed vertices
  nrecPV_ = 0;
  nrecPV1_spl_ = 0;
  nrecPV2_spl_ = 0;
  // number of tracks in the vertex
  fntrk_ = 0;
  fntrk_spl1_mct_ = 0;
  fntrk_spl2_mct_ = 0;
  fntrk_mct_ = 0;
  //pvtx position (x,y,z) residual and error
  for(int i = 0; i<3;i++)
    {
      fres_[i] = 0;
      ferror_[i] = 0;
      fres_spl1_mct_[i] = 0;
      ferror_spl1_mct_[i] = 0;
      fres_spl2_mct_[i] = 0;
      ferror_spl2_mct_[i] = 0;
      fres_mct_[i] = 0;
      ferror_mct_[i] = 0;
    }
}


Revision:	1.7
Committed:	Thu Oct 8 14:57:46 2009 UTC (15 years, 7 months ago) by yygao
Content type:	text/plain
Branch:	MAIN
CVS Tags:	CMSSW_3_1_2_20091015, CMSSW_3_1_2_Yanyan
Changes since 1.6:	+636 -224 lines
Log Message:	Update pvtxstudy with the two-vertices method