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$
//
//


// 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"


using namespace std;

PVStudy::PVStudy(const edm::ParameterSet& iConfig):
  nTrkMin_(10),nTrkMax_(150)
{
  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);

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

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

  h_nTrk           = fs->make<TH1F>("nTrk","Num of total tracks",300,0,300);
  h_nTrkPV         = fs->make<TH1F>("nTrkPV","Num of Tracks from PV",300,0,300);
  h_trkPt          = fs->make<TH1D>("trkPt","Pt of all tracks",100,0,100);
  h_trkPtPV        = fs->make<TH1D>("trkPtPV","Pt dist of tracks from PV",100,0,100);

  if (!realData_) {
    h_deltax       = fs->make<TH1D>("deltax","x-distance (Rec - Sim)",800,-0.04,0.04);
    h_deltay       = fs->make<TH1D>("deltay","y-distance (Rec - Sim)",800,-0.04,0.04);
    h_deltaz       = fs->make<TH1D>("deltaz","z-distance (Rec - Sim)",800,-0.04,0.04);
    h_pullx        = fs->make<TH1D>("pullx","x-pull (Rec - Sim)",1000,-25.,25.);
    h_pully        = fs->make<TH1D>("pully","y-pull (Rec - Sim)",1000,-25.,25.);
    h_pullz        = fs->make<TH1D>("pullz","z-pull (Rec - Sim)",1000,-25.,25.);
    // Summary of analysis:
    h_resx_nTrk    = fs->make<TH2D>("resx_ntrk","x-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
    h_resx_nTrk->SetMarkerStyle(21);
    h_resx_nTrk->GetYaxis()->SetTitle("#mum");
    h_resy_nTrk    = fs->make<TH2D>("resy_ntrk","y-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
    h_resy_nTrk->SetMarkerStyle(21);
    h_resy_nTrk->GetYaxis()->SetTitle("#mum");
    h_resz_nTrk    = fs->make<TH2D>("resz_ntrk","z-resolution (Rec - Sim) vs number of tracks",150,0,150,400,0.,200);
    h_resz_nTrk->SetMarkerStyle(21);
    h_resz_nTrk->GetYaxis()->SetTitle("#mum");
    h_pullx_nTrk   = fs->make<TH2D>("pullx_ntrk","x-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
    h_pullx_nTrk->SetMarkerStyle(21);
    h_pullx_nTrk->SetBit(TH1::kCanRebin);
    h_pullx_nTrk->GetYaxis()->SetTitle("cm");
    h_pully_nTrk   = fs->make<TH2D>("pully_ntrk","y-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
    h_pully_nTrk->SetMarkerStyle(21);
    h_pully_nTrk->SetBit(TH1::kCanRebin);
    h_pully_nTrk->GetYaxis()->SetTitle("cm");
    h_pullz_nTrk   = fs->make<TH2D>("pullz_ntrk","z-pull (Rec - Sim) vs number of tracks",150,0,150,100,0.,2.);
    h_pullz_nTrk->SetMarkerStyle(21);
    h_pullz_nTrk->SetBit(TH1::kCanRebin);
    h_pullz_nTrk->GetYaxis()->SetTitle("cm");

    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)",400,-0.02,0.02);
      h["resx_"+suffix]->GetXaxis()->SetTitle("cm");
      h["resy_"+suffix]  = fs->make<TH1D> (TString("deltay_"+suffix),"y-distance (Rec - Sim)",400,-0.02,0.02);
      h["resy_"+suffix]->GetXaxis()->SetTitle("cm");
      h["resz_"+suffix]  = fs->make<TH1D> (TString("deltaz_"+suffix),"z-distance (Rec - Sim)",400,-0.02,0.02);
      h["resz_"+suffix]->GetXaxis()->SetTitle("cm");
      h["pullx_"+suffix] = fs->make<TH1D> (TString("pullx_"+suffix),"x-pull (Rec - Sim)",200,-10.,10.);
      h["pullx_"+suffix]->GetXaxis()->SetTitle("cm");
      h["pully_"+suffix] = fs->make<TH1D> (TString("pully_"+suffix),"y-pull (Rec - Sim)",200,-10.,10.);
      h["pully_"+suffix]->GetXaxis()->SetTitle("cm");
      h["pullz_"+suffix] = fs->make<TH1D> (TString("pullz_"+suffix),"z-pull (Rec - Sim)",200,-10.,10.);
      h["pullz_"+suffix]->GetXaxis()->SetTitle("cm");
      suffix.clear();
    }
  }
}


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

}


