HTrack.cc

Go to the documentation of this file.

#include "Validation/RecoMuon/src/HTrack.h" 
#include "Validation/RecoMuon/src/Histograms.h" 
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
#include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
#include "SimDataFormats/Track/interface/SimTrack.h"

#include "DQMServices/Core/interface/DQMStore.h"
#include "DQMServices/Core/interface/MonitorElement.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "TFile.h"
#include "TDirectory.h"


using namespace std;

HTrack::HTrack(string dirName_,string name,string whereIs):
  theName(name.c_str()),where(whereIs.c_str()){

  dbe_ = edm::Service<DQMStore>().operator->();
  dbe_->cd();
  std::string dirName=dirName_;
  dirName+="/";
  dirName+=name.c_str();
  dirName+="_";
  dirName+=whereIs.c_str();

  dbe_->cd();
  dbe_->setCurrentFolder(dirName.c_str());
  
  hVariables = new HTrackVariables(dirName.c_str(),name,whereIs);


  dbe_->cd();
  dbe_->setCurrentFolder(dirName.c_str());
  string resName = dirName;
  resName+="/Resolution";
  hResolution = new HResolution(resName.c_str(),name+"_Res",whereIs);
  dbe_->cd();
  dbe_->setCurrentFolder(dirName.c_str());
  hTDRResolution = new HResolution(resName.c_str(),name+"_TDRRes",whereIs);

  dbe_->cd();
  dbe_->setCurrentFolder(dirName.c_str());
  string pullName = dirName;
  pullName+="/Pull";
  hPull = new HResolution(pullName.c_str(),name+"_Pull",whereIs);
  hTDRPull = new HResolution(pullName.c_str(),name+"_TDRPull",whereIs);
  
  
  doSubHisto = false;

  if(doSubHisto){
    dbe_->cd();
    dbe_->setCurrentFolder(dirName.c_str());
    string subName = dirName;
    subName+="/subHistos";
    // [5-10] GeV range
    hResolution_5_10 = new HResolution(subName.c_str(),name+"_Res_Pt_5_10",whereIs);
    hTDRResolution_5_10 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_5_10",whereIs);
    hPull_5_10 = new HResolution(subName.c_str(),name+"_Pull_Pt_5_10",whereIs);
    
    
    hResolution_10_40 = new HResolution(subName.c_str(),name+"_Res_Pt_10_40",whereIs);
    hTDRResolution_10_40 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_10_40",whereIs);
    hPull_10_40 = new HResolution(subName.c_str(),name+"_Pull_Pt_10_40",whereIs);
    
    
    hResolution_40_70 = new HResolution(subName.c_str(),name+"_Res_Pt_40_70",whereIs);
    hTDRResolution_40_70 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_40_70",whereIs);
    hPull_40_70 = new HResolution(subName.c_str(),name+"_Pull_Pt_40_70",whereIs);
    
    
    hResolution_70_100 = new HResolution(subName.c_str(),name+"_Res_Pt_70_100",whereIs);
    hTDRResolution_70_100 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_70_100",whereIs);
    hPull_70_100 = new HResolution(subName.c_str(),name+"_Pull_Pt_70_100",whereIs);
    
    
    hResolution_08 = new HResolution(subName.c_str(),name+"_Res_Eta_08",whereIs);
    hTDRResolution_08 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_08",whereIs);
    hPull_08 = new HResolution(subName.c_str(),name+"_Pull_Eta_08",whereIs);
    
    
    hResolution_08_12 = new HResolution(subName.c_str(),name+"_Res_Eta_08_12",whereIs);
    hTDRResolution_08_12 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_08_12",whereIs);
    hPull_08_12 = new HResolution(subName.c_str(),name+"_Pull_Eta_08_12",whereIs);
    
    
    hResolution_12_21 = new HResolution(subName.c_str(),name+"_Res_Eta_12_21",whereIs);
    hTDRResolution_12_21 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21",whereIs);
    hPull_12_21 = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21",whereIs);
    
    
    hResolution_12_24 = new HResolution(subName.c_str(),name+"_Res_Eta_12_24",whereIs);
    hTDRResolution_12_24 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24",whereIs);
    hPull_12_24 = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24",whereIs);

    
    hResolution_12_21_plus = new HResolution(subName.c_str(),name+"_Res_Eta_12_21_plus",whereIs);
    hTDRResolution_12_21_plus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21_plus",whereIs);
    hPull_12_21_plus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21_plus",whereIs);
    
    
    hResolution_12_24_plus = new HResolution(subName.c_str(),name+"_Res_Eta_12_24_plus",whereIs);
    hTDRResolution_12_24_plus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24_plus",whereIs);
    hPull_12_24_plus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24_plus",whereIs);
    
    
    hResolution_12_21_minus = new HResolution(subName.c_str(),name+"_Res_Eta_12_21_minus",whereIs);
    hTDRResolution_12_21_minus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21_minus",whereIs);
    hPull_12_21_minus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21_minus",whereIs);
    
    
    hResolution_12_24_minus = new HResolution(subName.c_str(),name+"_Res_Eta_12_24_minus",whereIs);
    hTDRResolution_12_24_minus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24_minus",whereIs);
    hPull_12_24_minus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24_minus",whereIs);
  }
}

