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"
//         Created:  Thu Aug 20 11:55:40 CDT 2009
// $Id: PVStudy.cc,v 1.3 2009/09/08 07:15: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):
  nTrkMin_(10),nTrkMax_(100)
{
  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");

  //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 subDir1 = fs->mkdir( "Others" );
  //   TFileDirectory subSubDir = subDir.mkdir( "mySubSubDirectory" );

  setRootStyle();

  h["nTrk"]           = subDir.make<TH1F>("nTrk","Num of total tracks",300,0,300);
  h["nTrkPV"]         = subDir.make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
  h["trkPt"]          = subDir.make<TH1D>("trkPt","Pt of all tracks",100,0,100);
  h["trkPtPV"]        = subDir.make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);

  if (!realData_) {
    h["deltax"]       = subDir.make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
    h["deltay"]       = subDir.make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
    h["deltaz"]       = subDir.make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
    h["pullx"]        = subDir.make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
    h["pully"]        = subDir.make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
    h["pullz"]        = subDir.make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);

    // Summary of analysis:
    if (analyze_) {
      h2["resx_nTrk"]= subDir.make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
      h2["resx_nTrk"]->SetMarkerStyle(21);
      h2["resx_nTrk"]->SetMarkerColor(4);
      h2["resx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resx_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["resy_nTrk"]= subDir.make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
      h2["resy_nTrk"]->SetMarkerStyle(21);
      h2["resy_nTrk"]->SetMarkerColor(4);
      h2["resy_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resy_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["resz_nTrk"]= subDir.make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
      h2["resz_nTrk"]->SetMarkerStyle(21);
      h2["resz_nTrk"]->SetMarkerColor(4);
      h2["resz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["resz_nTrk"]->GetYaxis()->SetTitle("#mum");
      h2["pullx_nTrk"]= subDir.make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
      h2["pullx_nTrk"]->SetMarkerStyle(21);
      h2["pullx_nTrk"]->SetMarkerColor(4);
      h2["pullx_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pullx_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["pully_nTrk"]= subDir.make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
      h2["pully_nTrk"]->SetMarkerStyle(21);
      h2["pully_nTrk"]->SetMarkerColor(4);
      h2["pully_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pully_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
      h2["pullz_nTrk"]= subDir.make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
      h2["pullz_nTrk"]->SetMarkerStyle(21);
      h2["pullz_nTrk"]->SetMarkerColor(4);
      h2["pullz_nTrk"]->SetBit(TH1::kCanRebin);
      h2["pullz_nTrk"]->GetXaxis()->SetTitle("Num of tracks");
    }
    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.);
      h["errPVx_"+suffix]  = subDir1.make<TH1D> (TString("errPVx_"+suffix),"X vertex error",100,0.,0.02);
      h["errPVy_"+suffix]  = subDir1.make<TH1D> (TString("errPVy_"+suffix),"Y vertex error",100,0.,0.02);
      h["errPVz_"+suffix]  = subDir1.make<TH1D> (TString("errPVz_"+suffix),"Z vertex error",100,0.,0.02);
      suffix.clear();
    }
  }
}


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

        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++;
              }
            }
          }
        }
        idx++;
      }
  }
  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 reco::TrackCollection;

  ftree_->SetBranchAddress("residual",&fres);
  ftree_->SetBranchAddress("error",&ferror);
  ftree_->SetBranchAddress("nTrk",&fntrk);

  Handle<TrackCollection> tracks;
  iEvent.getByLabel(tracksLabel_,tracks);

  Handle<reco::VertexCollection> recVtxs;
  iEvent.getByLabel(vtxSample_, recVtxs);

  if (!realData_) {
    Handle<SimVertexContainer> simVtxs;
    iEvent.getByLabel( simG4_, simVtxs);
    
    Handle<SimTrackContainer> simTrks;
    iEvent.getByLabel( simG4_, simTrks);
  }

  if(verbose_) printRecVtxs(recVtxs);

  try{
    iSetup.getData(pdt);
  }catch(const Exception&){
    cout << "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;

          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] );
          pvinfo.push_back(PVStudy::PVInfo(res(fres[0], fres[1], fres[2]),
                                           error(ferror[0], ferror[1], ferror[2]),
                                           nrectrk));
          // Fill histo according to its track multiplicity
          fillHisto(res(fres[0], fres[1], fres[2]),
                    error(ferror[0], ferror[1], ferror[2]),
                    nrectrk);

          ftree_->Fill();
        }
        else{  // no rec vertex found for this simvertex
          if(verbose_) {
            cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
          }
        }
      }
    }
  }// 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());
        }
      }
    }
    catch (...) {
      // exception thrown when trying to use linked track
      //          h["trkPtPV"]->Fill(0.);
    }
    
    h["nTrkPV"]->Fill(v->tracksSize());
    
  }
  
  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());
  }
  
}


// ------------ 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 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x()/sqt2, 1.);
//      if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y()/sqt2, 1.);
//      if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z()/sqt2, 1.);
        if ( ipv.res_.x() > 0 ) h2["resx_nTrk"]->Fill(ntrk,ipv.res_.x(), 1.);
        if ( ipv.res_.y() > 0 ) h2["resy_nTrk"]->Fill(ntrk,ipv.res_.y(), 1.);
        if ( ipv.res_.z() > 0 ) h2["resz_nTrk"]->Fill(ntrk,ipv.res_.z(), 1.);
        if ( ipv.pull_.x() > 0 ) h2["pullx_nTrk"]->Fill(ntrk,ipv.pull_.x(), 1.);
        if ( ipv.pull_.y() > 0 ) h2["pully_nTrk"]->Fill(ntrk,ipv.pull_.y(), 1.);
        if ( ipv.pull_.z() > 0 ) h2["pullz_nTrk"]->Fill(ntrk,ipv.pull_.z(), 1.);
      }
    }
  }
  file_->cd();
  ftree_->Write();
  file_->Close();
}

bool PVStudy::matchVertex(const PVStudy::simPrimaryVertex & vsim, 
                          const reco::Vertex & vrec)
{
  return (fabs(vsim.z-vrec.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){
  stringstream ss;
  ss << ntk;
  string suffix = 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(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) {
  map<int,double> data;
  data.clear();
  double fpar[2] = {-999,-999};
  double cm2um = 10000;
  stringstream ss;
  ss << ntk;
  string suffix = ss.str();

  if(analyze_) {
    for ( int i = 0; i < 6; ++i) {
      switch (i) {
      case 0:
        fit(h["resx_"+suffix], fpar);
        data[i] = fpar[0]*cm2um;
        break;
      case 1:
        fit(h["resy_"+suffix], fpar);
        data[i] = fpar[0]*cm2um;
        break;
      case 2:
        fit(h["resz_"+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[3] = {0.,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[2]) ) sigMax[2] = abs(sigMax[0]);
      }
      if ( (axis0->GetBinLowEdge(bi) >= 0) && (hdist->GetBinContent(bi) > 0) ) {
        sigMax[1] = axis0->GetBinUpEdge(bi);
        if ( abs(sigMax[1]) > abs(sigMax[2]) ) sigMax[2] = abs(sigMax[1]);
      }
    }
    if (verbose_) cout << "Fit sigMax = " << sqrt(2.)*sigMax[2] << 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[2]);
    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;
}

//define this as a plug-in
DEFINE_FWK_MODULE(PVStudy);
Revision:	1.4
Committed:	Tue Sep 8 21:48:16 2009 UTC (15 years, 8 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+154 -96 lines
Log Message:	update