//
// 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;
  
  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;
          }
          h_deltax->Fill( vsim->recVtx->x()-vsim->x );
          h_deltay->Fill( vsim->recVtx->y()-vsim->y );
          h_deltaz->Fill( vsim->recVtx->z()-vsim->z );
          h_pullx->Fill( (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError() );
          h_pully->Fill( (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError() );
          h_pullz->Fill( (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError() );
          pvanainfo.push_back(PVStudy::PVAnaInfo( vsim->recVtx->x()-vsim->x,
                                                  vsim->recVtx->y()-vsim->y,
                                                  vsim->recVtx->z()-vsim->z,
                                                  (vsim->recVtx->x()-vsim->x )/vsim->recVtx->xError(),
                                                  (vsim->recVtx->y()-vsim->y )/vsim->recVtx->yError(),
                                                  (vsim->recVtx->z()-vsim->z )/vsim->recVtx->zError(),
                                                  nrectrk )
                              );
        }
        else{  // no rec vertex found for this simvertex
          if(verbose_) {
            cout <<"primary not found " << vsim->x << " " << vsim->y << " " << vsim->z << " nGenTrk=" << vsim->nGenTrk << endl;
          }
        }
      }
    }
  }// 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 ( ipv.resx > 0 ) h_resx_nTrk->Fill(ntrk,ipv.resx/sqt2, 1.);
      if ( ipv.resy > 0 ) h_resy_nTrk->Fill(ntrk,ipv.resy/sqt2, 1.);
      if ( ipv.resz > 0 ) h_resz_nTrk->Fill(ntrk,ipv.resz/sqt2, 1.);
      if ( ipv.pullx > 0 ) h_pullx_nTrk->Fill(ntrk,ipv.pullx, 1.);
      if ( ipv.pully > 0 ) h_pully_nTrk->Fill(ntrk,ipv.pully, 1.);
      if ( ipv.pullz > 0 ) h_pullz_nTrk->Fill(ntrk,ipv.pullz, 1.);
    }
  }
}

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