double HTrack::pull(double rec,double sim, double sigmarec){
  return (rec-sim)/sigmarec;
}

double HTrack::resolution(double rec,double sim){
  return (rec-sim)/sim;
}


void HTrack::Fill(TrajectoryStateOnSurface &tsos){ 
  Fill(*tsos.freeState());
}

void HTrack::Fill(FreeTrajectoryState &fts){
  
  hVariables->Fill(fts.momentum().mag(),
           fts.momentum().perp(),
           fts.momentum().eta(),
           fts.momentum().phi(),
           fts.charge());  
}

void HTrack::FillDeltaR(double deltaR){
  hVariables->FillDeltaR(deltaR);
}


double HTrack::computeEfficiency(HTrackVariables *sim){
  return hVariables->computeEfficiency(sim);
}


void HTrack::computeResolutionAndPull(TrajectoryStateOnSurface& tsos, SimTrack& simTrack){
  computeResolutionAndPull(*tsos.freeState(),simTrack);
}


void HTrack::computeResolutionAndPull(FreeTrajectoryState& fts, SimTrack& simTrack){
  

  // Global Resolution
  computeResolution(fts,simTrack,hResolution); 
  computePull(fts,simTrack,hPull); 

  // TDR Resolution
  computeTDRResolution(fts,simTrack,hTDRResolution);
  // computeTDRPull(fts,simTracks,hTDRPull); 


  hVariables->Fill(sqrt(simTrack.momentum().Perp2()),
           simTrack.momentum().eta(),
           simTrack.momentum().phi()); //FIXME


  if(doSubHisto){
    // [5-10] GeV range
    if(sqrt(simTrack.momentum().Perp2()) <10 ){  
      computeResolution(fts,simTrack,hResolution_5_10);
      computeTDRResolution(fts,simTrack,hTDRResolution_5_10);
      computePull(fts,simTrack,hPull_5_10);
    }

    // [10-40] GeV range
    if(sqrt(simTrack.momentum().Perp2()) >=10 && sqrt(simTrack.momentum().Perp2()) < 40){
      computeResolution(fts,simTrack,hResolution_10_40); 
      computeTDRResolution(fts,simTrack,  hTDRResolution_10_40);
      computePull(fts,simTrack,hPull_10_40);
    }
  
    // [40-70] GeV range
    if(sqrt(simTrack.momentum().Perp2()) >=40 && sqrt(simTrack.momentum().Perp2()) < 70){
      computeResolution(fts,simTrack,hResolution_40_70); 
      computeTDRResolution(fts,simTrack,hTDRResolution_40_70);
      computePull(fts,simTrack,hPull_40_70); 
    }

    // [70-100] GeV range
    if(sqrt(simTrack.momentum().Perp2()) >= 70 && sqrt(simTrack.momentum().Perp2()) < 100){        
      computeResolution(fts,simTrack,hResolution_70_100); 
      computeTDRResolution(fts,simTrack,hTDRResolution_70_100);
      computePull(fts,simTrack,hPull_70_100); 
    }
  
    // eta range |eta|<0.8
    if(abs(simTrack.momentum().eta()) <= 0.8){
      computeResolution(fts,simTrack,hResolution_08); 
      computeTDRResolution(fts,simTrack,  hTDRResolution_08);
      computePull(fts,simTrack,hPull_08);
    }

    // eta range 0.8<|eta|<1.2
    if( abs(simTrack.momentum().eta()) > 0.8 && abs(simTrack.momentum().eta()) <= 1.2 ){
      computeResolution(fts,simTrack,hResolution_08_12); 
      computeTDRResolution(fts,simTrack,  hTDRResolution_08_12);
      computePull(fts,simTrack,hPull_08_12);
    }

    // eta range 1.2<|eta|<2.1
    if( abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.1 ){
      computeResolution(fts,simTrack,hResolution_12_21); 
      computeTDRResolution(fts,simTrack,  hTDRResolution_12_21);
      computePull(fts,simTrack,hPull_12_21);
    
      if(simTrack.momentum().eta() > 0){
    computeResolution(fts,simTrack,hResolution_12_21_plus); 
    computeTDRResolution(fts,simTrack,  hTDRResolution_12_21_plus);
    computePull(fts,simTrack,hPull_12_21_plus);
      }
      else{
    computeResolution(fts,simTrack,hResolution_12_21_minus); 
    computeTDRResolution(fts,simTrack,  hTDRResolution_12_21_minus);
    computePull(fts,simTrack,hPull_12_21_minus);
      }  
    }

    // eta range 1.2<|eta|<2.4
    if( abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.4 ){
      computeResolution(fts,simTrack,hResolution_12_24); 
      computeTDRResolution(fts,simTrack,  hTDRResolution_12_24);
      computePull(fts,simTrack,hPull_12_24);

      if(simTrack.momentum().eta() > 0){
    computeResolution(fts,simTrack,hResolution_12_24_plus); 
    computeTDRResolution(fts,simTrack,  hTDRResolution_12_24_plus);
    computePull(fts,simTrack,hPull_12_24_plus);
      }
      else{
    computeResolution(fts,simTrack,hResolution_12_24_minus); 
    computeTDRResolution(fts,simTrack,  hTDRResolution_12_24_minus);
    computePull(fts,simTrack,hPull_12_24_minus);
      }
    }
  }
}

