TrackParameterAnalyzer.cc

Go to the documentation of this file.

#include "Validation/RecoVertex/interface/TrackParameterAnalyzer.h"

//system includes
#include <memory>
#include <vector>

// core framework
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Version/interface/GetReleaseVersion.h"

// Hep MC stuff from CLHEP
#include "CLHEP/Vector/LorentzVector.h"

// track
#include "DataFormats/TrackReco/interface/Track.h"

// Root
#include <TH1.h>
#include <TH2.h>
#include <TFile.h>

//
//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
TrackParameterAnalyzer::TrackParameterAnalyzer(const edm::ParameterSet& iConfig)
  : edmSimVertexContainerToken_( consumes<edm::SimVertexContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
  , edmSimTrackContainerToken_( consumes<edm::SimTrackContainer>( iConfig.getParameter<edm::InputTag>( "simG4" ) ) )
  , recoTrackCollectionToken_( consumes<reco::TrackCollection>( edm::InputTag( iConfig.getUntrackedParameter<std::string>( "recoTrackProducer" ) ) ) )
  , outputFile_( iConfig.getUntrackedParameter<std::string>( "outputFile" ) )
  , simUnit_( 1.0 ) //  starting from  CMSSW_1_2_x, I think
  , verbose_( iConfig.getUntrackedParameter<bool>( "verbose", false ) ) {
   //now do whatever initialization is needed
  // open output file to store histograms}
  auto tversion = edm::getReleaseVersion();
  tversion = tversion.erase(tversion.size() - 1, 1).erase(0, 1);
  outputFile_  = tversion + "_" + outputFile_;
  rootFile_ = TFile::Open(outputFile_.c_str(),"RECREATE"); 
  if ( (edm::getReleaseVersion()).find("CMSSW_1_1_",0)!=std::string::npos){
    simUnit_=0.1;  // for use in  CMSSW_1_1_1 tutorial
  }
}


TrackParameterAnalyzer::~TrackParameterAnalyzer()
{
   // do anything here that needs to be done at destruction time
   // (e.g. close files, deallocate resources etc.)
  delete rootFile_;
}



//
// member functions
//
void TrackParameterAnalyzer::beginJob(){
  std::cout << " TrackParameterAnalyzer::beginJob  conversion from sim units to rec units is " << simUnit_ << std::endl;

  rootFile_->cd();
  h1_pull0_ = new TH1F("pull0","pull q/p",100,-25.,25.);
  h1_pull1_ = new TH1F("pull1","pull lambda",100,-25.,25.);
  h1_pull2_ = new TH1F("pull2","pull phi  ",100,-25.,25.);
  h1_pull3_ = new TH1F("pull3","pull dca  ",100,-25.,25.);
  h1_pull4_ = new TH1F("pull4","pull zdca ",100,-25.,25.);

  h1_res0_ = new TH1F("res0","res q/p",100,-0.1,0.1);
  h1_res1_ = new TH1F("res1","res lambda",100,-0.1,0.1);
  h1_res2_ = new TH1F("res2","res phi  ",100,-0.1,0.1);
  h1_res3_ = new TH1F("res3","res dca  ",100,-0.1,0.1);
  h1_res4_ = new TH1F("res4","res zdca ",100,-0.1,0.1);

  h1_Beff_  = new TH1F("Beff", "Beff",2000,-10.,10.);
  h2_dvsphi_ = new TH2F("dvsphi","dvsphi",360,-M_PI,M_PI,100,-0.1,0.1);
  h1_par0_ = new TH1F("par0","q/p",100,-0.1,0.1);
  h1_par1_ = new TH1F("par1","lambda",100, -M_PI/2.,M_PI/2.);
  h1_par2_ = new TH1F("par2","phi  ",100,-M_PI,M_PI);
  h1_par3_ = new TH1F("par3","dca  ",100,-0.1,0.1);
  h1_par4_ = new TH1F("par4","zdca ",1000,-10.,10.);
}


void TrackParameterAnalyzer::endJob() {
  rootFile_->cd();
  h1_pull0_->Write();
  h1_pull1_->Write();
  h1_pull2_->Write();
  h1_pull3_->Write();
  h1_pull4_->Write();

  h1_res0_->Write();
  h1_res1_->Write();
  h1_res2_->Write();
  h1_res3_->Write();
  h1_res4_->Write();

  h1_Beff_->Write();
  h2_dvsphi_->Write();
  h1_par0_->Write();
  h1_par1_->Write();
  h1_par2_->Write();
  h1_par3_->Write();
  h1_par4_->Write();
}

// helper function
bool TrackParameterAnalyzer::match(const ParameterVector  &a, 
                   const ParameterVector  &b){
  double dtheta=a(1)-b(1);
  double dphi  =a(2)-b(2);
  if (dphi>M_PI){ dphi-=M_2_PI; }else if(dphi<-M_PI){dphi+=M_2_PI;}
  return ( (fabs(dtheta)<0.02) && (fabs(dphi)<0.04) );
}

// ------------ method called to produce the data  ------------
void
TrackParameterAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using CLHEP::HepLorentzVector;

   const double fBfield=3.8;
  
   edm::Handle<edm::SimVertexContainer> simVtcs;
   iEvent.getByToken( edmSimVertexContainerToken_, simVtcs );
   if(verbose_){
     std::cout << "SimVertex " << simVtcs->size() << std::endl;
     for(edm::SimVertexContainer::const_iterator v=simVtcs->begin();
     v!=simVtcs->end(); ++v){
       std::cout << "simvtx "
         << std::setw(10) << std::setprecision(4)
         << v->position().x() << " "
         << v->position().y() << " "
         << v->position().z() << " "
         << v->parentIndex() << " "
         << v->noParent() << " "
         << std::endl;
     }
   }
   
   // get the simulated tracks, extract perigee parameters
   edm::Handle<edm::SimTrackContainer> simTrks;
   iEvent.getByToken( edmSimTrackContainerToken_, simTrks );
   
   if(verbose_){std::cout << "simtrks " << simTrks->size() << std::endl;}
   std::vector<ParameterVector > tsim;
   for(edm::SimTrackContainer::const_iterator t=simTrks->begin();
       t!=simTrks->end(); ++t){
     if (t->noVertex()){
       std::cout << "simtrk  has no vertex" << std::endl;
       return;
     }else{
       // get the vertex position
       HepLorentzVector v((*simVtcs)[t->vertIndex()].position().x(),
                          (*simVtcs)[t->vertIndex()].position().y(),
                          (*simVtcs)[t->vertIndex()].position().z(),
                          (*simVtcs)[t->vertIndex()].position().e());
       int pdgCode=t->type();

       if( pdgCode==-99 ){
     // such entries cause crashes, no idea what they are
     std::cout << "funny particle skipped  , code="  << pdgCode << std::endl;
       }else{
     double Q=0;
     if ((pdgCode==11)||(pdgCode==13)||(pdgCode==15)||(pdgCode==-211)||(pdgCode==-2212)||(pdgCode==321)){Q=-1;}
     else if((pdgCode==-11)||(pdgCode==-13)||(pdgCode==-15)||(pdgCode==211)||(pdgCode==2212)||(pdgCode==321)){Q=1;}
     else {
       std::cout << pdgCode << " " <<std::endl;
     }
     HepLorentzVector p(t->momentum().x(),t->momentum().y(),t->momentum().z(),t->momentum().e());
     if(verbose_){
       std::cout << "simtrk "
               << " gen=" << std::setw( 4) << t->genpartIndex()
               << " vtx=" << std::setw( 4) << t->vertIndex() 
               << " pdg=" << std::setw( 5) << t->type() 
               << " Q="   << std::setw( 3) << Q 
               << " pt="  << std::setw(6) << p.perp()
               << " vx="  << std::setw(6) << v.x() 
               << " vy="  << std::setw(6) << v.y() 
               << " vz="  << std::setw(6) << v.z()
               << std::endl;
     }
     if ( (Q != 0) && (p.perp()>0.1) ){
       double x0=v.x()*simUnit_;
       double y0=v.y()*simUnit_;
       double z0=v.z()*simUnit_;
       double kappa=-Q*0.002998*fBfield/p.perp();
       double D0=x0*sin(p.phi())-y0*cos(p.phi())-0.5*kappa*(x0*x0+y0*y0);
       double q=sqrt(1.-2.*kappa*D0);
       double s0=(x0*cos(p.phi())+y0*sin(p.phi()))/q;
       double s1;
       if (fabs(kappa*s0)>0.001){
         s1=asin(kappa*s0)/kappa;
       }else{
         double ks02=(kappa*s0)*(kappa*s0);
         s1=s0*(1.+ks02/6.+3./40.*ks02*ks02+5./112.*pow(ks02,3));
       }
       ParameterVector par;
       par[reco::TrackBase::i_qoverp] = Q/sqrt(p.perp2()+p.pz()*p.pz());
       par[reco::TrackBase::i_lambda] = M_PI/2.-p.theta();
       par[reco::TrackBase::i_phi] = p.phi()-asin(kappa*s0);
       par[reco::TrackBase::i_dxy] = 2.*D0/(1.+q);
       par[reco::TrackBase::i_dsz] = z0*sin(p.theta())-s1*cos(p.theta());
       tsim.push_back(par);
     }
       }
     }// has vertex
   }//for loop
   
   // simtrack parameters are in now tsim
   // loop over tracks and try to match them to simulated tracks


   edm::Handle<reco::TrackCollection> recTracks;
   iEvent.getByToken( recoTrackCollectionToken_, recTracks );

   for(reco::TrackCollection::const_iterator t=recTracks->begin();
       t!=recTracks->end(); ++t){
     reco::TrackBase::ParameterVector  p = t->parameters();
     reco::TrackBase::CovarianceMatrix c = t->covariance();
     if(verbose_){
       std::cout << "reco pars= " << p << std::endl;
     }
     for(std::vector<ParameterVector>::const_iterator s=tsim.begin();
     s!=tsim.end(); ++s){
       if (match(*s,p)){
     h1_pull0_->Fill((p(0)-(*s)(0))/sqrt(c(0,0)));
     h1_pull1_->Fill((p(1)-(*s)(1))/sqrt(c(1,1)));
     h1_pull2_->Fill((p(2)-(*s)(2))/sqrt(c(2,2)));
     h1_pull3_->Fill((p(3)-(*s)(3))/sqrt(c(3,3)));
     h1_pull4_->Fill((p(4)-(*s)(4))/sqrt(c(4,4)));

     h1_res0_->Fill(p(0)-(*s)(0));
     h1_res1_->Fill(p(1)-(*s)(1));
     h1_res2_->Fill(p(2)-(*s)(2));
     h1_res3_->Fill(p(3)-(*s)(3));
     h1_res4_->Fill(p(4)-(*s)(4));

     h1_Beff_->Fill(p(0)/(*s)(0)*fBfield);
     h2_dvsphi_->Fill(p(2),p(3));
     h1_par0_->Fill(p(0));
     h1_par1_->Fill(p(1));
     h1_par2_->Fill(p(2));
     h1_par3_->Fill(p(3));
     h1_par4_->Fill(p(4));
      }
     }
   }

}