PVStudy::PVAnaInfo PVStudy::GetPVAnaInfo(int ntk) {
  map<int,double> data;
  double xmin=0.,xmax=0.;
  double fpar[2];
  double cm2um = 10000;
  stringstream ss;
  ss << ntk;
  string suffix = ss.str();

  for ( int i = 0; i < 6; ++i) {
    switch(i){
    case 0:
    case 1:
    case 2:
      if (verbose_) cout << "Analyzing Resolution:" << endl;
      xmin = -0.01; xmax = 0.01;
      break;
    default:
      if (verbose_) cout << "Analyzing Pull:" << endl;
      xmin = -10.; xmax = 10;
      break;
    }
    TH1D *hdist = new TH1D("hdist","distribution to be fitted", 1000,xmin,xmax);

    for (std::vector<PVStudy::PVAnaInfo>::const_iterator ipvinfo = pvanainfo.begin();
         ipvinfo != pvanainfo.end(); ++ipvinfo){
      if ( verbose_ && i==0 ) {
        cout << "Num of tracks of the PV = " << ipvinfo->nTrk;
        cout << "; delta x = " << ipvinfo->resx << endl;
      }
      if (ipvinfo->nTrk == ntk) {
        switch(i) {
        case 0: if (verbose_) cout << "Matched resx" << endl; 
          hdist->Fill(ipvinfo->resx);
          h["resx_"+suffix]->Fill(ipvinfo->resx);
          break;
        case 1: if (verbose_) cout << "Matched resy" << endl;
          hdist->Fill(ipvinfo->resy);
          h["resy_"+suffix]->Fill(ipvinfo->resy);
          break;
        case 2: if (verbose_) cout << "Matched resz" << endl;
          hdist->Fill(ipvinfo->resz);
          h["resz_"+suffix]->Fill(ipvinfo->resz);
          break;
        case 3: if (verbose_) cout << "Matched pullx" << endl;
          hdist->Fill(ipvinfo->pullx);
          h["pullx_"+suffix]->Fill(ipvinfo->pullx);
          break;
        case 4: if (verbose_) cout << "Matched pully" << endl;
          hdist->Fill(ipvinfo->pully);
          h["pully_"+suffix]->Fill(ipvinfo->pully);
          break;
        case 5: if (verbose_) cout << "Matched pullz" << endl;
          hdist->Fill(ipvinfo->pullz);
          h["pullz_"+suffix]->Fill(ipvinfo->pullz);
          break;
        }
      }
      else {
        //      switch(i) {
        //      case 0: if (verbose_) cout << "bad resx" << endl; break;
        //      case 1: if (verbose_) cout << "bad resy" << endl; break;
        //      case 2: if (verbose_) cout << "bad resz" << endl; break;
        //      case 3: if (verbose_) cout << "bad pullx" << endl; break;
        //      case 4: if (verbose_) cout << "bad pully" << endl; break;
        //      case 5: if (verbose_) cout << "bad pullz" << endl; break;
        //      }
      }
    }
    TAxis *axis0 = hdist->GetXaxis();
    int nbin = axis0->GetLast();
    double nOF = hdist->GetBinContent(nbin+1);
    double nUF = hdist->GetBinContent(0);
    if ( verbose_ && i==0 ){
      cout << "Last bin = " << nbin;
      cout << "; hdist: Overflow Entries = " << nOF;
      cout << "; Underflow Entries = " << nUF;
      cout << "; hdist->GetEntries() = " << hdist->GetEntries() << endl;
    }
    if (hdist->GetEntries() - nOF - nUF >= 5) { // FIXME: for reasonable Gaussian fit
      if(i<3){
        hdist->Fit("gaus","Q0");
        TF1 *fgaus = hdist->GetFunction("gaus");
        if (verbose_) cout << "got function" << endl;
        fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
        fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
        delete fgaus;
      }
      else {
        hdist->Fit("gausn","Q0");
        TF1 *fgaus = hdist->GetFunction("gausn");
        if (verbose_) cout << "got function" << endl;
        fpar[0] = ((fgaus->GetParameter(1))?fgaus->GetParameter(1):-999);
        fpar[1] = ((fgaus->GetParameter(2))?fgaus->GetParameter(2):-999);
        delete fgaus;
      }
      if (verbose_) cout << "fpar[2] = (" << fpar[0] << "," << fpar[1] << ")" << endl;
      switch (i) {
      case 0:
        h["resx_"+suffix]->Fit("gaus","Q");
        data[i] = fpar[1]*cm2um;
        break;
      case 1:
        h["resy_"+suffix]->Fit("gaus","Q");
        data[i] = fpar[1]*cm2um;
        break;
      case 2:
        h["resz_"+suffix]->Fit("gaus","Q");
        data[i] = fpar[1]*cm2um;
        break;
      case 3:
        h["pullx_"+suffix]->Fit("gausn","Q");
        data[i] = fpar[1];
        break;
      case 4:
        h["pully_"+suffix]->Fit("gausn","Q");
        data[i] = fpar[1];
        break;
      case 5:
        h["pullz_"+suffix]->Fit("gausn","Q");
        data[i] = fpar[1];
        break;
      }
      if (verbose_ && data[i] > 0) cout << "data[" << i << "] = " << data[i] << (i<3?" micro m":" cm") << endl;
    }
    else {
      if ( verbose_) cout << "Empty Histo! No Fit!!!" << endl;
      data[i] = -999;
    }
    delete hdist;
  }

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

//define this as a plug-in
DEFINE_FWK_MODULE(PVStudy);
Revision:	1.1
Committed:	Thu Sep 3 17:31:08 2009 UTC (15 years, 8 months ago) by jengbou
Content type:	text/plain
Branch:	MAIN
Branch point for:	CMSSW_2_2_X
Log Message:	Add file