LhcTrackAnalyzer.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:    LhcTrackAnalyzer
// Class:      LhcTrackAnalyzer
// 
//

// updated to 25/2/2009 5.30 pm

//
//

// system include files
#include <memory>


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

#include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
#include "FWCore/Framework/interface/Event.h"
#include <FWCore/Framework/interface/ESHandle.h>
#include "FWCore/Framework/interface/MakerMacros.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"


#include "TH1F.h"
#include "TH2F.h"
#include "TFile.h"
#include "TROOT.h"
#include "TChain.h"
#include "TNtuple.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
#include <Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h>
#include <DataFormats/GeometrySurface/interface/Surface.h>
#include <DataFormats/GeometrySurface/interface/GloballyPositioned.h>
#include <Geometry/CommonDetUnit/interface/GeomDet.h>
#include "MagneticField/Engine/interface/MagneticField.h" 
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" 

#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "DataFormats/VertexReco/interface/VertexFwd.h"

// Constructor

LhcTrackAnalyzer::LhcTrackAnalyzer(const edm::ParameterSet& iConfig)

{
  //now do what ever initialization is needed
  debug_    = iConfig.getParameter<bool>       ("Debug");  
  TrackCollectionTag_      = iConfig.getParameter<edm::InputTag>("TrackCollectionTag");  
  PVtxCollectionTag_      = iConfig.getParameter<edm::InputTag>("PVtxCollectionTag");  
  filename_ = iConfig.getParameter<std::string>("OutputFileName");
}
   
// Destructor
LhcTrackAnalyzer::~LhcTrackAnalyzer()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
LhcTrackAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
  using namespace edm;
  using namespace reco;
  using namespace std;

  //=======================================================
  // Initialize Root-tuple variables
  //=======================================================

  SetVarToZero();
  
  //=======================================================
  // Retrieve the Track information
  //=======================================================
  
  Handle< TrackCollection>  trackCollectionHandle;
  iEvent.getByLabel(TrackCollectionTag_, trackCollectionHandle);
  Handle<VertexCollection>  vertexCollectionHandle;
  iEvent.getByLabel(PVtxCollectionTag_, vertexCollectionHandle);
  for(VertexCollection::const_iterator vtx = vertexCollectionHandle->begin();vtx!=vertexCollectionHandle->end(); ++vtx)
    {
      if(vtx==vertexCollectionHandle->begin()){
    if(vtx->isFake())goodvtx_=false;
    else goodvtx_=true;
      }
      else break;
    }



  goodbx_=true;
  // int bx = iEvent.bunchCrossing();
  //if (bx==51 || bx==2724) goodbx_=true;
  

 
  run_=iEvent.id().run();
  event_=iEvent.id().event();
 
  if(debug_)
    cout<<"LhcTrackAnalyzer::analyze() looping over "<< trackCollectionHandle->size()<< "tracks." << endl;    
  
  // unsigned int i = 0;   
  for(TrackCollection::const_iterator track = trackCollectionHandle->begin(); track!= trackCollectionHandle->end(); ++track)
    {
      if ( nTracks_ >= nMaxtracks_ ) {
    std::cout << " LhcTrackAnalyzer::analyze() : Warning - Run "<< run_<<" Event "<< event_<<"\tNumber of tracks: " <<  trackCollectionHandle->size()<< " , greater than " << nMaxtracks_ << std::endl;
      continue;
     
    }
      pt_[nTracks_]       = track->pt();
      eta_[nTracks_]      = track->eta();
      phi_[nTracks_]      = track->phi();
      chi2_[nTracks_]     = track->chi2();
      chi2ndof_[nTracks_] = track->normalizedChi2();
      charge_[nTracks_]   = track->charge();
      qoverp_[nTracks_]   = track->qoverp();
      dz_[nTracks_]       = track->dz();
      dxy_[nTracks_]      = track->dxy();
      xPCA_[nTracks_]     = track->vertex().x();
      yPCA_[nTracks_]     = track->vertex().y();
      zPCA_[nTracks_]     = track->vertex().z(); 
      validhits_[nTracks_][0]=track->numberOfValidHits();
      validhits_[nTracks_][1]=track->hitPattern().numberOfValidPixelBarrelHits();
      validhits_[nTracks_][2]=track->hitPattern().numberOfValidPixelEndcapHits();
      validhits_[nTracks_][3]=track->hitPattern().numberOfValidStripTIBHits();
      validhits_[nTracks_][4]=track->hitPattern().numberOfValidStripTIDHits();
      validhits_[nTracks_][5]=track->hitPattern().numberOfValidStripTOBHits();
      validhits_[nTracks_][6]=track->hitPattern().numberOfValidStripTECHits();



      int myalgo=-88;
      if(track->algo()==reco::TrackBase::undefAlgorithm)myalgo=0;
      if(track->algo()==reco::TrackBase::ctf)myalgo=1;
      if(track->algo()==reco::TrackBase::iter0)myalgo=4;
      if(track->algo()==reco::TrackBase::iter1)myalgo=5;
      if(track->algo()==reco::TrackBase::iter2)myalgo=6;
      if(track->algo()==reco::TrackBase::iter3)myalgo=7;
      if(track->algo()==reco::TrackBase::iter4)myalgo=8;
      if(track->algo()==reco::TrackBase::iter5)myalgo=9;
      if(track->algo()==reco::TrackBase::iter6)myalgo=10;
      if(track->algo()==reco::TrackBase::iter7)myalgo=11;
      trkAlgo_[nTracks_]  = myalgo;

      int myquality=-99;
      if(track->quality(reco::TrackBase::undefQuality))myquality=-1;
      if(track->quality(reco::TrackBase::loose))myquality=0;
      if(track->quality(reco::TrackBase::tight))myquality=1;
      if(track->quality(reco::TrackBase::highPurity))myquality=2;
      //if(track->quality(reco::TrackBase::confirmed))myquality=3;
      // if(track->quality(reco::TrackBase::goodIterative))myquality=4;
      // if(track->quality(reco::TrackBase::qualitySize))myquality=5;     
      trkQuality_[nTracks_]= myquality;

      if(track->quality(reco::TrackBase::highPurity))isHighPurity_[nTracks_]=1;
      else isHighPurity_[nTracks_]=0;
      nTracks_++;


    

     
    }//end loop on tracks

  for(int d=0;d<nTracks_;++d){
    if(abs(trkQuality_[d])>5)cout<<"MYQUALITY!!! " <<trkQuality_[d] <<" at track # "<<d<<"/"<< nTracks_<<endl;
  }

 


  rootTree_->Fill();
} 


