RecoMuon/src/HTrack.cc

0001 #include "Validation/RecoMuon/src/HTrack.h"
0002 #include "Validation/RecoMuon/src/Histograms.h"
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0005 #include "SimDataFormats/Track/interface/SimTrack.h"
0006
0007 //#include "DQMServices/Core/interface/DQMStore.h"
0008 #include "FWCore/ServiceRegistry/interface/Service.h"
0009
0010 #include "TFile.h"
0011 #include "TDirectory.h"
0012
0013 using namespace std;
0014
0015 HTrack::HTrack(DQMStore::IBooker& ibooker, string dirName_, string name, string whereIs)
0016     : theName(name.c_str()), where(whereIs.c_str()) {
0017   ibooker.cd();
0018   std::string dirName = dirName_;
0019   dirName += "/";
0020   dirName += name;
0021   dirName += "_";
0022   dirName += whereIs;
0023
0024   ibooker.setCurrentFolder(dirName);
0025   hVariables = new HTrackVariables(ibooker, dirName, name, whereIs);
0026
0027   ibooker.cd();
0028   string resName = dirName;
0029   resName += "/Resolution";
0030   hResolution = new HResolution(ibooker, resName, name + "_Res", whereIs);
0031   ibooker.cd();
0032   ibooker.setCurrentFolder(dirName);
0033   hTDRResolution = new HResolution(ibooker, resName, name + "_TDRRes", whereIs);
0034
0035   ibooker.cd();
0036   ibooker.setCurrentFolder(dirName);
0037   string pullName = dirName;
0038   pullName += "/Pull";
0039   hPull = new HResolution(ibooker, pullName, name + "_Pull", whereIs);
0040   hTDRPull = new HResolution(ibooker, pullName, name + "_TDRPull", whereIs);
0041
0042   doSubHisto = false;
0043
0044   if (doSubHisto) {
0045     ibooker.cd();
0046     ibooker.setCurrentFolder(dirName);
0047     string subName = dirName;
0048     subName += "/subHistos";
0049     // [5-10] GeV range
0050     hResolution_5_10 = new HResolution(ibooker, subName, name + "_Res_Pt_5_10", whereIs);
0051     hTDRResolution_5_10 = new HResolution(ibooker, subName, name + "_TDRRes_Pt_5_10", whereIs);
0052     hPull_5_10 = new HResolution(ibooker, subName, name + "_Pull_Pt_5_10", whereIs);
0053
0054     hResolution_10_40 = new HResolution(ibooker, subName, name + "_Res_Pt_10_40", whereIs);
0055     hTDRResolution_10_40 = new HResolution(ibooker, subName, name + "_TDRRes_Pt_10_40", whereIs);
0056     hPull_10_40 = new HResolution(ibooker, subName, name + "_Pull_Pt_10_40", whereIs);
0057
0058     hResolution_40_70 = new HResolution(ibooker, subName, name + "_Res_Pt_40_70", whereIs);
0059     hTDRResolution_40_70 = new HResolution(ibooker, subName, name + "_TDRRes_Pt_40_70", whereIs);
0060     hPull_40_70 = new HResolution(ibooker, subName, name + "_Pull_Pt_40_70", whereIs);
0061
0062     hResolution_70_100 = new HResolution(ibooker, subName, name + "_Res_Pt_70_100", whereIs);
0063     hTDRResolution_70_100 = new HResolution(ibooker, subName, name + "_TDRRes_Pt_70_100", whereIs);
0064     hPull_70_100 = new HResolution(ibooker, subName, name + "_Pull_Pt_70_100", whereIs);
0065
0066     hResolution_08 = new HResolution(ibooker, subName, name + "_Res_Eta_08", whereIs);
0067     hTDRResolution_08 = new HResolution(ibooker, subName, name + "_TDRRes_Eta_08", whereIs);
0068     hPull_08 = new HResolution(ibooker, subName, name + "_Pull_Eta_08", whereIs);
0069
0070     hResolution_08_12 = new HResolution(ibooker, subName, name + "_Res_Eta_08_12", whereIs);
0071     hTDRResolution_08_12 = new HResolution(ibooker, subName, name + "_TDRRes_Eta_08_12", whereIs);
0072     hPull_08_12 = new HResolution(ibooker, subName, name + "_Pull_Eta_08_12", whereIs);
0073
0074     hResolution_12_21 = new HResolution(ibooker, subName, name + "_Res_Eta_12_21", whereIs);
0075     hTDRResolution_12_21 = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_21", whereIs);
0076     hPull_12_21 = new HResolution(ibooker, subName, name + "_Pull_Eta_12_21", whereIs);
0077
0078     hResolution_12_24 = new HResolution(ibooker, subName, name + "_Res_Eta_12_24", whereIs);
0079     hTDRResolution_12_24 = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_24", whereIs);
0080     hPull_12_24 = new HResolution(ibooker, subName, name + "_Pull_Eta_12_24", whereIs);
0081
0082     hResolution_12_21_plus = new HResolution(ibooker, subName, name + "_Res_Eta_12_21_plus", whereIs);
0083     hTDRResolution_12_21_plus = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_21_plus", whereIs);
0084     hPull_12_21_plus = new HResolution(ibooker, subName, name + "_Pull_Eta_12_21_plus", whereIs);
0085
0086     hResolution_12_24_plus = new HResolution(ibooker, subName, name + "_Res_Eta_12_24_plus", whereIs);
0087     hTDRResolution_12_24_plus = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_24_plus", whereIs);
0088     hPull_12_24_plus = new HResolution(ibooker, subName, name + "_Pull_Eta_12_24_plus", whereIs);
0089
0090     hResolution_12_21_minus = new HResolution(ibooker, subName, name + "_Res_Eta_12_21_minus", whereIs);
0091     hTDRResolution_12_21_minus = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_21_minus", whereIs);
0092     hPull_12_21_minus = new HResolution(ibooker, subName, name + "_Pull_Eta_12_21_minus", whereIs);
0093
0094     hResolution_12_24_minus = new HResolution(ibooker, subName, name + "_Res_Eta_12_24_minus", whereIs);
0095     hTDRResolution_12_24_minus = new HResolution(ibooker, subName, name + "_TDRRes_Eta_12_24_minus", whereIs);
0096     hPull_12_24_minus = new HResolution(ibooker, subName, name + "_Pull_Eta_12_24_minus", whereIs);
0097   }
0098 }
0099
0100 double HTrack::pull(double rec, double sim, double sigmarec) { return (rec - sim) / sigmarec; }
0101
0102 double HTrack::resolution(double rec, double sim) { return (rec - sim) / sim; }
0103
0104 void HTrack::Fill(TrajectoryStateOnSurface& tsos) { Fill(*tsos.freeState()); }
0105
0106 void HTrack::Fill(const FreeTrajectoryState& fts) {
0107   hVariables->Fill(
0108       fts.momentum().mag(), fts.momentum().perp(), fts.momentum().eta(), fts.momentum().phi(), fts.charge());
0109 }
0110
0111 void HTrack::FillDeltaR(double deltaR) { hVariables->FillDeltaR(deltaR); }
0112
0113 double HTrack::computeEfficiency(HTrackVariables* sim, DQMStore::IBooker& ibooker) {
0114   return hVariables->computeEfficiency(sim, ibooker);
0115 }
0116
0117 void HTrack::computeResolutionAndPull(TrajectoryStateOnSurface& tsos, SimTrack& simTrack) {
0118   computeResolutionAndPull(*tsos.freeState(), simTrack);
0119 }
0120
0121 void HTrack::computeResolutionAndPull(const FreeTrajectoryState& fts, SimTrack& simTrack) {
0122   // Global Resolution
0123   computeResolution(fts, simTrack, hResolution);
0124   computePull(fts, simTrack, hPull);
0125
0126   // TDR Resolution
0127   computeTDRResolution(fts, simTrack, hTDRResolution);
0128   // computeTDRPull(fts,simTracks,hTDRPull);
0129
0130   hVariables->Fill(sqrt(simTrack.momentum().Perp2()), simTrack.momentum().eta(),
0131                    simTrack.momentum().phi());  //FIXME
0132
0133   if (doSubHisto) {
0134     // [5-10] GeV range
0135     if (sqrt(simTrack.momentum().Perp2()) < 10) {
0136       computeResolution(fts, simTrack, hResolution_5_10);
0137       computeTDRResolution(fts, simTrack, hTDRResolution_5_10);
0138       computePull(fts, simTrack, hPull_5_10);
0139     }
0140
0141     // [10-40] GeV range
0142     if (sqrt(simTrack.momentum().Perp2()) >= 10 && sqrt(simTrack.momentum().Perp2()) < 40) {
0143       computeResolution(fts, simTrack, hResolution_10_40);
0144       computeTDRResolution(fts, simTrack, hTDRResolution_10_40);
0145       computePull(fts, simTrack, hPull_10_40);
0146     }
0147
0148     // [40-70] GeV range
0149     if (sqrt(simTrack.momentum().Perp2()) >= 40 && sqrt(simTrack.momentum().Perp2()) < 70) {
0150       computeResolution(fts, simTrack, hResolution_40_70);
0151       computeTDRResolution(fts, simTrack, hTDRResolution_40_70);
0152       computePull(fts, simTrack, hPull_40_70);
0153     }
0154
0155     // [70-100] GeV range
0156     if (sqrt(simTrack.momentum().Perp2()) >= 70 && sqrt(simTrack.momentum().Perp2()) < 100) {
0157       computeResolution(fts, simTrack, hResolution_70_100);
0158       computeTDRResolution(fts, simTrack, hTDRResolution_70_100);
0159       computePull(fts, simTrack, hPull_70_100);
0160     }
0161
0162     // eta range |eta|<0.8
0163     if (abs(simTrack.momentum().eta()) <= 0.8) {
0164       computeResolution(fts, simTrack, hResolution_08);
0165       computeTDRResolution(fts, simTrack, hTDRResolution_08);
0166       computePull(fts, simTrack, hPull_08);
0167     }
0168
0169     // eta range 0.8<|eta|<1.2
0170     if (abs(simTrack.momentum().eta()) > 0.8 && abs(simTrack.momentum().eta()) <= 1.2) {
0171       computeResolution(fts, simTrack, hResolution_08_12);
0172       computeTDRResolution(fts, simTrack, hTDRResolution_08_12);
0173       computePull(fts, simTrack, hPull_08_12);
0174     }
0175
0176     // eta range 1.2<|eta|<2.1
0177     if (abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.1) {
0178       computeResolution(fts, simTrack, hResolution_12_21);
0179       computeTDRResolution(fts, simTrack, hTDRResolution_12_21);
0180       computePull(fts, simTrack, hPull_12_21);
0181
0182       if (simTrack.momentum().eta() > 0) {
0183         computeResolution(fts, simTrack, hResolution_12_21_plus);
0184         computeTDRResolution(fts, simTrack, hTDRResolution_12_21_plus);
0185         computePull(fts, simTrack, hPull_12_21_plus);
0186       } else {
0187         computeResolution(fts, simTrack, hResolution_12_21_minus);
0188         computeTDRResolution(fts, simTrack, hTDRResolution_12_21_minus);
0189         computePull(fts, simTrack, hPull_12_21_minus);
0190       }
0191     }
0192
0193     // eta range 1.2<|eta|<2.4
0194     if (abs(simTrack.momentum().eta()) > 1.2 && abs(simTrack.momentum().eta()) <= 2.4) {
0195       computeResolution(fts, simTrack, hResolution_12_24);
0196       computeTDRResolution(fts, simTrack, hTDRResolution_12_24);
0197       computePull(fts, simTrack, hPull_12_24);
0198
0199       if (simTrack.momentum().eta() > 0) {
0200         computeResolution(fts, simTrack, hResolution_12_24_plus);
0201         computeTDRResolution(fts, simTrack, hTDRResolution_12_24_plus);
0202         computePull(fts, simTrack, hPull_12_24_plus);
0203       } else {
0204         computeResolution(fts, simTrack, hResolution_12_24_minus);
0205         computeTDRResolution(fts, simTrack, hTDRResolution_12_24_minus);
0206         computePull(fts, simTrack, hPull_12_24_minus);
0207       }
0208     }
0209   }
0210 }
0211
0212 void HTrack::computeResolution(const FreeTrajectoryState& fts, SimTrack& simTrack, HResolution* hReso) {
0213   hReso->Fill(simTrack.momentum().mag(),
0214               sqrt(simTrack.momentum().Perp2()),
0215               simTrack.momentum().eta(),
0216               simTrack.momentum().phi(),
0217               resolution(fts.momentum().mag(), simTrack.momentum().mag()),
0218               resolution(fts.momentum().perp(), sqrt(simTrack.momentum().Perp2())),
0219               fts.momentum().eta() - simTrack.momentum().eta(),
0220               fts.momentum().phi() - simTrack.momentum().phi(),
0221               fts.charge() + simTrack.type() / abs(simTrack.type()));  // FIXME
0222 }
0223
0224 void HTrack::computeTDRResolution(const FreeTrajectoryState& fts, SimTrack& simTrack, HResolution* hReso) {
0225   int simCharge = -simTrack.type() / abs(simTrack.type());
0226
0227   double invSimP = (simTrack.momentum().mag() == 0) ? 0 : simTrack.momentum().mag();
0228   double signedInverseRecMom = (simTrack.momentum().mag() == 0) ? 0 : fts.signedInverseMomentum();
0229
0230   hReso->Fill(simTrack.momentum().mag(),
0231               sqrt(simTrack.momentum().Perp2()),
0232               simTrack.momentum().eta(),
0233               simTrack.momentum().phi(),
0234               resolution(signedInverseRecMom, simCharge * invSimP),
0235               resolution(fts.charge() / fts.momentum().perp(), simCharge / sqrt(simTrack.momentum().Perp2())));
0236 }
0237
0238 void HTrack::computePull(const FreeTrajectoryState& fts, SimTrack& simTrack, HResolution* hReso) {
0239   // x,y,z, px,py,pz
0240   AlgebraicSymMatrix66 const errors = fts.cartesianError().matrix();
0241
0242   double partialPterror = errors[3][3] * pow(fts.momentum().x(), 2) + errors[4][4] * pow(fts.momentum().y(), 2);
0243
0244   // sqrt( (px*spx)^2 + (py*spy)^2 ) / pt
0245   double pterror = sqrt(partialPterror) / fts.momentum().perp();
0246
0247   // sqrt( (px*spx)^2 + (py*spy)^2 + (pz*spz)^2 ) / p
0248   double perror = sqrt(partialPterror + errors[5][5] * pow(fts.momentum().z(), 2)) / fts.momentum().mag();
0249
0250   double phierror = sqrt(fts.curvilinearError().matrix()[2][2]);
0251
0252   double etaerror = sqrt(fts.curvilinearError().matrix()[1][1]) * abs(sin(fts.momentum().theta()));
0253
0254   hReso->Fill(simTrack.momentum().mag(),
0255               sqrt(simTrack.momentum().Perp2()),
0256               simTrack.momentum().eta(),
0257               simTrack.momentum().phi(),
0258               pull(fts.momentum().mag(), simTrack.momentum().mag(), perror),
0259               pull(fts.momentum().perp(), sqrt(simTrack.momentum().Perp2()), pterror),
0260               pull(fts.momentum().eta(), simTrack.momentum().eta(), etaerror),
0261               pull(fts.momentum().phi(), simTrack.momentum().phi(), phierror),
0262               pull(fts.charge(), -simTrack.type() / abs(simTrack.type()), 1.));  // FIXME
0263 }