void HTrack::computeResolution(FreeTrajectoryState &fts,
                   SimTrack& simTrack,
                   HResolution* hReso){


  hReso->Fill(simTrack.momentum().mag(),
          sqrt(simTrack.momentum().Perp2()),
          simTrack.momentum().eta(),
          simTrack.momentum().phi(),
          resolution(fts.momentum().mag(),simTrack.momentum().mag()),
          resolution(fts.momentum().perp(),sqrt(simTrack.momentum().Perp2())),
          fts.momentum().eta()-simTrack.momentum().eta(),
          fts.momentum().phi()-simTrack.momentum().phi(),
          fts.charge()+simTrack.type()/abs(simTrack.type())); // FIXME
}

void HTrack::computeTDRResolution(FreeTrajectoryState &fts,
                  SimTrack& simTrack,
                  HResolution* hReso){

  int simCharge = -simTrack.type()/ abs(simTrack.type());

  double invSimP =  (simTrack.momentum().mag() == 0 ) ? 0 : simTrack.momentum().mag();
  double signedInverseRecMom = (simTrack.momentum().mag() == 0 ) ? 0 : fts.signedInverseMomentum();

  hReso->Fill(simTrack.momentum().mag(),
          sqrt(simTrack.momentum().Perp2()),
          simTrack.momentum().eta(),
          simTrack.momentum().phi(),
          resolution(signedInverseRecMom , simCharge*invSimP ),
          resolution(fts.charge()/fts.momentum().perp(),simCharge/sqrt(simTrack.momentum().Perp2())));
}

void HTrack::computePull(FreeTrajectoryState &fts,
             SimTrack& simTrack,
             HResolution* hReso){

  // x,y,z, px,py,pz
  AlgebraicSymMatrix66 const & errors = fts.cartesianError().matrix();
  
  double partialPterror = errors[3][3]*pow(fts.momentum().x(),2) + errors[4][4]*pow(fts.momentum().y(),2);
  
  // sqrt( (px*spx)^2 + (py*spy)^2 ) / pt
  double pterror = sqrt(partialPterror)/fts.momentum().perp();
  
  // sqrt( (px*spx)^2 + (py*spy)^2 + (pz*spz)^2 ) / p
  double perror = sqrt(partialPterror+errors[5][5]*pow(fts.momentum().z(),2))/fts.momentum().mag();

  double phierror = sqrt(fts.curvilinearError().matrix()[2][2]);

  double etaerror = sqrt(fts.curvilinearError().matrix()[1][1])*abs(sin(fts.momentum().theta()));


  hReso->Fill(simTrack.momentum().mag(),
          sqrt(simTrack.momentum().Perp2()),
          simTrack.momentum().eta(),
          simTrack.momentum().phi(),
          pull(fts.momentum().mag(),simTrack.momentum().mag(),perror),
          pull(fts.momentum().perp(),sqrt(simTrack.momentum().Perp2()),pterror),
          pull(fts.momentum().eta(),simTrack.momentum().eta(),etaerror),
          pull(fts.momentum().phi(),simTrack.momentum().phi(),phierror),
          pull(fts.charge() , -simTrack.type()/ abs(simTrack.type()), 1.)); // FIXME
}