PVStudy/plugins/PVEffAnalyzer.cc


// -*- C++ -*-
//
// Package:    PVEffAnalyzer
// Class:      PVEffAnalyzer
// 
/**\class PVEffAnalyzer PVEffAnalyzer.cc UserCode/PVEffAnalyzer/plugins/PVEffAnalyzer.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: PVEffAnalyzer.cc,v 1.1 2010/04/19 16:07:17 yygao 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 "CommonTools/UtilAlgos/interface/TFileService.h"
#include "UserCode/PVStudy/interface/PVEffAnalyzer.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>

// BeamSpot 
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
// For the clusters
#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
#include "TrackingTools/Records/interface/TransientTrackRecord.h"
#include "TrackingTools/TransientTrack/interface/TransientTrack.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 

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

using namespace std;
using namespace edm;
using namespace reco;

typedef math::XYZTLorentzVectorF LorentzVector;
typedef math::XYZPoint Point;

PVEffAnalyzer::PVEffAnalyzer(const edm::ParameterSet& iConfig)
  :theTrackClusterizer(iConfig.getParameter<edm::ParameterSet>("TkClusParameters"))
{
  //=======================================================================
  // 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);
  histoFileName_             = iConfig.getUntrackedParameter<std::string> ("histoFileName");
  nTrkMin_                   = iConfig.getUntrackedParameter<int>("nTrkMin");
  nTrkMax_                   = iConfig.getUntrackedParameter<int>("nTrkMax");
  zsigncut_                  = iConfig.getUntrackedParameter<double>("zsigncut");
  analyze_                   = iConfig.getUntrackedParameter<bool>("analyze",false);
  bsSrc = iConfig.getParameter< edm::InputTag >("beamSpot");
  reqCluster_               = iConfig.getUntrackedParameter<bool>("reqCluster",false);
    
  // Specify the data mode vector
  if(realData_) datamode.push_back(0);
  else {
    datamode.push_back(0);
    datamode.push_back(1);
  }
  
  theFile = new TFile(histoFileName_.c_str(), "RECREATE"); 
  theFile->mkdir("Summary"); 
  theFile->cd();

  // Book MC only plots
  if (!realData_) {
    h_gen = new PVEffHistograms();
    h_gen->Init("generator");
  }
  
  h_summary = new PVEffHistograms();
  //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 = "_mct";
      break;
    }
    h_summary->Init("summary",suffix, nTrkMin_, nTrkMax_);
  }

}

PVEffAnalyzer::~PVEffAnalyzer()
{
  // do anything here that needs to be done at desctruction time
  // (e.g. close files, deallocate resources etc.)
  theFile->cd();
  theFile->cd("Summary");
  h_summary->Save();
  if (!realData_) 
    h_gen->Save();
  theFile->Close();
}

//
// member functions
//
std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::getSimPVs(const Handle<HepMCProduct> evtMC, std::string suffix="")
{
 std::vector<PVEffAnalyzer::simPrimaryVertex> simpv;
  const HepMC::GenEvent* evt=evtMC->GetEvent();
  if (evt) {
    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 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 ) {
        //        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(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++;
            }
          }
        }
      }//loop MC vertices daughters
    }//loop MC vertices
  }
  return simpv;
}

std::vector<PVEffAnalyzer::simPrimaryVertex> PVEffAnalyzer::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<PVEffAnalyzer::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 );
          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())));
          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<PVEffAnalyzer::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){
            PVEffAnalyzer::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
PVEffAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;
  using namespace reco;
  
  //========================================================================
  // Step 0: Prepare root variables and get information from the Event
  //========================================================================
  
  edm::LogInfo("Debug")<<"[PVEffAnalyzer]"<<endl;

  // ====== 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") << "[PVEffAnalyzer] 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") << "[PVEffAnalyzer] 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") << "[PVEffAnalyzer] splitTrackCollection2 cannot be found -> using empty collection of same type." <<endl;
  }
  
  // ======= PrimaryVertexCollection
  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") << "[PVEffAnalyzer] 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") << "[PVEffAnalyzer] 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") << "[PVEffAnalyzer] splitVertexCollection2 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")<<"[PVEffAnalyzer] End accessing the track, beamSpot, primary 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&){
    edm::LogInfo("Debug") << "[PVEffAnalyzer] Some problem occurred with the particle data table. This may not work !." <<endl;
  }
  
  //setUpVectors(RealData, nTrkMin_, nTrkMax_ ) ;
  
  // ======= Analyze MC efficiency
  if (!realData_) { 
    bool MC=false;  
    Handle<HepMCProduct> evtMC;
    iEvent.getByLabel("generator",evtMC);
    if (!evtMC.isValid()) {
      MC=false;
      edm::LogInfo("Debug") << "[PVEffAnalyzer] no HepMCProduct found"<< endl;
    } else {
      edm::LogInfo("Debug") << "[PVEffAnalyzer] generator HepMCProduct found"<< endl;
      MC=true;
    }
    if(MC){
      TString suffix = "_mct";
      // make a list of primary vertices:
      std::vector<simPrimaryVertex> simpv;
      simpv=getSimPVs(evtMC,"");
      //   simpv=getSimPVs(evtMC, simVtxs, simTrks);
      h_gen->Fill1d("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;
        if(verbose_ && simpv.size() > 1 )
          std::cout<<"Simulated Vertex # " << isimpv << ": ptsq = " << vsim->ptsq<<std::endl;
      }
      
      // Just analyze the first vertex in the collection 
      simPrimaryVertex vsim = *simpv.begin();
      h_summary->Fill1d(TString("denom_ntrack"+suffix), int(vsim.nGenTrk));
      if ( isAssoVertex(vsim, *vertexColl->begin(), zsigncut_)  )
        h_summary->Fill1d(TString("numer_ntrack"+suffix), int(vsim.nGenTrk));
    }
  } // End of Analyzing MC Efficiency
  
  // Get the Builder for TrackClusters
  edm::ESHandle<TransientTrackBuilder> theB;
  iSetup.get<TransientTrackRecord>().get("TransientTrackBuilder",theB); 
  
  std::vector< std::vector<reco::TransientTrack> > clusters1 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection1Handle));
  std::vector< std::vector<reco::TransientTrack> > clusters2 =  theTrackClusterizer.clusterize((*theB).build(splitTrackCollection2Handle));
  
  bool withValidCluster_tag = false;
  bool withValidCluster_probe = false;

  for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters1.begin(); iclus != clusters1.end(); iclus++) {
    if((*iclus).size()>1)  withValidCluster_probe = true;
  }
  for (vector< vector<reco::TransientTrack> >::const_iterator iclus = clusters2.begin(); iclus != clusters2.end(); iclus++) {
    if((*iclus).size()>1)  withValidCluster_tag = true;
  }

  // ======= Analyze efficiency with split method
  if( reqCluster_ ) {
    if ( !isGoodSplitEvent( *vertexColl->begin()) || !withValidCluster_tag || !withValidCluster_probe ) return;
  }
  else {
    if ( !isGoodSplitEvent( *vertexColl->begin())) return;
  }
  if ( isGoodTagVertex( *splitVertexColl2->begin(), *vertexColl->begin(), zsigncut_) ) {
    if(verbose_) {
      std::cout<<"splitTrackColl1->size() = " << int(splitTrackColl1->size()) << std::endl;
      std::cout<<"splitTrackColl2->size() = " << int(splitTrackColl2->size()) << std::endl;
    }
    h_summary->Fill1d("denom_ntrack", int(splitTrackColl1->size()));
    if ( isGoodProbeVertex( *splitVertexColl1->begin(), *vertexColl->begin(), zsigncut_) ) 
      h_summary->Fill1d("numer_ntrack", int(splitTrackColl1->size()) );
  }
}

  
// ------------ method called once each job just before starting event loop  ------------
void 
PVEffAnalyzer::beginJob()
{
}

// ------------ method called once each job just after ending the event loop  ------------
void 
PVEffAnalyzer::endJob() {
  edm::LogInfo("Analysis") << "[endJob] Analyzing PV info" << endl;
  
  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 = "_mct";
      break;
    }
    if (analyze_)
      MakeEff(h_summary->ReadHisto1D(TString("numer_ntrack"+suffix)), h_summary->ReadHisto1D(TString("denom_ntrack"+suffix)), h_summary->ReadHisto1D(TString("eff_ntrack"+suffix)), false, 1);
  }    
}


bool PVEffAnalyzer::isResonance(const HepMC::GenParticle * p){
  double ctau=(pdt->particle( abs(p->pdg_id()) ))->lifetime();
  edm::LogInfo("Debug") << "[isResonance] isResonance   " << p->pdg_id() << " " << ctau << endl;
  return  ctau >0 && ctau <1e-6;
}

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

bool PVEffAnalyzer::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 PVEffAnalyzer::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 PVEffAnalyzer::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 PVEffAnalyzer::printRecVtx(const reco::Vertex & v){

  cout << "#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 PVEffAnalyzer::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;
  }
}

// Select based on the highest SumPtSq Vertex 
// Require at least 6 tracks with track weight > ~  0.8
// vtx.ndof = 2*Sum(trackWeight) - 3 

bool PVEffAnalyzer::isGoodSplitEvent( const reco::Vertex & org_vtx) 
//, const reco::Vertex & tag_vtx, const reco::Vertex & probe_vtx)
{
  return (  org_vtx.tracksSize() >= 6 && org_vtx.ndof() >= 6 );
            //      tag_vtx.tracksSize() > = 2 && probe_vtx.tracksSize() > = 2      ) ;
}

// Comparing the tag vertex with the unsplit vertex position in z
bool PVEffAnalyzer::isGoodTagVertex( const reco::Vertex & tag_vtx, const reco::Vertex & org_vtx, double & zsign_cut) 
{
  if(tag_vtx.isValid () && !tag_vtx.isFake()) {
    float zsign = TMath::Abs(tag_vtx.z() - org_vtx.z() ) / TMath::Max(tag_vtx.zError(), org_vtx.zError() );
    return (zsign < zsign_cut) ;
  }
  else
    return false;
}

// Comparing the probe vertex with the unsplit vertex position in z
bool PVEffAnalyzer::isGoodProbeVertex( const reco::Vertex & probe_vtx, const reco::Vertex & org_vtx, double & zsign_cut) 
{
  if(probe_vtx.isValid () && !probe_vtx.isFake()) {
    float zsign = TMath::Abs(probe_vtx.z() - org_vtx.z() ) / TMath::Max(probe_vtx.zError(), org_vtx.zError() );
    return (zsign < zsign_cut) ;
  }
  else
    return false;
}


bool PVEffAnalyzer::isAssoVertex(const PVEffAnalyzer::simPrimaryVertex & vsim, 
                                 const reco::Vertex & vrec, double & zsign_cut)
{
  if(vrec.isValid() && !vrec.isFake() )
    return (TMath::Abs(vsim.z-vrec.z()) < zsign_cut * vrec.zError() );
  else 
    return false;
}


void PVEffAnalyzer::MakeEff(TH1D* numer, TH1D* denom, TH1D* eff,  //TGraphAsymmErrors* & gr_eff,
                            const bool rebin, const Float_t n ) {
  eff->Divide( numer, denom, 1,1,"B" );
  
  if( rebin ) {
    std::vector<Double_t> binEdges;
    Int_t bin = denom->GetNbinsX() + 1;
    binEdges.push_back(denom->GetBinLowEdge(bin));
    Float_t nEntries = 0;
    while (bin > 1) {
      bin --;
      nEntries = denom->GetBinContent(bin);
      while (nEntries < n && bin > 1) {
        bin --;
        nEntries += denom->GetBinContent(bin);
      }
      binEdges.push_back(denom->GetBinLowEdge(bin));
    }
    
    Double_t *array = new Double_t[binEdges.size()];
    
    Int_t j = 0;
    for (Int_t i = binEdges.size(); i > 0; --i) {
      array[j] = binEdges[i - 1];
      ++j;
    }
  }
  
  for (int i = 1; i<eff->GetNbinsX()+1; i++) {
    if(numer->GetBinContent(i) == 0 || denom->GetBinContent(i) == 0 ) continue;
    float error = eff->GetBinContent(i)*sqrt(pow(numer->GetBinError(i)/numer->GetBinContent(i),2)
                                             +pow(denom->GetBinError(i)/denom->GetBinContent(i),2));
    eff->SetBinError(i, error);
  }
}
Revision:	1.2
Committed:	Tue Apr 20 09:51:02 2010 UTC (15 years ago) by yygao
Content type:	text/plain
Branch:	MAIN
Changes since 1.1:	+40 -31 lines
Log Message:	add the pre-requirement of cluster