CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/Validation/RecoMuon/src/HTrack.cc

Go to the documentation of this file.
00001 #include "Validation/RecoMuon/src/HTrack.h" 
00002 #include "Validation/RecoMuon/src/Histograms.h" 
00003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00004 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00005 #include "SimDataFormats/Track/interface/SimTrack.h"
00006 
00007 #include "DQMServices/Core/interface/DQMStore.h"
00008 #include "DQMServices/Core/interface/MonitorElement.h"
00009 #include "FWCore/ServiceRegistry/interface/Service.h"
00010 
00011 #include "TFile.h"
00012 #include "TDirectory.h"
00013 
00014 
00015 using namespace std;
00016 
00017 HTrack::HTrack(string dirName_,string name,string whereIs):
00018   theName(name.c_str()),where(whereIs.c_str()){
00019 
00020   dbe_ = edm::Service<DQMStore>().operator->();
00021   dbe_->cd();
00022   std::string dirName=dirName_;
00023   dirName+="/";
00024   dirName+=name.c_str();
00025   dirName+="_";
00026   dirName+=whereIs.c_str();
00027 
00028   dbe_->cd();
00029   dbe_->setCurrentFolder(dirName.c_str());
00030   
00031   hVariables = new HTrackVariables(dirName.c_str(),name,whereIs);
00032 
00033 
00034   dbe_->cd();
00035   dbe_->setCurrentFolder(dirName.c_str());
00036   string resName = dirName;
00037   resName+="/Resolution";
00038   hResolution = new HResolution(resName.c_str(),name+"_Res",whereIs);
00039   dbe_->cd();
00040   dbe_->setCurrentFolder(dirName.c_str());
00041   hTDRResolution = new HResolution(resName.c_str(),name+"_TDRRes",whereIs);
00042 
00043   dbe_->cd();
00044   dbe_->setCurrentFolder(dirName.c_str());
00045   string pullName = dirName;
00046   pullName+="/Pull";
00047   hPull = new HResolution(pullName.c_str(),name+"_Pull",whereIs);
00048   hTDRPull = new HResolution(pullName.c_str(),name+"_TDRPull",whereIs);
00049   
00050   
00051   doSubHisto = false;
00052 
00053   if(doSubHisto){
00054     dbe_->cd();
00055     dbe_->setCurrentFolder(dirName.c_str());
00056     string subName = dirName;
00057     subName+="/subHistos";
00058     // [5-10] GeV range
00059     hResolution_5_10 = new HResolution(subName.c_str(),name+"_Res_Pt_5_10",whereIs);
00060     hTDRResolution_5_10 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_5_10",whereIs);
00061     hPull_5_10 = new HResolution(subName.c_str(),name+"_Pull_Pt_5_10",whereIs);
00062     
00063     
00064     hResolution_10_40 = new HResolution(subName.c_str(),name+"_Res_Pt_10_40",whereIs);
00065     hTDRResolution_10_40 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_10_40",whereIs);
00066     hPull_10_40 = new HResolution(subName.c_str(),name+"_Pull_Pt_10_40",whereIs);
00067     
00068     
00069     hResolution_40_70 = new HResolution(subName.c_str(),name+"_Res_Pt_40_70",whereIs);
00070     hTDRResolution_40_70 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_40_70",whereIs);
00071     hPull_40_70 = new HResolution(subName.c_str(),name+"_Pull_Pt_40_70",whereIs);
00072     
00073     
00074     hResolution_70_100 = new HResolution(subName.c_str(),name+"_Res_Pt_70_100",whereIs);
00075     hTDRResolution_70_100 = new HResolution(subName.c_str(),name+"_TDRRes_Pt_70_100",whereIs);
00076     hPull_70_100 = new HResolution(subName.c_str(),name+"_Pull_Pt_70_100",whereIs);
00077     
00078     
00079     hResolution_08 = new HResolution(subName.c_str(),name+"_Res_Eta_08",whereIs);
00080     hTDRResolution_08 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_08",whereIs);
00081     hPull_08 = new HResolution(subName.c_str(),name+"_Pull_Eta_08",whereIs);
00082     
00083         
00084     hResolution_08_12 = new HResolution(subName.c_str(),name+"_Res_Eta_08_12",whereIs);
00085     hTDRResolution_08_12 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_08_12",whereIs);
00086     hPull_08_12 = new HResolution(subName.c_str(),name+"_Pull_Eta_08_12",whereIs);
00087     
00088     
00089     hResolution_12_21 = new HResolution(subName.c_str(),name+"_Res_Eta_12_21",whereIs);
00090     hTDRResolution_12_21 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21",whereIs);
00091     hPull_12_21 = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21",whereIs);
00092     
00093     
00094     hResolution_12_24 = new HResolution(subName.c_str(),name+"_Res_Eta_12_24",whereIs);
00095     hTDRResolution_12_24 = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24",whereIs);
00096     hPull_12_24 = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24",whereIs);
00097 
00098     
00099     hResolution_12_21_plus = new HResolution(subName.c_str(),name+"_Res_Eta_12_21_plus",whereIs);
00100     hTDRResolution_12_21_plus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21_plus",whereIs);
00101     hPull_12_21_plus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21_plus",whereIs);
00102     
00103     
00104     hResolution_12_24_plus = new HResolution(subName.c_str(),name+"_Res_Eta_12_24_plus",whereIs);
00105     hTDRResolution_12_24_plus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24_plus",whereIs);
00106     hPull_12_24_plus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24_plus",whereIs);
00107     
00108     
00109     hResolution_12_21_minus = new HResolution(subName.c_str(),name+"_Res_Eta_12_21_minus",whereIs);
00110     hTDRResolution_12_21_minus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_21_minus",whereIs);
00111     hPull_12_21_minus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_21_minus",whereIs);
00112     
00113     
00114     hResolution_12_24_minus = new HResolution(subName.c_str(),name+"_Res_Eta_12_24_minus",whereIs);
00115     hTDRResolution_12_24_minus = new HResolution(subName.c_str(),name+"_TDRRes_Eta_12_24_minus",whereIs);
00116     hPull_12_24_minus = new HResolution(subName.c_str(),name+"_Pull_Eta_12_24_minus",whereIs);
00117   }
00118 }
00119 
00120 double HTrack::pull(double rec,double sim, double sigmarec){
00121   return (rec-sim)/sigmarec;
00122 }
00123 
00124 double HTrack::resolution(double rec,double sim){
00125   return (rec-sim)/sim;
00126 }
00127 
00128 
00129 void HTrack::Fill(TrajectoryStateOnSurface &tsos){ 
00130   Fill(*tsos.freeState());
00131 }
00132 
00133 void HTrack::Fill(FreeTrajectoryState &fts){
00134   
00135   hVariables->Fill(fts.momentum().mag(),
00136                    fts.momentum().perp(),
00137                    fts.momentum().eta(),
00138                    fts.momentum().phi(),
00139                    fts.charge());  
00140 }
00141 
00142 void HTrack::FillDeltaR(double deltaR){
00143   hVariables->FillDeltaR(deltaR);
00144 }
00145 
00146 
00147 double HTrack::computeEfficiency(HTrackVariables *sim){
00148   return hVariables->computeEfficiency(sim);
00149 }
00150 
00151 
00152 void HTrack::computeResolutionAndPull(TrajectoryStateOnSurface& tsos, SimTrack& simTrack){
00153   computeResolutionAndPull(*tsos.freeState(),simTrack);
00154 }
00155 
00156 
00157 void HTrack::computeResolutionAndPull(FreeTrajectoryState& fts, SimTrack& simTrack){
00158   
00159 
00160   // Global Resolution
00161   computeResolution(fts,simTrack,hResolution); 
00162   computePull(fts,simTrack,hPull); 
00163 
00164   // TDR Resolution
00165   computeTDRResolution(fts,simTrack,hTDRResolution);
00166   // computeTDRPull(fts,simTracks,hTDRPull); 
00167 
00168 
00169   hVariables->Fill(sqrt(simTrack.momentum().Perp2()),
00170                    simTrack.momentum().eta(),
00171                    simTrack.momentum().phi()); //FIXME
00172 
00173 
00174   if(doSubHisto){
00175     // [5-10] GeV range
00176     if(sqrt(simTrack.momentum().Perp2()) <10 ){  
00177       computeResolution(fts,simTrack,hResolution_5_10);
00178       computeTDRResolution(fts,simTrack,hTDRResolution_5_10);
00179       computePull(fts,simTrack,hPull_5_10);
00180     }
00181 
00182     // [10-40] GeV range
00183     if(sqrt(simTrack.momentum().Perp2()) >=10 && sqrt(simTrack.momentum().Perp2()) < 40){
00184       computeResolution(fts,simTrack,hResolution_10_40); 
00185       computeTDRResolution(fts,simTrack,  hTDRResolution_10_40);
00186       computePull(fts,simTrack,hPull_10_40);
00187     }
00188   
00189     // [40-70] GeV range
00190     if(sqrt(simTrack.momentum().Perp2()) >=40 && sqrt(simTrack.momentum().Perp2()) < 70){
00191       computeResolution(fts,simTrack,hResolution_40_70); 
00192       computeTDRResolution(fts,simTrack,hTDRResolution_40_70);
00193       computePull(fts,simTrack,hPull_40_70); 
00194     }
00195 
00196     // [70-100] GeV range
00197     if(sqrt(simTrack.momentum().Perp2()) >= 70 && sqrt(simTrack.momentum().Perp2()) < 100){            
00198       computeResolution(fts,simTrack,hResolution_70_100); 
00199       computeTDRResolution(fts,simTrack,hTDRResolution_70_100);
00200       computePull(fts,simTrack,hPull_70_100); 
00201     }
00202   
00203     // eta range |eta|<0.8
00204     if(abs(simTrack.momentum().eta()) <= 0.8){
00205       computeResolution(fts,simTrack,hResolution_08); 
00206       computeTDRResolution(fts,simTrack,  hTDRResolution_08);
00207       computePull(fts,simTrack,hPull_08);
00208     }
00209 
00210     // eta range 0.8<|eta|<1.2
00211     if( abs(simTrack.momentum().eta()) > 0.8 && abs(simTrack.momentum().eta()) <= 1.2 ){
00212       computeResolution(fts,simTrack,hResolution_08_12); 
00213       computeTDRResolution(fts,simTrack,  hTDRResolution_08_12);
00214       computePull(fts,simTrack,hPull_08_12);
00215     }
00216 
00217     // eta range 1.2<|eta|<2.1
00218     if( abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.1 ){
00219       computeResolution(fts,simTrack,hResolution_12_21); 
00220       computeTDRResolution(fts,simTrack,  hTDRResolution_12_21);
00221       computePull(fts,simTrack,hPull_12_21);
00222     
00223       if(simTrack.momentum().eta() > 0){
00224         computeResolution(fts,simTrack,hResolution_12_21_plus); 
00225         computeTDRResolution(fts,simTrack,  hTDRResolution_12_21_plus);
00226         computePull(fts,simTrack,hPull_12_21_plus);
00227       }
00228       else{
00229         computeResolution(fts,simTrack,hResolution_12_21_minus); 
00230         computeTDRResolution(fts,simTrack,  hTDRResolution_12_21_minus);
00231         computePull(fts,simTrack,hPull_12_21_minus);
00232       }  
00233     }
00234 
00235     // eta range 1.2<|eta|<2.4
00236     if( abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.4 ){
00237       computeResolution(fts,simTrack,hResolution_12_24); 
00238       computeTDRResolution(fts,simTrack,  hTDRResolution_12_24);
00239       computePull(fts,simTrack,hPull_12_24);
00240 
00241       if(simTrack.momentum().eta() > 0){
00242         computeResolution(fts,simTrack,hResolution_12_24_plus); 
00243         computeTDRResolution(fts,simTrack,  hTDRResolution_12_24_plus);
00244         computePull(fts,simTrack,hPull_12_24_plus);
00245       }
00246       else{
00247         computeResolution(fts,simTrack,hResolution_12_24_minus); 
00248         computeTDRResolution(fts,simTrack,  hTDRResolution_12_24_minus);
00249         computePull(fts,simTrack,hPull_12_24_minus);
00250       }
00251     }
00252   }
00253 }
00254 
00255 void HTrack::computeResolution(FreeTrajectoryState &fts,
00256                                SimTrack& simTrack,
00257                                HResolution* hReso){
00258 
00259 
00260   hReso->Fill(simTrack.momentum().mag(),
00261               sqrt(simTrack.momentum().Perp2()),
00262               simTrack.momentum().eta(),
00263               simTrack.momentum().phi(),
00264               resolution(fts.momentum().mag(),simTrack.momentum().mag()),
00265               resolution(fts.momentum().perp(),sqrt(simTrack.momentum().Perp2())),
00266               fts.momentum().eta()-simTrack.momentum().eta(),
00267               fts.momentum().phi()-simTrack.momentum().phi(),
00268               fts.charge()+simTrack.type()/abs(simTrack.type())); // FIXME
00269 }
00270 
00271 void HTrack::computeTDRResolution(FreeTrajectoryState &fts,
00272                                   SimTrack& simTrack,
00273                                   HResolution* hReso){
00274 
00275   int simCharge = -simTrack.type()/ abs(simTrack.type());
00276 
00277   double invSimP =  (simTrack.momentum().mag() == 0 ) ? 0 : simTrack.momentum().mag();
00278   double signedInverseRecMom = (simTrack.momentum().mag() == 0 ) ? 0 : fts.signedInverseMomentum();
00279 
00280   hReso->Fill(simTrack.momentum().mag(),
00281               sqrt(simTrack.momentum().Perp2()),
00282               simTrack.momentum().eta(),
00283               simTrack.momentum().phi(),
00284               resolution(signedInverseRecMom , simCharge*invSimP ),
00285               resolution(fts.charge()/fts.momentum().perp(),simCharge/sqrt(simTrack.momentum().Perp2())));
00286 }
00287 
00288 void HTrack::computePull(FreeTrajectoryState &fts,
00289                          SimTrack& simTrack,
00290                          HResolution* hReso){
00291 
00292   // x,y,z, px,py,pz
00293   AlgebraicSymMatrix66 const & errors = fts.cartesianError().matrix();
00294   
00295   double partialPterror = errors[3][3]*pow(fts.momentum().x(),2) + errors[4][4]*pow(fts.momentum().y(),2);
00296   
00297   // sqrt( (px*spx)^2 + (py*spy)^2 ) / pt
00298   double pterror = sqrt(partialPterror)/fts.momentum().perp();
00299   
00300   // sqrt( (px*spx)^2 + (py*spy)^2 + (pz*spz)^2 ) / p
00301   double perror = sqrt(partialPterror+errors[5][5]*pow(fts.momentum().z(),2))/fts.momentum().mag();
00302 
00303   double phierror = sqrt(fts.curvilinearError().matrix()[2][2]);
00304 
00305   double etaerror = sqrt(fts.curvilinearError().matrix()[1][1])*abs(sin(fts.momentum().theta()));
00306 
00307 
00308   hReso->Fill(simTrack.momentum().mag(),
00309               sqrt(simTrack.momentum().Perp2()),
00310               simTrack.momentum().eta(),
00311               simTrack.momentum().phi(),
00312               pull(fts.momentum().mag(),simTrack.momentum().mag(),perror),
00313               pull(fts.momentum().perp(),sqrt(simTrack.momentum().Perp2()),pterror),
00314               pull(fts.momentum().eta(),simTrack.momentum().eta(),etaerror),
00315               pull(fts.momentum().phi(),simTrack.momentum().phi(),phierror),
00316               pull(fts.charge() , -simTrack.type()/ abs(simTrack.type()), 1.)); // FIXME
00317 }