// ------------ method called once each job before begining the event loop  ------------
void LhcTrackAnalyzer::beginJob()
{
  edm::LogInfo("beginJob") << "Begin Job" << std::endl;
  // Define TTree for output
  rootFile_ = new TFile(filename_.c_str(),"recreate");
  rootTree_ = new TTree("tree","Lhc Track tree");
  
  // Track Paramters 
  rootTree_->Branch("run",&run_,"run/I");
  rootTree_->Branch("event",&event_,"event/I");
  rootTree_->Branch("goodbx",&goodbx_,"goodbx/O");
  rootTree_->Branch("goodvtx",&goodvtx_,"goodvtx/O");
  rootTree_->Branch("nTracks",&nTracks_,"nTracks/I");
  rootTree_->Branch("pt",&pt_,"pt[nTracks]/D");
  rootTree_->Branch("eta",&eta_,"eta[nTracks]/D");
  rootTree_->Branch("phi",&phi_,"phi[nTracks]/D");
  rootTree_->Branch("chi2",&chi2_,"chi2[nTracks]/D");
  rootTree_->Branch("chi2ndof",&chi2ndof_,"chi2ndof[nTracks]/D");
  rootTree_->Branch("charge",&charge_,"charge[nTracks]/I");
  rootTree_->Branch("qoverp",&qoverp_,"qoverp[nTracks]/D");
  rootTree_->Branch("dz",&dz_,"dz[nTracks]/D");
  rootTree_->Branch("dxy",&dxy_,"dxy[nTracks]/D");
  rootTree_->Branch("xPCA",&xPCA_,"xPCA[nTracks]/D");
  rootTree_->Branch("yPCA",&yPCA_,"yPCA[nTracks]/D");
  rootTree_->Branch("zPCA",&zPCA_,"zPCA[nTracks]/D");
 rootTree_->Branch("isHighPurity",&isHighPurity_,"isHighPurity[nTracks]/I");
  rootTree_->Branch("trkQuality",&trkQuality_,"trkQuality[nTracks]/I");
  rootTree_->Branch("trkAlgo",&trkAlgo_,"trkAlgo[nTracks]/I");
  rootTree_->Branch("nValidHits",&validhits_,"nValidHits[nTracks][7]/I");
 
  
}

// ------------ method called once each job just after ending the event loop  ------------
void LhcTrackAnalyzer::endJob() 
{
   if ( rootFile_ ) {
     rootFile_->Write();
     rootFile_->Close();
   }

   

}

void LhcTrackAnalyzer::SetVarToZero() {
  run_=-1;
  event_=-99;
  nTracks_ = 0;
  for ( int i=0; i<nMaxtracks_; ++i ) {
    pt_[i]        = 0;
    eta_[i]       = 0;
    phi_[i]       = 0;
    chi2_[i]      = 0;
    chi2ndof_[i]  = 0;
    charge_[i]    = 0;
    qoverp_[i]    = 0;
    dz_[i]        = 0;
    dxy_[i]       = 0;
    xPCA_[i]        = 0;
    yPCA_[i]        = 0;
    zPCA_[i]        = 0;
    trkQuality_[i]  = 0;
    trkAlgo_[i]     = -1;
    isHighPurity_[i]=-3;
    for(int j=0;j<7;j++){
      validhits_[nTracks_][j]=-1*j;
    }
  }


}

//define this as a plug-in
DEFINE_FWK_MODULE(LhcTrackAnalyzer);