00001 #include <iostream>
00002 #include <sstream>
00003 #include <istream>
00004 #include <fstream>
00005 #include <iomanip>
00006 #include <string>
00007 #include <cmath>
00008 #include <functional>
00009 #include <stdlib.h>
00010 #include <string.h>
00011
00012 #include "HLTrigger/HLTanalyzers/interface/HLTMuon.h"
00013
00014 HLTMuon::HLTMuon() {
00015 evtCounter=0;
00016
00017
00018 _Monte=false;
00019 _Debug=false;
00020 }
00021
00022
00023 void HLTMuon::setup(const edm::ParameterSet& pSet, TTree* HltTree) {
00024
00025 edm::ParameterSet myEmParams = pSet.getParameter<edm::ParameterSet>("RunParameters") ;
00026 std::vector<std::string> parameterNames = myEmParams.getParameterNames() ;
00027
00028 for ( std::vector<std::string>::iterator iParam = parameterNames.begin();
00029 iParam != parameterNames.end(); iParam++ ){
00030 if ( (*iParam) == "Monte" ) _Monte = myEmParams.getParameter<bool>( *iParam );
00031 else if ( (*iParam) == "Debug" ) _Debug = myEmParams.getParameter<bool>( *iParam );
00032 }
00033
00034 const int kMaxMuon = 10000;
00035 muonpt = new float[kMaxMuon];
00036 muonphi = new float[kMaxMuon];
00037 muoneta = new float[kMaxMuon];
00038 muonet = new float[kMaxMuon];
00039 muone = new float[kMaxMuon];
00040 muonchi2NDF = new float[kMaxMuon];
00041 muoncharge = new float[kMaxMuon];
00042 muonTrkIsoR03 = new float[kMaxMuon];
00043 muonECalIsoR03 = new float[kMaxMuon];
00044 muonHCalIsoR03 = new float[kMaxMuon];
00045 muonD0 = new float[kMaxMuon];
00046 muontype = new int[kMaxMuon];
00047 muonNValidTrkHits = new int[kMaxMuon];
00048 muonNValidMuonHits = new int[kMaxMuon];
00049 const int kMaxMuonL2 = 500;
00050 muonl2pt = new float[kMaxMuonL2];
00051 muonl2phi = new float[kMaxMuonL2];
00052 muonl2eta = new float[kMaxMuonL2];
00053 muonl2dr = new float[kMaxMuonL2];
00054 muonl2dz = new float[kMaxMuonL2];
00055 muonl2vtxz = new float[kMaxMuonL2];
00056 muonl2chg = new int[kMaxMuonL2];
00057 muonl2pterr = new float[kMaxMuonL2];
00058 muonl2iso = new int[kMaxMuonL2];
00059 muonl2nhits = new int[kMaxMuonL2];
00060 muonl2nchambers = new int[kMaxMuonL2];
00061 muonl2nstat = new int[kMaxMuonL2];
00062 muonl21idx = new int[kMaxMuonL2];
00063 const int kMaxMuonL3 = 500;
00064 muonl3pt = new float[kMaxMuonL3];
00065 muonl3phi = new float[kMaxMuonL3];
00066 muonl3eta = new float[kMaxMuonL3];
00067 muonl3dr = new float[kMaxMuonL3];
00068 muonl3dz = new float[kMaxMuonL3];
00069 muonl3vtxz = new float[kMaxMuonL3];
00070 muonl3chg = new int[kMaxMuonL3];
00071 muonl3pterr = new float[kMaxMuonL3];
00072 muonl3iso = new int[kMaxMuonL3];
00073 muonl3trk10iso = new int[kMaxMuonL3];
00074 muonl3nhits = new int[kMaxMuonL3];
00075 muonl3normchi2 = new float[kMaxMuonL3];
00076 muonl3ntrackerhits = new int[kMaxMuonL3];
00077 muonl3nmuonhits = new int[kMaxMuonL3];
00078 muonl32idx = new int[kMaxMuonL3];
00079 muonl3globalpt = new float[kMaxMuonL3];
00080 muonl3globaleta = new float[kMaxMuonL3];
00081 muonl3globalphi = new float[kMaxMuonL3];
00082 muonl3globaldr = new float[kMaxMuonL3];
00083 muonl3globaldz = new float[kMaxMuonL3];
00084 muonl3globalvtxz = new float[kMaxMuonL3];
00085 muonl3globalchg = new int[kMaxMuonL3];
00086 muonl3global2idx = new int[kMaxMuonL3];
00087 const int kMaxTrackerMuon = 500;
00088 trackermuonpt = new float[kMaxTrackerMuon];
00089 trackermuonphi = new float[kMaxTrackerMuon];
00090 trackermuoneta = new float[kMaxTrackerMuon];
00091 trackermuonchg = new int[kMaxTrackerMuon];
00092 trackermuonnhits = new int[kMaxTrackerMuon];
00093 const int kMaxOniaPixel = 500;
00094 oniaPixelpt = new float[kMaxOniaPixel];
00095 oniaPixelphi = new float[kMaxOniaPixel];
00096 oniaPixeleta = new float[kMaxOniaPixel];
00097 oniaPixeldr = new float[kMaxOniaPixel];
00098 oniaPixeldz = new float[kMaxOniaPixel];
00099 oniaPixelchg = new int[kMaxOniaPixel];
00100 oniaPixelHits = new int[kMaxOniaPixel];
00101 oniaPixelNormChi2 = new float[kMaxOniaPixel];
00102 const int kMaxTrackPixel = 500;
00103 oniaTrackpt = new float[kMaxTrackPixel];
00104 oniaTrackphi = new float[kMaxTrackPixel];
00105 oniaTracketa = new float[kMaxTrackPixel];
00106 oniaTrackdr = new float[kMaxTrackPixel];
00107 oniaTrackdz = new float[kMaxTrackPixel];
00108 oniaTrackchg = new int[kMaxTrackPixel];
00109 oniaTrackHits = new int[kMaxTrackPixel];
00110 oniaTrackNormChi2 = new float[kMaxTrackPixel];
00111 const int kMaxMuonL2NoVtx = 500;
00112 muonl2novtxpt = new float[kMaxMuonL2NoVtx];
00113 muonl2novtxphi = new float[kMaxMuonL2NoVtx];
00114 muonl2novtxeta = new float[kMaxMuonL2NoVtx];
00115 muonl2novtxdr = new float[kMaxMuonL2NoVtx];
00116 muonl2novtxdz = new float[kMaxMuonL2NoVtx];
00117 muonl2novtxchg = new int[kMaxMuonL2NoVtx];
00118 muonl2novtxpterr = new float[kMaxMuonL2NoVtx];
00119 muonl2novtxnhits = new int[kMaxMuonL2NoVtx];
00120 muonl2novtxnchambers = new int[kMaxMuonL2NoVtx];
00121 muonl2novtx1idx = new int[kMaxMuonL2NoVtx];
00122 const int kMaxDiMu = 500;
00123 dimudca = new float[kMaxDiMu];
00124 dimu1st = new int[kMaxDiMu];
00125 dimu2nd = new int[kMaxDiMu];
00126 const int kMaxDiMuVtx = 500;
00127 dimuvtx1st = new int[kMaxDiMuVtx];
00128 dimuvtx2nd = new int[kMaxDiMuVtx];
00129 dimuvtxchi2 = new float[kMaxDiMuVtx];
00130 dimuvtxr = new float[kMaxDiMuVtx];
00131 dimuvtxrsig = new float[kMaxDiMuVtx];
00132 dimuvtxroversig = new float[kMaxDiMuVtx];
00133 dimuvtxcosalpha = new float[kMaxDiMuVtx];
00134 dimuvtxmu2dipmax = new float[kMaxDiMuVtx];
00135 dimuvtxmu2dipmin = new float[kMaxDiMuVtx];
00136 dimuvtxmu2dipsigmax = new float[kMaxDiMuVtx];
00137 dimuvtxmu2dipsigmin = new float[kMaxDiMuVtx];
00138
00139
00140 const int kMaxPfmuon = 10000;
00141 pfmuonpt = new float[kMaxPfmuon];
00142 pfmuonphi = new float[kMaxPfmuon];
00143 pfmuoneta = new float[kMaxPfmuon];
00144 pfmuonet = new float[kMaxPfmuon];
00145 pfmuone = new float[kMaxPfmuon];
00146 pfmuoncharge = new float[kMaxPfmuon];
00147
00148
00149 HltTree->Branch("NrecoMuon",&nmuon,"NrecoMuon/I");
00150 HltTree->Branch("recoMuonPt",muonpt,"recoMuonPt[NrecoMuon]/F");
00151 HltTree->Branch("recoMuonPhi",muonphi,"recoMuonPhi[NrecoMuon]/F");
00152 HltTree->Branch("recoMuonEta",muoneta,"recoMuonEta[NrecoMuon]/F");
00153 HltTree->Branch("recoMuonEt",muonet,"recoMuonEt[NrecoMuon]/F");
00154 HltTree->Branch("recoMuonE",muone,"recoMuonE[NrecoMuon]/F");
00155 HltTree->Branch("recoMuonChi2NDF", muonchi2NDF, "recoMuonChi2NDF[NrecoMuon]/F");
00156 HltTree->Branch("recoMuonCharge", muoncharge , "recoMuonCharge[NrecoMuon]/F");
00157 HltTree->Branch("recoMuonTrkIsoR03", muonTrkIsoR03 , "recoMuonTrkIsoR03[NrecoMuon]/F");
00158 HltTree->Branch("recoMuonECalIsoR03", muonECalIsoR03 , "recoMuonECalIsoR03[NrecoMuon]/F");
00159 HltTree->Branch("recoMuonHCalIsoR03", muonHCalIsoR03 , "recoMuonHCalIsoR03[NrecoMuon]/F");
00160 HltTree->Branch("recoMuonD0", muonD0 , "recoMuonD0[NrecoMuon]/F");
00161 HltTree->Branch("recoMuonType", muontype , "recoMuonType[NrecoMuon]/I");
00162 HltTree->Branch("recoMuonNValidTrkHits", muonNValidTrkHits, "recoMuonNValidTrkHits[NrecoMuon]/I");
00163 HltTree->Branch("recoMuonNValidMuonHits", muonNValidMuonHits,"recoMuonNValidMuonHits[NrecoMuon]/I");
00164
00165 HltTree->Branch("NohMuL2",&nmu2cand,"NohMuL2/I");
00166 HltTree->Branch("ohMuL2Pt",muonl2pt,"ohMuL2Pt[NohMuL2]/F");
00167 HltTree->Branch("ohMuL2Phi",muonl2phi,"ohMuL2Phi[NohMuL2]/F");
00168 HltTree->Branch("ohMuL2Eta",muonl2eta,"ohMuL2Eta[NohMuL2]/F");
00169 HltTree->Branch("ohMuL2Chg",muonl2chg,"ohMuL2Chg[NohMuL2]/I");
00170 HltTree->Branch("ohMuL2PtErr",muonl2pterr,"ohMuL2PtErr[NohMuL2]/F");
00171 HltTree->Branch("ohMuL2Iso",muonl2iso,"ohMuL2Iso[NohMuL2]/I");
00172 HltTree->Branch("ohMuL2Dr",muonl2dr,"ohMuL2Dr[NohMuL2]/F");
00173 HltTree->Branch("ohMuL2Dz",muonl2dz,"ohMuL2Dz[NohMuL2]/F");
00174 HltTree->Branch("ohMuL2VtxZ",muonl2vtxz,"ohMuL2VtxZ[NohMuL2]/F");
00175 HltTree->Branch("ohMuL2Nhits",muonl2nhits,"ohMuL2Nhits[NohMuL2]/I");
00176 HltTree->Branch("ohMuL2Nchambers",muonl2nchambers,"ohMuL2Nchambers[NohMuL2]/I");
00177 HltTree->Branch("ohMuL2Nstat",muonl2nstat,"ohMuL2Nstat[NohMuL2]/I");
00178 HltTree->Branch("ohMuL2L1idx",muonl21idx,"ohMuL2L1idx[NohMuL2]/I");
00179 HltTree->Branch("NohMuL3",&nmu3cand,"NohMuL3/I");
00180 HltTree->Branch("ohMuL3Pt",muonl3pt,"ohMuL3Pt[NohMuL3]/F");
00181 HltTree->Branch("ohMuL3Phi",muonl3phi,"ohMuL3Phi[NohMuL3]/F");
00182 HltTree->Branch("ohMuL3Eta",muonl3eta,"ohMuL3Eta[NohMuL3]/F");
00183 HltTree->Branch("ohMuL3Chg",muonl3chg,"ohMuL3Chg[NohMuL3]/I");
00184 HltTree->Branch("ohMuL3PtErr",muonl3pterr,"ohMuL3PtErr[NohMuL3]/F");
00185 HltTree->Branch("ohMuL3Iso",muonl3iso,"ohMuL3Iso[NohMuL3]/I");
00186 HltTree->Branch("ohMuL3Trk10Iso",muonl3trk10iso,"ohMuL3Trk10Iso[NohMuL3]/I");
00187 HltTree->Branch("ohMuL3Dr",muonl3dr,"ohMuL3Dr[NohMuL3]/F");
00188 HltTree->Branch("ohMuL3Dz",muonl3dz,"ohMuL3Dz[NohMuL3]/F");
00189 HltTree->Branch("ohMuL3VtxZ",muonl3vtxz,"ohMuL3VtxZ[NohMuL3]/F");
00190 HltTree->Branch("ohMuL3Nhits",muonl3nhits,"ohMuL3Nhits[NohMuL3]/I");
00191 HltTree->Branch("ohMuL3NormChi2", muonl3normchi2, "ohMuL3NormChi2[NohMuL3]/F");
00192 HltTree->Branch("ohMuL3Ntrackerhits", muonl3ntrackerhits, "ohMuL3Ntrackerhits[NohMuL3]/I");
00193 HltTree->Branch("ohMuL3Nmuonhits", muonl3nmuonhits, "ohMuL3Nmuonhits[NohMuL3]/I");
00194 HltTree->Branch("ohMuL3L2idx",muonl32idx,"ohMuL3L2idx[NohMuL3]/I");
00195 HltTree->Branch("ohMuL3globalPt",muonl3globalpt,"ohMuL3globalPt[NohMuL3]/F");
00196 HltTree->Branch("ohMuL3globalEta",muonl3globaleta,"ohMuL3globalEta[NohMuL3]/F");
00197 HltTree->Branch("ohMuL3globalPhi",muonl3globalphi,"ohMuL3globalPhi[NohMuL3]/F");
00198 HltTree->Branch("ohMuL3globalDr",muonl3globaldr,"ohMuL3globalDr[NohMuL3]/F");
00199 HltTree->Branch("ohMuL3globalDz",muonl3globaldz,"ohMuL3globalDz[NohMuL3]/F");
00200 HltTree->Branch("ohMuL3globalVtxZ",muonl3vtxz,"ohMuL3globalVtxZ[NohMuL3]/F");
00201 HltTree->Branch("ohMuL3globalL2idx",muonl3global2idx,"ohMuL3globalL2idx[NohMuL3]/I");
00202
00203 HltTree->Branch("NohOniaPixel",&nOniaPixelCand,"NohOniaPixel/I");
00204 HltTree->Branch("ohOniaPixelPt",oniaPixelpt,"ohOniaPixelPt[NohOniaPixel]/F");
00205 HltTree->Branch("ohOniaPixelPhi",oniaPixelphi,"ohOniaPixelPhi[NohOniaPixel]/F");
00206 HltTree->Branch("ohOniaPixelEta",oniaPixeleta,"ohOniaPixelEta[NohOniaPixel]/F");
00207 HltTree->Branch("ohOniaPixelChg",oniaPixelchg,"ohOniaPixelChg[NohOniaPixel]/I");
00208 HltTree->Branch("ohOniaPixelDr",oniaPixeldr,"ohOniaPixelDr[NohOniaPixel]/F");
00209 HltTree->Branch("ohOniaPixelDz",oniaPixeldz,"ohOniaPixelDz[NohOniaPixel]/F");
00210 HltTree->Branch("ohOniaPixelHits",oniaPixelHits,"ohOniaPixelHits[NohOniaPixel]/I");
00211 HltTree->Branch("ohOniaPixelNormChi2",oniaPixelNormChi2,"ohOniaPixelNormChi2[NohOniaPixel]/F");
00212 HltTree->Branch("NohOniaTrack",&nOniaTrackCand,"NohOniaTrack/I");
00213 HltTree->Branch("ohOniaTrackPt",oniaTrackpt,"ohOniaTrackPt[NohOniaTrack]/F");
00214 HltTree->Branch("ohOniaTrackPhi",oniaTrackphi,"ohOniaTrackPhi[NohOniaTrack]/F");
00215 HltTree->Branch("ohOniaTrackEta",oniaTracketa,"ohOniaTrackEta[NohOniaTrack]/F");
00216 HltTree->Branch("ohOniaTrackChg",oniaTrackchg,"ohOniaTrackChg[NohOniaTrack]/I");
00217 HltTree->Branch("ohOniaTrackDr",oniaTrackdr,"ohOniaTrackDr[NohOniaTrack]/F");
00218 HltTree->Branch("ohOniaTrackDz",oniaTrackdz,"ohOniaTrackDz[NohOniaTrack]/F");
00219 HltTree->Branch("ohOniaTrackHits",oniaTrackHits,"ohOniaTrackHits[NohOniaTrack]/I");
00220 HltTree->Branch("ohOniaTrackNormChi2",oniaTrackNormChi2,"ohOniaTrackNormChi2[NohOniaTrack]/F");
00221 HltTree->Branch("NohMuL2NoVtx",&nmu2cand,"NohMuL2NoVtx/I");
00222 HltTree->Branch("ohMuL2NoVtxPt",muonl2novtxpt,"ohMuL2NoVtxPt[NohMuL2NoVtx]/F");
00223 HltTree->Branch("ohMuL2NoVtxPhi",muonl2novtxphi,"ohMuL2NoVtxPhi[NohMuL2NoVtx]/F");
00224 HltTree->Branch("ohMuL2NoVtxEta",muonl2novtxeta,"ohMuL2NoVtxEta[NohMuL2NoVtx]/F");
00225 HltTree->Branch("ohMuL2NoVtxChg",muonl2novtxchg,"ohMuL2NoVtxChg[NohMuL2NoVtx]/I");
00226 HltTree->Branch("ohMuL2NoVtxPtErr",muonl2novtxpterr,"ohMuL2NoVtxPtErr[NohMuL2NoVtx]/F");
00227 HltTree->Branch("ohMuL2NoVtxDr",muonl2novtxdr,"ohMuL2NoVtxDr[NohMuL2NoVtx]/F");
00228 HltTree->Branch("ohMuL2NoVtxDz",muonl2novtxdz,"ohMuL2NoVtxDz[NohMuL2NoVtx]/F");
00229 HltTree->Branch("ohMuL2NoVtxNhits",muonl2novtxnhits,"ohMuL2NoVtxNhits[NohMuL2NoVtx]/I");
00230 HltTree->Branch("ohMuL2NoVtxNchambers",muonl2novtxnchambers,"ohMuL2NoVtxNchambers[NohMuL2NoVtx]/I");
00231 HltTree->Branch("ohMuL2NoVtxL1idx",muonl2novtx1idx,"ohMuL2NoVtxL1idx[NohMuL2NoVtx]/I");
00232 HltTree->Branch("NohDiMu",&nDiMu,"NohDiMu/I");
00233 HltTree->Branch("ohDiMuDCA",dimudca,"ohDiMuDCA[NohDiMu]/F");
00234 HltTree->Branch("ohDiMu1st",dimu1st,"ohDiMu1st[NohDiMu]/I");
00235 HltTree->Branch("ohDiMu2nd",dimu2nd,"ohDiMu2nd[NohDiMu]/I");
00236 HltTree->Branch("NohDiMuVtx",&nDiMuVtx,"NohDiMuVtx/I");
00237 HltTree->Branch("ohDiMuVtx1st",dimuvtx1st,"ohDiMuVtx1st[NohDiMuVtx]/I");
00238 HltTree->Branch("ohDiMuVtx2nd",dimuvtx2nd,"ohDiMuVtx2nd[NohDiMuVtx]/I");
00239 HltTree->Branch("ohDiMuVtxChi2",dimuvtxchi2,"ohDiMuVtxChi2[NohDiMuVtx]/F");
00240 HltTree->Branch("ohDiMuVtxR",dimuvtxr,"ohDiMuVtxR[NohDiMuVtx]/F");
00241 HltTree->Branch("ohDiMuVtxRSig",dimuvtxrsig,"ohDiMuVtxRSig[NohDiMuVtx]/F");
00242 HltTree->Branch("ohDiMuVtxROverSig",dimuvtxroversig,"ohDiMuVtxROverSig[NohDiMuVtx]/F");
00243 HltTree->Branch("ohDiMuVtxCosAlpha",dimuvtxcosalpha,"ohDiMuVtxCosAlpha[NohDiMuVtx]/F");
00244 HltTree->Branch("ohDiMuVtxMu2DIpMax",dimuvtxmu2dipmax,"ohDiMuVtxMu2DIpMax[NohDiMuVtx]/F");
00245 HltTree->Branch("ohDiMuVtxMu2DIpMin",dimuvtxmu2dipmin,"ohDiMuVtxMu2DIpMin[NohDiMuVtx]/F");
00246 HltTree->Branch("ohDiMuVtxMu2DIpSigMax",dimuvtxmu2dipsigmax,"ohDiMuVtxMu2DIpSigMax[NohDiMuVtx]/F");
00247 HltTree->Branch("ohDiMuVtxMu2DIpSigMin",dimuvtxmu2dipsigmin,"ohDiMuVtxMu2DIpSigMin[NohDiMuVtx]/F");
00248 HltTree->Branch("NohTrackerMuon",&ntrackermuoncand,"NohTrackerMuon/I");
00249 HltTree->Branch("ohTrackerMuonPt",trackermuonpt,"ohTrackerMuonPt[NohTrackerMuon]/F");
00250 HltTree->Branch("ohTrackerMuonPhi",trackermuonphi,"ohTrackerMuonPhi[NohTrackerMuon]/F");
00251 HltTree->Branch("ohTrackerMuonEta",trackermuoneta,"ohTrackerMuonEta[NohTrackerMuon]/F");
00252 HltTree->Branch("ohTrackerMuonChg",trackermuonchg,"ohTrackerMuonChg[NohTrackerMuon]/I");
00253 HltTree->Branch("ohTrackerMuonNhits",trackermuonnhits,"ohTrackerMuonNhits[NohTrackerMuon]/I");
00254
00255 HltTree->Branch("NpfMuon",&npfmuon,"NpfMuon/I");
00256 HltTree->Branch("pfMuonPt",pfmuonpt,"pfMuonPt[NpfMuon]/F");
00257 HltTree->Branch("pfMuonPhi",pfmuonphi,"pfMuonPhi[NpfMuon]/F");
00258 HltTree->Branch("pfMuonEta",pfmuoneta,"pfMuonEta[NpfMuon]/F");
00259 HltTree->Branch("pfMuonEt",pfmuonet,"pfMuonEt[NpfMuon]/F");
00260 HltTree->Branch("pfMuonE",pfmuone,"pfMuonE[NpfMuon]/F");
00261 HltTree->Branch("pfMuonCharge", pfmuoncharge , "pfMuonCharge[NpfMuon]/F");
00262
00263
00264 }
00265
00266
00267 void HLTMuon::analyze(const edm::Handle<reco::MuonCollection> & Muon,
00268 const edm::Handle<reco::PFCandidateCollection> & pfMuon,
00269 const edm::Handle<l1extra::L1MuonParticleCollection> & MuCands1,
00270 const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands2,
00271 const edm::Handle<edm::ValueMap<bool> > & isoMap2,
00272 const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands3,
00273 const edm::Handle<edm::ValueMap<bool> > & isoMap3,
00274 const edm::Handle<edm::ValueMap<bool> > & isoTrk10Map3,
00275 const edm::Handle<reco::RecoChargedCandidateCollection> & oniaPixelCands,
00276 const edm::Handle<reco::RecoChargedCandidateCollection> & oniaTrackCands,
00277 const edm::Handle<reco::VertexCollection> & DiMuVtxCands3,
00278 const edm::Handle<reco::RecoChargedCandidateCollection> & MuNoVtxCands2,
00279 const edm::Handle<reco::MuonCollection> & trkmucands,
00280 const edm::ESHandle<MagneticField> & theMagField,
00281 const edm::Handle<reco::BeamSpot> & recoBeamSpotHandle,
00282 TTree* HltTree) {
00283
00284 reco::BeamSpot::Point BSPosition(0,0,0);
00285 BSPosition = recoBeamSpotHandle->position();
00286 const GlobalPoint theBeamSpot = GlobalPoint(recoBeamSpotHandle->position().x(),
00287 recoBeamSpotHandle->position().y(),
00288 recoBeamSpotHandle->position().z());
00289 reco::BeamSpot vtxBS = *recoBeamSpotHandle;
00290
00291
00292
00293 if (Muon.isValid()) {
00294 reco::MuonCollection mymuons;
00295 mymuons = * Muon;
00296 std::sort(mymuons.begin(),mymuons.end(),PtGreater());
00297 nmuon = mymuons.size();
00298 typedef reco::MuonCollection::const_iterator muiter;
00299 int imu=0;
00300 for (muiter i=mymuons.begin(); i!=mymuons.end(); i++)
00301 {
00302 muonpt[imu] = i->pt();
00303 muonphi[imu] = i->phi();
00304 muoneta[imu] = i->eta();
00305 muonet[imu] = i->et();
00306 muone[imu] = i->energy();
00307 muontype[imu] = i->type();
00308 muoncharge[imu] = i->charge();
00309 muonTrkIsoR03[imu] = i->isolationR03().sumPt;
00310 muonECalIsoR03[imu] = i->isolationR03().emEt;
00311 muonHCalIsoR03[imu] = i->isolationR03().hadEt;
00312
00313
00314 if (i->globalTrack().isNonnull())
00315 {
00316 muonchi2NDF[imu] = i->globalTrack()->normalizedChi2();
00317 muonD0[imu] = i->globalTrack()->dxy(BSPosition);
00318 }
00319 else
00320 {
00321 muonchi2NDF[imu] = -99.;
00322 muonD0[imu] = -99.;}
00323
00324 if (i->innerTrack().isNonnull()) muonNValidTrkHits[imu] = i->innerTrack()->numberOfValidHits();
00325 else muonNValidTrkHits[imu] = -99;
00326
00327 if (i->isGlobalMuon()!=0) muonNValidMuonHits[imu] = i->globalTrack()->hitPattern().numberOfValidMuonHits();
00328 else muonNValidMuonHits[imu] = -99;
00329
00330 imu++;
00331 }
00332 }
00333 else {nmuon = 0;}
00334
00335 l1extra::L1MuonParticleCollection myMucands1;
00336 myMucands1 = * MuCands1;
00337
00338 std::sort(myMucands1.begin(),myMucands1.end(),PtGreater());
00339
00341
00342
00343 reco::RecoChargedCandidateCollection myMucands2;
00344 if (MuCands2.isValid()) {
00345
00346 myMucands2 = * MuCands2;
00347 std::sort(myMucands2.begin(),myMucands2.end(),PtGreater());
00348 nmu2cand = myMucands2.size();
00349 typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00350 int imu2c=0;
00351 for (cand i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00352 reco::TrackRef tk = i->get<reco::TrackRef>();
00353
00354 muonl2pt[imu2c] = tk->pt();
00355
00356 muonl2eta[imu2c] = tk->eta();
00357 muonl2phi[imu2c] = tk->phi();
00358
00359
00360
00361
00362 muonl2dr[imu2c] = fabs(tk->dxy(BSPosition));
00363
00364
00365
00366 muonl2dz[imu2c] = tk->dz(BSPosition);
00367 muonl2vtxz[imu2c] = tk->dz();
00368 muonl2nhits[imu2c] = tk->numberOfValidHits();
00369 muonl2nchambers[imu2c] = validChambers(tk);
00370 muonl2nstat[imu2c] = tk->hitPattern().muonStationsWithAnyHits();
00371
00372
00373
00374
00375
00376
00377
00378 double l2_err0 = tk->error(0);
00379 double l2_abspar0 = fabs(tk->parameter(0));
00380
00381
00382
00383
00384
00385
00386
00387
00388
00389
00390
00391
00392 muonl2pterr[imu2c] = l2_err0/l2_abspar0;
00393 muonl2chg[imu2c] = tk->charge();
00394
00395 if (isoMap2.isValid()){
00396
00397 edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap2)[tk];
00398 muonl2iso[imu2c] = muon1IsIsolated;
00399 }
00400 else {muonl2iso[imu2c] = -999;}
00401
00402 l1extra::L1MuonParticleRef l1;
00403 int il2 = 0;
00404
00405 l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00406 il2++;
00407 int imu1idx = 0;
00408 if (MuCands1.isValid()) {
00409 typedef l1extra::L1MuonParticleCollection::const_iterator candl1;
00410 for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) {
00411 if((j->pt() == l1->pt()) &&
00412 (j->eta() == l1->eta()) &&
00413 (j->phi() == l1->phi()) &&
00414 (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00415 {break;}
00416
00417
00418 imu1idx++;
00419 }
00420 }
00421 else {imu1idx = -999;}
00422 muonl21idx[imu2c] = imu1idx;
00423
00424 imu2c++;
00425 }
00426 }
00427 else {nmu2cand = 0;}
00428
00429
00430 reco::RecoChargedCandidateCollection myMucands3;
00431 if (MuCands3.isValid()) {
00432 int k = 0;
00433 myMucands3 = * MuCands3;
00434 std::sort(myMucands3.begin(),myMucands3.end(),PtGreater());
00435 nmu3cand = myMucands3.size();
00436 typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00437 int imu3c=0;
00438 int idimuc=0;
00439 for (cand i=myMucands3.begin(); i!=myMucands3.end(); i++) {
00440 reco::TrackRef tk = i->get<reco::TrackRef>();
00441
00442 reco::RecoChargedCandidateRef candref = reco::RecoChargedCandidateRef(MuCands3,k);
00443
00444 reco::TrackRef staTrack;
00445 typedef reco::MuonTrackLinksCollection::const_iterator l3muon;
00446 int il3 = 0;
00447
00448 staTrack = tk->seedRef().castTo<edm::Ref< L3MuonTrajectorySeedCollection> >()->l2Track();
00449 il3++;
00450 int imu2idx = 0;
00451 if (MuCands2.isValid()) {
00452 typedef reco::RecoChargedCandidateCollection::const_iterator candl2;
00453 for (candl2 i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00454 reco::TrackRef tkl2 = i->get<reco::TrackRef>();
00455 if ( tkl2 == staTrack ) {break;}
00456 imu2idx++;
00457 }
00458 }
00459 else {imu2idx = -999;}
00460 muonl3global2idx[imu3c] = imu2idx;
00461 muonl32idx[imu3c] = imu2idx;
00462
00463 muonl3globalpt[imu3c] = tk->pt();
00464 muonl3pt[imu3c] = candref->pt();
00465
00466 muonl3globaleta[imu3c] = tk->eta();
00467 muonl3globalphi[imu3c] = tk->phi();
00468 muonl3eta[imu3c] = candref->eta();
00469 muonl3phi[imu3c] = candref->phi();
00470
00471
00472
00473
00474 muonl3dr[imu3c] = fabs(tk->dxy(BSPosition));
00475
00476
00477
00478 muonl3dz[imu3c] = tk->dz(BSPosition);
00479 muonl3vtxz[imu3c] = tk->dz();
00480 muonl3nhits[imu3c] = tk->numberOfValidHits();
00481
00482
00483
00484
00485
00486
00487
00488 double l3_err0 = tk->error(0);
00489 double l3_abspar0 = fabs(tk->parameter(0));
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502 muonl3pterr[imu3c] = l3_err0/l3_abspar0;
00503 muonl3globalchg[imu3c] = tk->charge();
00504 muonl3chg[imu3c] = candref->charge();
00505
00506 muonl3normchi2[imu3c] = tk->normalizedChi2();
00507 muonl3ntrackerhits[imu3c] = tk->hitPattern().numberOfValidTrackerHits();
00508 muonl3nmuonhits[imu3c] = tk->hitPattern().numberOfValidMuonHits();
00509
00510 if (isoMap3.isValid()){
00511
00512 edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap3)[tk];
00513 muonl3iso[imu3c] = muon1IsIsolated;
00514 }
00515 else {muonl3iso[imu3c] = -999;}
00516
00517 if (isoTrk10Map3.isValid()){
00518
00519 edm::ValueMap<bool> ::value_type muon1IsTrk10Isolated = (*isoTrk10Map3)[tk];
00520 muonl3trk10iso[imu3c] = muon1IsTrk10Isolated;
00521 }
00522 else {muonl3trk10iso[imu3c] =-999;}
00523
00524
00525 int imu3c2nd = imu3c + 1;
00526
00527 for (cand j=i; j!=myMucands3.end(); j++) if (i!=j) {
00528 reco::TrackRef tk2nd = j->get<reco::TrackRef>();
00529 reco::TransientTrack transMu1(*tk, &(*theMagField) );
00530 reco::TransientTrack transMu2(*tk2nd, &(*theMagField) );
00531 TrajectoryStateClosestToPoint mu1TS = transMu1.impactPointTSCP();
00532 TrajectoryStateClosestToPoint mu2TS = transMu2.impactPointTSCP();
00533 if (mu1TS.isValid() && mu2TS.isValid()) {
00534 ClosestApproachInRPhi cApp;
00535 cApp.calculate(mu1TS.theState(), mu2TS.theState());
00536 if (cApp.status()) {
00537 dimudca[idimuc] = cApp.distance();
00538 dimu1st[idimuc] = imu3c;
00539 dimu2nd[idimuc] = imu3c2nd;
00540 idimuc++;
00541 }
00542 }
00543 imu3c2nd++;
00544 }
00545
00546 imu3c++;
00547 k++;
00548 }
00549 nDiMu = idimuc;
00550 }
00551
00552 else {nmu3cand = 0; nDiMu = 0;}
00553
00554 reco::VertexCollection myDimuvtxcands3;
00555 if (DiMuVtxCands3.isValid()) {
00556 myDimuvtxcands3 = * DiMuVtxCands3;
00557 nDiMuVtx = myDimuvtxcands3.size();
00558 typedef reco::VertexCollection::const_iterator cand;
00559 int idimu3c=0;
00560 for (cand ivtx = myDimuvtxcands3.begin(); ivtx != myDimuvtxcands3.end(); ++ivtx) {
00561 dimuvtxchi2[idimu3c] = ivtx->normalizedChi2();
00562 reco::Vertex::trackRef_iterator trackIt = ivtx->tracks_begin();
00563 reco::TrackRef vertextkRef1 = (*trackIt).castTo<reco::TrackRef>();
00564 ++trackIt;
00565 reco::TrackRef vertextkRef2 = (*trackIt).castTo<reco::TrackRef>();
00566 dimuvtx2nd[idimu3c] = -1; dimuvtx1st[idimu3c] = -1;
00567 for (int j=0 ; j<nmu3cand ; j++){
00568 if(fabs(muonl3pt[j] - vertextkRef1->pt()) < 0.0001 && fabs(muonl3eta[j] - vertextkRef1->eta()) < 0.0001 && fabs(muonl3phi[j] - vertextkRef1->phi()) < 0.0001) dimuvtx1st[idimu3c] = j;
00569 if(fabs(muonl3pt[j] - vertextkRef2->pt()) < 0.0001 && fabs(muonl3eta[j] - vertextkRef2->eta()) < 0.0001 && fabs(muonl3phi[j] - vertextkRef2->phi()) < 0.0001) dimuvtx2nd[idimu3c] = j;
00570 }
00571 math::XYZVector pperp(vertextkRef1->px() + vertextkRef2->px(),
00572 vertextkRef1->py() + vertextkRef2->py(),
00573 0.);
00574 reco::Vertex::Point vpoint = ivtx->position();
00575 GlobalPoint vtxPos (vpoint.x(), vpoint.y(), vpoint.z());
00576 reco::Vertex::Error verr = ivtx->error();
00577 GlobalError vtxErr (verr.At(0,0),verr.At(1,0),verr.At(1,1),verr.At(2,0),verr.At(2,1),verr.At(2,2));
00578 GlobalPoint vtxDisFromBS(-1*((vtxBS.x0() - vtxPos.x()) + (vtxPos.z() - vtxBS.z0())*vtxBS.dxdz()),
00579 -1*((vtxBS.y0() - vtxPos.y()) + (vtxPos.z() - vtxBS.z0())*vtxBS.dydz()), 0.0);
00580 dimuvtxr[idimu3c] = vtxDisFromBS.perp();
00581 dimuvtxrsig[idimu3c] = sqrt(vtxErr.rerr(vtxDisFromBS));
00582 dimuvtxroversig[idimu3c] = dimuvtxr[idimu3c]/dimuvtxrsig[idimu3c];
00583 reco::Vertex::Point vperp(vtxDisFromBS.x(),vtxDisFromBS.y(),0.);
00584 dimuvtxcosalpha[idimu3c] = vperp.Dot(pperp)/(vperp.R()*pperp.R());
00585 float mu1ip = -1.0;
00586 float mu2ip = -1.0;
00587 float mu1ipsig = -1.0;
00588 float mu2ipsig = -1.0;
00589 reco::TransientTrack transMu1(*vertextkRef1, &(*theMagField) );
00590 TrajectoryStateClosestToPoint trajMu1BS = transMu1.trajectoryStateClosestToPoint(theBeamSpot);
00591 if(trajMu1BS.isValid()){
00592 mu1ip = fabs(trajMu1BS.perigeeParameters().transverseImpactParameter());
00593 if(trajMu1BS.hasError()) mu1ipsig = mu1ip/trajMu1BS.perigeeError().transverseImpactParameterError();
00594 }
00595 reco::TransientTrack transMu2(*vertextkRef2, &(*theMagField) );
00596 TrajectoryStateClosestToPoint trajMu2BS = transMu2.trajectoryStateClosestToPoint(theBeamSpot);
00597 if(trajMu2BS.isValid()){
00598 mu2ip = fabs(trajMu2BS.perigeeParameters().transverseImpactParameter());
00599 if(trajMu2BS.hasError()) mu2ipsig = mu2ip/trajMu2BS.perigeeError().transverseImpactParameterError();
00600 }
00601 dimuvtxmu2dipmax[idimu3c] = fmax(mu1ip,mu2ip);
00602 dimuvtxmu2dipmin[idimu3c] = fmin(mu1ip,mu2ip);
00603 dimuvtxmu2dipsigmax[idimu3c] = fmax(mu1ipsig,mu2ipsig);
00604 dimuvtxmu2dipsigmin[idimu3c] = fmin(mu1ipsig,mu2ipsig);
00605 }
00606
00607
00608 }
00609 else {nDiMuVtx = 0;}
00610
00611
00612
00613 reco::RecoChargedCandidateCollection muNoVtxMucands2;
00614 if (MuNoVtxCands2.isValid()) {
00615 muNoVtxMucands2 = * MuNoVtxCands2;
00616 std::sort(muNoVtxMucands2.begin(),muNoVtxMucands2.end(),PtGreater());
00617 nmu2cand = muNoVtxMucands2.size();
00618 typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00619 int imu2c=0;
00620 for (cand i=muNoVtxMucands2.begin(); i!=muNoVtxMucands2.end(); i++) {
00621 reco::TrackRef tk = i->get<reco::TrackRef>();
00622
00623 muonl2novtxpt[imu2c] = tk->pt();
00624 muonl2novtxeta[imu2c] = tk->eta();
00625 muonl2novtxphi[imu2c] = tk->phi();
00626 muonl2novtxdr[imu2c] = fabs(tk->dxy(BSPosition));
00627 muonl2novtxdz[imu2c] = tk->dz(BSPosition);
00628 muonl2novtxnhits[imu2c] = tk->numberOfValidHits();
00629 muonl2novtxnchambers[imu2c] = validChambers(tk);
00630
00631 double l2_err0 = tk->error(0);
00632 double l2_abspar0 = fabs(tk->parameter(0));
00633
00634 muonl2novtxpterr[imu2c] = l2_err0/l2_abspar0;
00635 muonl2novtxchg[imu2c] = tk->charge();
00636
00637 l1extra::L1MuonParticleRef l1;
00638 int il2 = 0;
00639
00640 l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00641 il2++;
00642 int imu1idx = 0;
00643 if (MuCands1.isValid()) {
00644 typedef l1extra::L1MuonParticleCollection::const_iterator candl1;
00645 for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) {
00646 if((j->pt() == l1->pt()) &&
00647 (j->eta() == l1->eta()) &&
00648 (j->phi() == l1->phi()) &&
00649 (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00650 {break;}
00651 imu1idx++;
00652 }
00653 }
00654 else {imu1idx = -999;}
00655 muonl2novtx1idx[imu2c] = imu1idx;
00656
00657 imu2c++;
00658 }
00659 }
00660 else {nmu2cand = 0;}
00661
00662
00663
00664
00665 reco::RecoChargedCandidateCollection myOniaPixelCands;
00666 if (oniaPixelCands.isValid()) {
00667 myOniaPixelCands = * oniaPixelCands;
00668 std::sort(myOniaPixelCands.begin(),myOniaPixelCands.end(),PtGreater());
00669 nOniaPixelCand = myOniaPixelCands.size();
00670 typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00671 int ic=0;
00672 for (cand i=myOniaPixelCands.begin(); i!=myOniaPixelCands.end(); i++) {
00673 reco::TrackRef tk = i->get<reco::TrackRef>();
00674
00675 oniaPixelpt[ic] = tk->pt();
00676 oniaPixeleta[ic] = tk->eta();
00677 oniaPixelphi[ic] = tk->phi();
00678 oniaPixeldr[ic] = tk->dxy(BSPosition);
00679 oniaPixeldz[ic] = tk->dz(BSPosition);
00680 oniaPixelchg[ic] = tk->charge();
00681 oniaPixelHits[ic] = tk->numberOfValidHits();
00682 oniaPixelNormChi2[ic] = tk->normalizedChi2();
00683
00684 ic++;
00685 }
00686 }
00687 else {nOniaPixelCand = 0;}
00688
00689
00690 reco::RecoChargedCandidateCollection myOniaTrackCands;
00691 if (oniaTrackCands.isValid()) {
00692 myOniaTrackCands = * oniaTrackCands;
00693 std::sort(myOniaTrackCands.begin(),myOniaTrackCands.end(),PtGreater());
00694 nOniaTrackCand = myOniaTrackCands.size();
00695 typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00696 int ic=0;
00697 for (cand i=myOniaTrackCands.begin(); i!=myOniaTrackCands.end(); i++) {
00698 reco::TrackRef tk = i->get<reco::TrackRef>();
00699
00700 oniaTrackpt[ic] = tk->pt();
00701 oniaTracketa[ic] = tk->eta();
00702 oniaTrackphi[ic] = tk->phi();
00703 oniaTrackdr[ic] = tk->dxy(BSPosition);
00704 oniaTrackdz[ic] = tk->dz(BSPosition);
00705 oniaTrackchg[ic] = tk->charge();
00706 oniaTrackHits[ic] = tk->numberOfValidHits();
00707 oniaTrackNormChi2[ic] = tk->normalizedChi2();
00708
00709 ic++;
00710 }
00711 }
00712 else {nOniaTrackCand = 0;}
00713
00714
00715
00716 if(trkmucands.isValid()) {
00717 int itrackermuc=0;
00718 for ( unsigned int i=0; i<trkmucands->size(); ++i ){
00719 const reco::Muon& muon(trkmucands->at(i));
00720 if (muon.isTrackerMuon()) {
00721 trackermuonpt[itrackermuc] = muon.pt();
00722 trackermuoneta[itrackermuc] = muon.eta();
00723 trackermuonphi[itrackermuc] = muon.phi();
00724 trackermuonchg[itrackermuc] = muon.charge();
00725 if ( !muon.innerTrack().isNull() ){
00726 trackermuonnhits[itrackermuc] = muon.innerTrack()->numberOfValidHits();
00727 }
00728 itrackermuc++;
00729 }
00730 }
00731 ntrackermuoncand=itrackermuc;
00732 }
00733 else {ntrackermuoncand = 0;}
00734
00736
00737 if (pfMuon.isValid()) {
00738 reco::PFCandidateCollection mypfmuons;
00739 mypfmuons = * pfMuon;
00740 std::sort(mypfmuons.begin(),mypfmuons.end(),PtGreater());
00741 npfmuon = mypfmuons.size();
00742 typedef reco::PFCandidateCollection::const_iterator muiter;
00743 int ipfmu=0;
00744 for (muiter i=mypfmuons.begin(); i!=mypfmuons.end(); i++)
00745 {
00746 pfmuonpt[ipfmu] = i->pt();
00747 pfmuonphi[ipfmu] = i->phi();
00748 pfmuoneta[ipfmu] = i->eta();
00749 pfmuonet[ipfmu] = i->et();
00750 pfmuone[ipfmu] = i->energy();
00751 pfmuoncharge[ipfmu] = i->charge();
00752
00753 ipfmu++;
00754 }
00755 }
00756 else {npfmuon = 0;}
00757
00758 }
00759
00760 int HLTMuon::validChambers(const reco::TrackRef & track)
00761 {
00762
00763 std::map<uint32_t,int> DTchambers;
00764 std::map<uint32_t,int> CSCchambers;
00765
00766 for (trackingRecHit_iterator hit = track->recHitsBegin(); hit != track->recHitsEnd(); ++hit) {
00767 if( !((*hit)->isValid()) ) continue;
00768
00769 DetId id = (*hit)->geographicalId();
00770
00771 if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::DT) {
00772
00773 uint32_t index = DTChamberId(id).rawId();
00774
00775 if (DTchambers.find(index) == DTchambers.end()) {
00776 DTchambers[index] = 0;
00777 }
00778 DTchambers[index]++;
00779 }
00780
00781 else if (id.det() == DetId::Muon && id.subdetId() == MuonSubdetId::CSC) {
00782
00783 CSCDetId id2(id);
00784 uint32_t index = CSCDetId(id2.endcap(), id2.station(), id2.ring(), id2.chamber(), 0);
00785
00786 if (CSCchambers.find(index) == CSCchambers.end()) {
00787 CSCchambers[index] = 0;
00788 }
00789 CSCchambers[index]++;
00790 }
00791 }
00792
00793
00794 int validChambers = 0;
00795
00796 int minDThits = 1;
00797 int minCSChits = 1;
00798
00799 for (std::map<uint32_t,int>::const_iterator iter = DTchambers.begin(); iter != DTchambers.end(); ++iter) {
00800 if (iter->second >= minDThits) {
00801 validChambers++;
00802 }
00803 }
00804 for (std::map<uint32_t,int>::const_iterator iter = CSCchambers.begin(); iter != CSCchambers.end(); ++iter) {
00805 if (iter->second >= minCSChits) {
00806 validChambers++;
00807 }
00808 }
00809 return validChambers;
00810 }