Reference Manual

Go to the documentation of this file.
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   //set parameter defaults 
00018   _Monte=false;
00019   _Debug=false;
00020 }
00021 
00022 /*  Setup the analysis to put the branch-variables into the tree. */
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   const int kMaxTrackerMuon = 500;
00080   trackermuonpt = new float[kMaxTrackerMuon];
00081   trackermuonphi = new float[kMaxTrackerMuon];
00082   trackermuoneta = new float[kMaxTrackerMuon];
00083   trackermuonchg = new int[kMaxTrackerMuon];
00084   trackermuonnhits = new int[kMaxTrackerMuon];
00085   const int kMaxOniaPixel = 500;
00086   oniaPixelpt = new float[kMaxOniaPixel];
00087   oniaPixelphi = new float[kMaxOniaPixel];
00088   oniaPixeleta = new float[kMaxOniaPixel];
00089   oniaPixeldr = new float[kMaxOniaPixel];
00090   oniaPixeldz = new float[kMaxOniaPixel];
00091   oniaPixelchg = new int[kMaxOniaPixel];
00092   oniaPixelHits = new int[kMaxOniaPixel];
00093   oniaPixelNormChi2 = new float[kMaxOniaPixel];
00094   const int kMaxTrackPixel = 500;
00095   oniaTrackpt = new float[kMaxTrackPixel];
00096   oniaTrackphi = new float[kMaxTrackPixel];
00097   oniaTracketa = new float[kMaxTrackPixel];
00098   oniaTrackdr = new float[kMaxTrackPixel];
00099   oniaTrackdz = new float[kMaxTrackPixel];
00100   oniaTrackchg = new int[kMaxTrackPixel];
00101   oniaTrackHits = new int[kMaxTrackPixel];
00102   oniaTrackNormChi2 = new float[kMaxTrackPixel];
00103   const int kMaxMuonL2NoVtx = 500; 
00104   muonl2novtxpt = new float[kMaxMuonL2NoVtx]; 
00105   muonl2novtxphi = new float[kMaxMuonL2NoVtx]; 
00106   muonl2novtxeta = new float[kMaxMuonL2NoVtx]; 
00107   muonl2novtxdr = new float[kMaxMuonL2NoVtx]; 
00108   muonl2novtxdz = new float[kMaxMuonL2NoVtx]; 
00109   muonl2novtxchg = new int[kMaxMuonL2NoVtx]; 
00110   muonl2novtxpterr = new float[kMaxMuonL2NoVtx]; 
00111   muonl2novtxnhits = new int[kMaxMuonL2NoVtx];
00112   muonl2novtxnchambers = new int[kMaxMuonL2NoVtx];
00113   muonl2novtx1idx = new int[kMaxMuonL2NoVtx]; 
00114   const int kMaxDiMu = 500;
00115   dimudca = new float[kMaxDiMu];
00116   dimu1st = new int[kMaxDiMu];
00117   dimu2nd = new int[kMaxDiMu];
00118   const int kMaxDiMuVtx = 500;
00119   dimuvtx1st = new int[kMaxDiMuVtx];
00120   dimuvtx2nd = new int[kMaxDiMuVtx];
00121   dimuvtxchi2 = new float[kMaxDiMuVtx];
00122   dimuvtxr = new float[kMaxDiMuVtx];
00123   dimuvtxrsig = new float[kMaxDiMuVtx];
00124   dimuvtxroversig = new float[kMaxDiMuVtx];
00125   dimuvtxcosalpha = new float[kMaxDiMuVtx];
00126   dimuvtxmu2dipmax = new float[kMaxDiMuVtx];
00127   dimuvtxmu2dipmin = new float[kMaxDiMuVtx];
00128   dimuvtxmu2dipsigmax = new float[kMaxDiMuVtx];
00129   dimuvtxmu2dipsigmin = new float[kMaxDiMuVtx];
00130 
00131   // Muon-specific branches of the tree 
00132   HltTree->Branch("NrecoMuon",&nmuon,"NrecoMuon/I");
00133   HltTree->Branch("recoMuonPt",muonpt,"recoMuonPt[NrecoMuon]/F");
00134   HltTree->Branch("recoMuonPhi",muonphi,"recoMuonPhi[NrecoMuon]/F");
00135   HltTree->Branch("recoMuonEta",muoneta,"recoMuonEta[NrecoMuon]/F");
00136   HltTree->Branch("recoMuonEt",muonet,"recoMuonEt[NrecoMuon]/F");
00137   HltTree->Branch("recoMuonE",muone,"recoMuonE[NrecoMuon]/F");
00138   HltTree->Branch("recoMuonChi2NDF",        muonchi2NDF,       "recoMuonChi2NDF[NrecoMuon]/F");
00139   HltTree->Branch("recoMuonCharge",         muoncharge  ,      "recoMuonCharge[NrecoMuon]/F");
00140   HltTree->Branch("recoMuonTrkIsoR03",      muonTrkIsoR03 ,    "recoMuonTrkIsoR03[NrecoMuon]/F");
00141   HltTree->Branch("recoMuonECalIsoR03",     muonECalIsoR03 ,   "recoMuonECalIsoR03[NrecoMuon]/F");
00142   HltTree->Branch("recoMuonHCalIsoR03",     muonHCalIsoR03 ,   "recoMuonHCalIsoR03[NrecoMuon]/F");
00143   HltTree->Branch("recoMuonD0",             muonD0 ,           "recoMuonD0[NrecoMuon]/F");
00144   HltTree->Branch("recoMuonType",           muontype       ,   "recoMuonType[NrecoMuon]/I");
00145   HltTree->Branch("recoMuonNValidTrkHits",  muonNValidTrkHits, "recoMuonNValidTrkHits[NrecoMuon]/I");
00146   HltTree->Branch("recoMuonNValidMuonHits", muonNValidMuonHits,"recoMuonNValidMuonHits[NrecoMuon]/I");
00147 
00148   HltTree->Branch("NohMuL2",&nmu2cand,"NohMuL2/I");
00149   HltTree->Branch("ohMuL2Pt",muonl2pt,"ohMuL2Pt[NohMuL2]/F");
00150   HltTree->Branch("ohMuL2Phi",muonl2phi,"ohMuL2Phi[NohMuL2]/F");
00151   HltTree->Branch("ohMuL2Eta",muonl2eta,"ohMuL2Eta[NohMuL2]/F");
00152   HltTree->Branch("ohMuL2Chg",muonl2chg,"ohMuL2Chg[NohMuL2]/I");
00153   HltTree->Branch("ohMuL2PtErr",muonl2pterr,"ohMuL2PtErr[NohMuL2]/F");
00154   HltTree->Branch("ohMuL2Iso",muonl2iso,"ohMuL2Iso[NohMuL2]/I");
00155   HltTree->Branch("ohMuL2Dr",muonl2dr,"ohMuL2Dr[NohMuL2]/F");
00156   HltTree->Branch("ohMuL2Dz",muonl2dz,"ohMuL2Dz[NohMuL2]/F");
00157   HltTree->Branch("ohMuL2VtxZ",muonl2vtxz,"ohMuL2VtxZ[NohMuL2]/F");
00158   HltTree->Branch("ohMuL2Nhits",muonl2nhits,"ohMuL2Nhits[NohMuL2]/I");
00159   HltTree->Branch("ohMuL2Nchambers",muonl2nchambers,"ohMuL2Nchambers[NohMuL2]/I");   
00160   HltTree->Branch("ohMuL2Nstat",muonl2nstat,"ohMuL2Nstat[NohMuL2]/I");   
00161   HltTree->Branch("ohMuL2L1idx",muonl21idx,"ohMuL2L1idx[NohMuL2]/I");   
00162   HltTree->Branch("NohMuL3",&nmu3cand,"NohMuL3/I");
00163   HltTree->Branch("ohMuL3Pt",muonl3pt,"ohMuL3Pt[NohMuL3]/F");
00164   HltTree->Branch("ohMuL3Phi",muonl3phi,"ohMuL3Phi[NohMuL3]/F");
00165   HltTree->Branch("ohMuL3Eta",muonl3eta,"ohMuL3Eta[NohMuL3]/F");
00166   HltTree->Branch("ohMuL3Chg",muonl3chg,"ohMuL3Chg[NohMuL3]/I");
00167   HltTree->Branch("ohMuL3PtErr",muonl3pterr,"ohMuL3PtErr[NohMuL3]/F");
00168   HltTree->Branch("ohMuL3Iso",muonl3iso,"ohMuL3Iso[NohMuL3]/I");
00169   HltTree->Branch("ohMuL3Trk10Iso",muonl3trk10iso,"ohMuL3Trk10Iso[NohMuL3]/I");
00170   HltTree->Branch("ohMuL3Dr",muonl3dr,"ohMuL3Dr[NohMuL3]/F");
00171   HltTree->Branch("ohMuL3Dz",muonl3dz,"ohMuL3Dz[NohMuL3]/F");
00172   HltTree->Branch("ohMuL3VtxZ",muonl3vtxz,"ohMuL3VtxZ[NohMuL3]/F");
00173   HltTree->Branch("ohMuL3Nhits",muonl3nhits,"ohMuL3Nhits[NohMuL3]/I");    
00174   HltTree->Branch("ohMuL3NormChi2", muonl3normchi2, "ohMuL3NormChi2[NohMuL3]/F");
00175   HltTree->Branch("ohMuL3Ntrackerhits", muonl3ntrackerhits, "ohMuL3Ntrackerhits[NohMuL3]/I"); 
00176   HltTree->Branch("ohMuL3Nmuonhits", muonl3nmuonhits, "ohMuL3Nmuonhits[NohMuL3]/I"); 
00177   HltTree->Branch("ohMuL3L2idx",muonl32idx,"ohMuL3L2idx[NohMuL3]/I");
00178   HltTree->Branch("NohOniaPixel",&nOniaPixelCand,"NohOniaPixel/I");
00179   HltTree->Branch("ohOniaPixelPt",oniaPixelpt,"ohOniaPixelPt[NohOniaPixel]/F");
00180   HltTree->Branch("ohOniaPixelPhi",oniaPixelphi,"ohOniaPixelPhi[NohOniaPixel]/F");
00181   HltTree->Branch("ohOniaPixelEta",oniaPixeleta,"ohOniaPixelEta[NohOniaPixel]/F");
00182   HltTree->Branch("ohOniaPixelChg",oniaPixelchg,"ohOniaPixelChg[NohOniaPixel]/I");
00183   HltTree->Branch("ohOniaPixelDr",oniaPixeldr,"ohOniaPixelDr[NohOniaPixel]/F");
00184   HltTree->Branch("ohOniaPixelDz",oniaPixeldz,"ohOniaPixelDz[NohOniaPixel]/F");
00185   HltTree->Branch("ohOniaPixelHits",oniaPixelHits,"ohOniaPixelHits[NohOniaPixel]/I");
00186   HltTree->Branch("ohOniaPixelNormChi2",oniaPixelNormChi2,"ohOniaPixelNormChi2[NohOniaPixel]/F");
00187   HltTree->Branch("NohOniaTrack",&nOniaTrackCand,"NohOniaTrack/I");
00188   HltTree->Branch("ohOniaTrackPt",oniaTrackpt,"ohOniaTrackPt[NohOniaTrack]/F");
00189   HltTree->Branch("ohOniaTrackPhi",oniaTrackphi,"ohOniaTrackPhi[NohOniaTrack]/F");
00190   HltTree->Branch("ohOniaTrackEta",oniaTracketa,"ohOniaTrackEta[NohOniaTrack]/F");
00191   HltTree->Branch("ohOniaTrackChg",oniaTrackchg,"ohOniaTrackChg[NohOniaTrack]/I");
00192   HltTree->Branch("ohOniaTrackDr",oniaTrackdr,"ohOniaTrackDr[NohOniaTrack]/F");
00193   HltTree->Branch("ohOniaTrackDz",oniaTrackdz,"ohOniaTrackDz[NohOniaTrack]/F");
00194   HltTree->Branch("ohOniaTrackHits",oniaTrackHits,"ohOniaTrackHits[NohOniaTrack]/I");
00195   HltTree->Branch("ohOniaTrackNormChi2",oniaTrackNormChi2,"ohOniaTrackNormChi2[NohOniaTrack]/F");
00196   HltTree->Branch("NohMuL2NoVtx",&nmu2cand,"NohMuL2NoVtx/I"); 
00197   HltTree->Branch("ohMuL2NoVtxPt",muonl2novtxpt,"ohMuL2NoVtxPt[NohMuL2NoVtx]/F"); 
00198   HltTree->Branch("ohMuL2NoVtxPhi",muonl2novtxphi,"ohMuL2NoVtxPhi[NohMuL2NoVtx]/F"); 
00199   HltTree->Branch("ohMuL2NoVtxEta",muonl2novtxeta,"ohMuL2NoVtxEta[NohMuL2NoVtx]/F"); 
00200   HltTree->Branch("ohMuL2NoVtxChg",muonl2novtxchg,"ohMuL2NoVtxChg[NohMuL2NoVtx]/I"); 
00201   HltTree->Branch("ohMuL2NoVtxPtErr",muonl2novtxpterr,"ohMuL2NoVtxPtErr[NohMuL2NoVtx]/F"); 
00202   HltTree->Branch("ohMuL2NoVtxDr",muonl2novtxdr,"ohMuL2NoVtxDr[NohMuL2NoVtx]/F"); 
00203   HltTree->Branch("ohMuL2NoVtxDz",muonl2novtxdz,"ohMuL2NoVtxDz[NohMuL2NoVtx]/F"); 
00204   HltTree->Branch("ohMuL2NoVtxNhits",muonl2novtxnhits,"ohMuL2NoVtxNhits[NohMuL2NoVtx]/I");
00205   HltTree->Branch("ohMuL2NoVtxNchambers",muonl2novtxnchambers,"ohMuL2NoVtxNchambers[NohMuL2NoVtx]/I");  
00206   HltTree->Branch("ohMuL2NoVtxL1idx",muonl2novtx1idx,"ohMuL2NoVtxL1idx[NohMuL2NoVtx]/I");   
00207   HltTree->Branch("NohDiMu",&nDiMu,"NohDiMu/I");    
00208   HltTree->Branch("ohDiMuDCA",dimudca,"ohDiMuDCA[NohDiMu]/F");    
00209   HltTree->Branch("ohDiMu1st",dimu1st,"ohDiMu1st[NohDiMu]/I");    
00210   HltTree->Branch("ohDiMu2nd",dimu2nd,"ohDiMu2nd[NohDiMu]/I");    
00211   HltTree->Branch("NohDiMuVtx",&nDiMuVtx,"NohDiMuVtx/I");    
00212   HltTree->Branch("ohDiMuVtx1st",dimuvtx1st,"ohDiMuVtx1st[NohDiMuVtx]/I");    
00213   HltTree->Branch("ohDiMuVtx2nd",dimuvtx2nd,"ohDiMuVtx2nd[NohDiMuVtx]/I");    
00214   HltTree->Branch("ohDiMuVtxChi2",dimuvtxchi2,"ohDiMuVtxChi2[NohDiMuVtx]/F");    
00215   HltTree->Branch("ohDiMuVtxR",dimuvtxr,"ohDiMuVtxR[NohDiMuVtx]/F");    
00216   HltTree->Branch("ohDiMuVtxRSig",dimuvtxrsig,"ohDiMuVtxRSig[NohDiMuVtx]/F");    
00217   HltTree->Branch("ohDiMuVtxROverSig",dimuvtxroversig,"ohDiMuVtxROverSig[NohDiMuVtx]/F");    
00218   HltTree->Branch("ohDiMuVtxCosAlpha",dimuvtxcosalpha,"ohDiMuVtxCosAlpha[NohDiMuVtx]/F");    
00219   HltTree->Branch("ohDiMuVtxMu2DIpMax",dimuvtxmu2dipmax,"ohDiMuVtxMu2DIpMax[NohDiMuVtx]/F");    
00220   HltTree->Branch("ohDiMuVtxMu2DIpMin",dimuvtxmu2dipmin,"ohDiMuVtxMu2DIpMin[NohDiMuVtx]/F");    
00221   HltTree->Branch("ohDiMuVtxMu2DIpSigMax",dimuvtxmu2dipsigmax,"ohDiMuVtxMu2DIpSigMax[NohDiMuVtx]/F");    
00222   HltTree->Branch("ohDiMuVtxMu2DIpSigMin",dimuvtxmu2dipsigmin,"ohDiMuVtxMu2DIpSigMin[NohDiMuVtx]/F");    
00223   HltTree->Branch("NohTrackerMuon",&ntrackermuoncand,"NohTrackerMuon/I");
00224   HltTree->Branch("ohTrackerMuonPt",trackermuonpt,"ohTrackerMuonPt[NohTrackerMuon]/F");
00225   HltTree->Branch("ohTrackerMuonPhi",trackermuonphi,"ohTrackerMuonPhi[NohTrackerMuon]/F");
00226   HltTree->Branch("ohTrackerMuonEta",trackermuoneta,"ohTrackerMuonEta[NohTrackerMuon]/F");
00227   HltTree->Branch("ohTrackerMuonChg",trackermuonchg,"ohTrackerMuonChg[NohTrackerMuon]/I");
00228   HltTree->Branch("ohTrackerMuonNhits",trackermuonnhits,"ohTrackerMuonNhits[NohTrackerMuon]/I");
00229 }
00230 
00231 /* **Analyze the event** */
00232 void HLTMuon::analyze(const edm::Handle<reco::MuonCollection>                 & Muon,
00233                       const edm::Handle<l1extra::L1MuonParticleCollection>    & MuCands1, 
00234                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands2,
00235                       const edm::Handle<edm::ValueMap<bool> >                 & isoMap2,
00236                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands3,
00237                       const edm::Handle<edm::ValueMap<bool> >                 & isoMap3,
00238                       const edm::Handle<edm::ValueMap<bool> >                 & isoTrk10Map3,
00239                       const edm::Handle<reco::RecoChargedCandidateCollection> & oniaPixelCands,
00240                       const edm::Handle<reco::RecoChargedCandidateCollection> & oniaTrackCands,
00241                       const edm::Handle<reco::VertexCollection> & DiMuVtxCands3,
00242                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuNoVtxCands2, 
00243                       const edm::Handle<reco::MuonCollection>                 & trkmucands,
00244                       const edm::ESHandle<MagneticField> & theMagField,
00245                       const edm::Handle<reco::BeamSpot> & recoBeamSpotHandle,
00246                       TTree* HltTree) {
00247 
00248   reco::BeamSpot::Point BSPosition(0,0,0);
00249   BSPosition = recoBeamSpotHandle->position();
00250   const GlobalPoint theBeamSpot = GlobalPoint(recoBeamSpotHandle->position().x(),
00251                                               recoBeamSpotHandle->position().y(),
00252                                               recoBeamSpotHandle->position().z());
00253   reco::BeamSpot vtxBS = *recoBeamSpotHandle;
00254 
00255   //std::cout << " Beginning HLTMuon " << std::endl;
00256 
00257   if (Muon.isValid()) {
00258     reco::MuonCollection mymuons;
00259     mymuons = * Muon;
00260     std::sort(mymuons.begin(),mymuons.end(),PtGreater());
00261     nmuon = mymuons.size();
00262     typedef reco::MuonCollection::const_iterator muiter;
00263     int imu=0;
00264     for (muiter i=mymuons.begin(); i!=mymuons.end(); i++) 
00265       {
00266         muonpt[imu]         = i->pt();
00267         muonphi[imu]        = i->phi();
00268         muoneta[imu]        = i->eta();
00269         muonet[imu]         = i->et();
00270         muone[imu]          = i->energy(); 
00271         muontype[imu]       = i->type();
00272         muoncharge[imu]     = i->charge(); 
00273         muonTrkIsoR03[imu]  = i->isolationR03().sumPt;
00274         muonECalIsoR03[imu] = i->isolationR03().emEt;
00275         muonHCalIsoR03[imu] = i->isolationR03().hadEt;
00276 
00277 
00278         if (i->globalTrack().isNonnull())
00279           {
00280             muonchi2NDF[imu] = i->globalTrack()->normalizedChi2();
00281             muonD0[imu] = i->globalTrack()->dxy(BSPosition);
00282           }
00283         else 
00284           {
00285             muonchi2NDF[imu] = -99.;
00286             muonD0[imu] = -99.;}
00287 
00288         if (i->innerTrack().isNonnull()) muonNValidTrkHits[imu] = i->innerTrack()->numberOfValidHits();
00289         else muonNValidTrkHits[imu] = -99;
00290 
00291         if (i->isGlobalMuon()!=0) muonNValidMuonHits[imu] = i->globalTrack()->hitPattern().numberOfValidMuonHits();
00292         else muonNValidMuonHits[imu] = -99;
00293 
00294         imu++;
00295       }
00296   }
00297   else {nmuon = 0;}
00298 
00299   l1extra::L1MuonParticleCollection myMucands1; 
00300   myMucands1 = * MuCands1; 
00301   //  reco::RecoChargedCandidateCollection myMucands1;
00302   std::sort(myMucands1.begin(),myMucands1.end(),PtGreater()); 
00303 
00305 
00306   // Dealing with L2 muons
00307   reco::RecoChargedCandidateCollection myMucands2;
00308   if (MuCands2.isValid()) {
00309     //     reco::RecoChargedCandidateCollection myMucands2;
00310     myMucands2 = * MuCands2;
00311     std::sort(myMucands2.begin(),myMucands2.end(),PtGreater());
00312     nmu2cand = myMucands2.size();
00313     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00314     int imu2c=0;
00315     for (cand i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00316       reco::TrackRef tk = i->get<reco::TrackRef>();
00317 
00318       muonl2pt[imu2c] = tk->pt();
00319       // eta (we require |eta|<2.5 in all filters
00320       muonl2eta[imu2c] = tk->eta();
00321       muonl2phi[imu2c] = tk->phi();
00322 
00323       // Dr (transverse distance to (0,0,0))
00324       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00325       // However, we use |dr|<200 microns at L3, which it probably too tough for LHC startup
00326       muonl2dr[imu2c] = fabs(tk->dxy(BSPosition));
00327 
00328       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00329       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00330       muonl2dz[imu2c] = tk->dz(BSPosition);
00331       muonl2vtxz[imu2c] = tk->dz();
00332       muonl2nhits[imu2c] = tk->numberOfValidHits();
00333       muonl2nchambers[imu2c] = validChambers(tk);
00334       muonl2nstat[imu2c] = tk->hitPattern().muonStationsWithAnyHits();
00335 
00336       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00337       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00338       // Baseline cuts (HLT exercise):
00339       //                Relaxed Single muon:  ptLx>16 GeV
00340       //                Isolated Single muon: ptLx>11 GeV
00341       //                Relaxed Double muon: ptLx>3 GeV
00342       double l2_err0 = tk->error(0); // error on q/p
00343       double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00344       //       double ptLx = tk->pt();
00345       // convert 50% efficiency threshold to 90% efficiency threshold
00346       // For L2 muons: nsigma_Pt_ = 3.9
00347       //       double nsigma_Pt_ = 3.9;
00348       // For L3 muons: nsigma_Pt_ = 2.2
00349       // these are the old TDR values for nsigma_Pt_
00350       // We know that these values are slightly smaller for CMSSW
00351       // But as quoted above, we want to get rid of this gymnastics in the future
00352       //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00353 
00354       // Charge
00355       // We use the charge in some dimuon paths
00356       muonl2pterr[imu2c] = l2_err0/l2_abspar0;
00357       muonl2chg[imu2c] = tk->charge();
00358 
00359       if (isoMap2.isValid()){
00360         // Isolation flag (this is a bool value: true => isolated)
00361         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap2)[tk];
00362         muonl2iso[imu2c] = muon1IsIsolated;
00363       }
00364       else {muonl2iso[imu2c] = -999;}
00365 
00366       //JH
00367       l1extra::L1MuonParticleRef l1; 
00368       int il2 = 0; 
00369       //find the corresponding L1 
00370       l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00371       il2++; 
00372       int imu1idx = 0; 
00373       if (MuCands1.isValid()) { 
00374         typedef l1extra::L1MuonParticleCollection::const_iterator candl1; 
00375         for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) { 
00376           if((j->pt() == l1->pt()) &&
00377              (j->eta() == l1->eta()) &&
00378              (j->phi() == l1->phi()) &&
00379              (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00380             {break;}
00381           //      std::cout << << std::endl;
00382           //          if ( tkl1 == l1 ) {break;} 
00383           imu1idx++; 
00384         } 
00385       } 
00386       else {imu1idx = -999;} 
00387       muonl21idx[imu2c] = imu1idx; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00388       //end JH
00389 
00390       imu2c++;
00391     }
00392   }
00393   else {nmu2cand = 0;}
00394 
00395   // Dealing with L3 muons
00396   reco::RecoChargedCandidateCollection myMucands3;
00397   if (MuCands3.isValid()) {
00398     myMucands3 = * MuCands3;
00399     std::sort(myMucands3.begin(),myMucands3.end(),PtGreater());
00400     nmu3cand = myMucands3.size();
00401     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00402     int imu3c=0;
00403     int idimuc=0;
00404     for (cand i=myMucands3.begin(); i!=myMucands3.end(); i++) {
00405       reco::TrackRef tk = i->get<reco::TrackRef>();
00406 
00407       reco::TrackRef staTrack;
00408       typedef reco::MuonTrackLinksCollection::const_iterator l3muon;
00409       int il3 = 0;
00410       //find the corresponding L2 track
00411       staTrack = tk->seedRef().castTo<edm::Ref< L3MuonTrajectorySeedCollection> >()->l2Track();
00412       il3++;
00413       int imu2idx = 0;
00414       if (MuCands2.isValid()) {
00415         typedef reco::RecoChargedCandidateCollection::const_iterator candl2;
00416         for (candl2 i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00417           reco::TrackRef tkl2 = i->get<reco::TrackRef>();
00418           if ( tkl2 == staTrack ) {break;}
00419           imu2idx++;
00420         }
00421       }
00422       else {imu2idx = -999;}
00423       muonl32idx[imu3c] = imu2idx; // Index of the L2 muon having matched with the L3 muon with index imu3c
00424 
00425       muonl3pt[imu3c] = tk->pt();
00426       // eta (we require |eta|<2.5 in all filters
00427       muonl3eta[imu3c] = tk->eta();
00428       muonl3phi[imu3c] = tk->phi();
00429 
00430       //       // Dr (transverse distance to (0,0,0))
00431       //       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00432       //       // However, we use |dr|<300 microns at L3, which it probably too tough for LHC startup
00433       muonl3dr[imu3c] = fabs(tk->dxy(BSPosition));
00434 
00435       //       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00436       //       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00437       muonl3dz[imu3c] = tk->dz(BSPosition);
00438       muonl3vtxz[imu3c] = tk->dz();
00439       muonl3nhits[imu3c] = tk->numberOfValidHits();  
00440 
00441       //       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00442       //       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00443       //       // Baseline cuts (HLT exercise):
00444       //       //                Relaxed Single muon:  ptLx>16 GeV
00445       //       //                Isolated Single muon: ptLx>11 GeV
00446       //       //                Relaxed Double muon: ptLx>3 GeV
00447       double l3_err0 = tk->error(0); // error on q/p
00448       double l3_abspar0 = fabs(tk->parameter(0)); // |q/p|
00449       // //       double ptLx = tk->pt();
00450       //       // convert 50% efficiency threshold to 90% efficiency threshold
00451       //       // For L2 muons: nsigma_Pt_ = 3.9
00452       //       // For L3 muons: nsigma_Pt_ = 2.2
00453       // //       double nsigma_Pt_ = 2.2;
00454       //       // these are the old TDR values for nsigma_Pt_
00455       //       // We know that these values are slightly smaller for CMSSW
00456       //       // But as quoted above, we want to get rid of this gymnastics in the future
00457       // //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00458 
00459       // Charge
00460       // We use the charge in some dimuon paths
00461       muonl3pterr[imu3c] = l3_err0/l3_abspar0;
00462       muonl3chg[imu3c] = tk->charge();
00463 
00464       muonl3normchi2[imu3c] = tk->normalizedChi2();
00465       muonl3ntrackerhits[imu3c] = tk->hitPattern().numberOfValidTrackerHits();
00466       muonl3nmuonhits[imu3c] = tk->hitPattern().numberOfValidMuonHits();
00467 
00468       if (isoMap3.isValid()){
00469         // Isolation flag (this is a bool value: true => isolated)
00470         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap3)[tk];
00471         muonl3iso[imu3c] = muon1IsIsolated;
00472       }
00473       else {muonl3iso[imu3c] = -999;}
00474 
00475       if (isoTrk10Map3.isValid()){
00476         // Isolation flag (this is a bool value: true => isolated) 
00477         edm::ValueMap<bool> ::value_type muon1IsTrk10Isolated = (*isoTrk10Map3)[tk];
00478         muonl3trk10iso[imu3c] = muon1IsTrk10Isolated;
00479       }
00480       else {muonl3trk10iso[imu3c] =-999;}
00481 
00482       //Check DCA for muon combinations
00483       int imu3c2nd = imu3c + 1;// This will be the index in the hltTree for the 2nd muon of the dimuon combination
00484 
00485       for (cand j=i; j!=myMucands3.end(); j++) if (i!=j) {//Loop over all L3 muons from the one we are already treating
00486         reco::TrackRef tk2nd = j->get<reco::TrackRef>();
00487         reco::TransientTrack transMu1(*tk, &(*theMagField) );
00488         reco::TransientTrack transMu2(*tk2nd, &(*theMagField) );
00489         TrajectoryStateClosestToPoint mu1TS = transMu1.impactPointTSCP();
00490         TrajectoryStateClosestToPoint mu2TS = transMu2.impactPointTSCP();
00491         if (mu1TS.isValid() && mu2TS.isValid()) {
00492           ClosestApproachInRPhi cApp;
00493           cApp.calculate(mu1TS.theState(), mu2TS.theState());
00494           if (cApp.status()) {
00495             dimudca[idimuc] = cApp.distance();//Save the DCA
00496             dimu1st[idimuc] = imu3c;//Save which is the index in the hltTree for the 1st muon
00497             dimu2nd[idimuc] = imu3c2nd;//Save which is the index in the hltTree for the 2nd muon
00498             idimuc++;
00499           }
00500         }
00501         imu3c2nd++;
00502       }
00503 
00504       imu3c++;
00505 
00506     }
00507     nDiMu = idimuc;
00508   }
00509 
00510   else {nmu3cand = 0;  nDiMu = 0;}
00511   // Dealing with dimu vertices
00512   reco::VertexCollection myDimuvtxcands3;
00513   if (DiMuVtxCands3.isValid()) {
00514     myDimuvtxcands3 = * DiMuVtxCands3;
00515     nDiMuVtx = myDimuvtxcands3.size();
00516     typedef reco::VertexCollection::const_iterator cand;
00517     int idimu3c=0;
00518     for (cand ivtx = myDimuvtxcands3.begin(); ivtx != myDimuvtxcands3.end(); ++ivtx) {
00519       dimuvtxchi2[idimu3c] = ivtx->normalizedChi2();
00520       reco::Vertex::trackRef_iterator trackIt = ivtx->tracks_begin();
00521       reco::TrackRef vertextkRef1 = (*trackIt).castTo<reco::TrackRef>();
00522       ++trackIt;
00523       reco::TrackRef vertextkRef2 = (*trackIt).castTo<reco::TrackRef>();
00524       dimuvtx2nd[idimu3c] = -1; dimuvtx1st[idimu3c] = -1;
00525       for (int j=0 ; j<nmu3cand ; j++){
00526         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; 
00527         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; 
00528       }
00529       math::XYZVector pperp(vertextkRef1->px() + vertextkRef2->px(), 
00530                             vertextkRef1->py() + vertextkRef2->py(), 
00531                             0.);
00532       reco::Vertex::Point vpoint = ivtx->position();
00533       GlobalPoint vtxPos (vpoint.x(), vpoint.y(), vpoint.z());
00534       reco::Vertex::Error verr = ivtx->error();
00535       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));
00536       GlobalPoint vtxDisFromBS(-1*((vtxBS.x0() - vtxPos.x()) + (vtxPos.z() - vtxBS.z0())*vtxBS.dxdz()),
00537                                -1*((vtxBS.y0() - vtxPos.y()) + (vtxPos.z() - vtxBS.z0())*vtxBS.dydz()), 0.0);
00538       dimuvtxr[idimu3c] = vtxDisFromBS.perp();
00539       dimuvtxrsig[idimu3c] = sqrt(vtxErr.rerr(vtxDisFromBS));
00540       dimuvtxroversig[idimu3c] = dimuvtxr[idimu3c]/dimuvtxrsig[idimu3c];
00541       reco::Vertex::Point vperp(vtxDisFromBS.x(),vtxDisFromBS.y(),0.);
00542       dimuvtxcosalpha[idimu3c] = vperp.Dot(pperp)/(vperp.R()*pperp.R());
00543       float mu1ip = -1.0;
00544       float mu2ip = -1.0;
00545       float mu1ipsig = -1.0;
00546       float mu2ipsig = -1.0;
00547       reco::TransientTrack transMu1(*vertextkRef1, &(*theMagField) );
00548       TrajectoryStateClosestToPoint trajMu1BS = transMu1.trajectoryStateClosestToPoint(theBeamSpot);
00549       if(trajMu1BS.isValid()){
00550         mu1ip = fabs(trajMu1BS.perigeeParameters().transverseImpactParameter());
00551         if(trajMu1BS.hasError()) mu1ipsig = mu1ip/trajMu1BS.perigeeError().transverseImpactParameterError();
00552       }
00553       reco::TransientTrack transMu2(*vertextkRef2, &(*theMagField) );
00554       TrajectoryStateClosestToPoint trajMu2BS = transMu2.trajectoryStateClosestToPoint(theBeamSpot);
00555       if(trajMu2BS.isValid()){
00556         mu2ip = fabs(trajMu2BS.perigeeParameters().transverseImpactParameter());
00557         if(trajMu2BS.hasError()) mu2ipsig = mu2ip/trajMu2BS.perigeeError().transverseImpactParameterError();
00558       }
00559       dimuvtxmu2dipmax[idimu3c] = fmax(mu1ip,mu2ip);
00560       dimuvtxmu2dipmin[idimu3c] = fmin(mu1ip,mu2ip);
00561       dimuvtxmu2dipsigmax[idimu3c] = fmax(mu1ipsig,mu2ipsig);
00562       dimuvtxmu2dipsigmin[idimu3c] = fmin(mu1ipsig,mu2ipsig);
00563     }
00564 
00565 
00566   }
00567   else {nDiMuVtx = 0;}
00568 
00569 
00570   // Dealing with L2 no-Vertex muons
00571   reco::RecoChargedCandidateCollection muNoVtxMucands2;
00572   if (MuNoVtxCands2.isValid()) {
00573     muNoVtxMucands2 = * MuNoVtxCands2;
00574     std::sort(muNoVtxMucands2.begin(),muNoVtxMucands2.end(),PtGreater());
00575     nmu2cand = muNoVtxMucands2.size();
00576     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00577     int imu2c=0;
00578     for (cand i=muNoVtxMucands2.begin(); i!=muNoVtxMucands2.end(); i++) {
00579       reco::TrackRef tk = i->get<reco::TrackRef>();
00580 
00581       muonl2novtxpt[imu2c] = tk->pt();
00582       muonl2novtxeta[imu2c] = tk->eta();
00583       muonl2novtxphi[imu2c] = tk->phi();
00584       muonl2novtxdr[imu2c] = fabs(tk->dxy(BSPosition));
00585       muonl2novtxdz[imu2c] = tk->dz(BSPosition);
00586       muonl2novtxnhits[imu2c] = tk->numberOfValidHits();
00587       muonl2novtxnchambers[imu2c] = validChambers(tk);
00588 
00589       double l2_err0 = tk->error(0); // error on q/p
00590       double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00591 
00592       muonl2novtxpterr[imu2c] = l2_err0/l2_abspar0;
00593       muonl2novtxchg[imu2c] = tk->charge();
00594 
00595       l1extra::L1MuonParticleRef l1; 
00596       int il2 = 0; 
00597       //find the corresponding L1 
00598       l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00599       il2++; 
00600       int imu1idx = 0; 
00601       if (MuCands1.isValid()) { 
00602         typedef l1extra::L1MuonParticleCollection::const_iterator candl1; 
00603         for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) { 
00604           if((j->pt() == l1->pt()) &&
00605              (j->eta() == l1->eta()) &&
00606              (j->phi() == l1->phi()) &&
00607              (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00608             {break;}
00609           imu1idx++; 
00610         } 
00611       } 
00612       else {imu1idx = -999;} 
00613       muonl2novtx1idx[imu2c] = imu1idx; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00614 
00615       imu2c++;
00616     }
00617   }
00618   else {nmu2cand = 0;}
00619 
00620 
00621 
00622   // Dealing with Onia Pixel tracks
00623   reco::RecoChargedCandidateCollection myOniaPixelCands;
00624   if (oniaPixelCands.isValid()) {
00625     myOniaPixelCands = * oniaPixelCands;
00626     std::sort(myOniaPixelCands.begin(),myOniaPixelCands.end(),PtGreater());
00627     nOniaPixelCand = myOniaPixelCands.size();
00628     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00629     int ic=0;
00630     for (cand i=myOniaPixelCands.begin(); i!=myOniaPixelCands.end(); i++) {
00631       reco::TrackRef tk = i->get<reco::TrackRef>();
00632 
00633       oniaPixelpt[ic] = tk->pt();
00634       oniaPixeleta[ic] = tk->eta();
00635       oniaPixelphi[ic] = tk->phi();
00636       oniaPixeldr[ic] = tk->dxy(BSPosition);
00637       oniaPixeldz[ic] = tk->dz(BSPosition);
00638       oniaPixelchg[ic] = tk->charge();
00639       oniaPixelHits[ic] = tk->numberOfValidHits();
00640       oniaPixelNormChi2[ic] = tk->normalizedChi2();
00641 
00642       ic++;
00643     }
00644   }
00645   else {nOniaPixelCand = 0;}
00646 
00647   // Dealing with Onia Tracks
00648   reco::RecoChargedCandidateCollection myOniaTrackCands;
00649   if (oniaTrackCands.isValid()) {
00650     myOniaTrackCands = * oniaTrackCands;
00651     std::sort(myOniaTrackCands.begin(),myOniaTrackCands.end(),PtGreater());
00652     nOniaTrackCand = myOniaTrackCands.size();
00653     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00654     int ic=0;
00655     for (cand i=myOniaTrackCands.begin(); i!=myOniaTrackCands.end(); i++) {
00656       reco::TrackRef tk = i->get<reco::TrackRef>();
00657 
00658       oniaTrackpt[ic] = tk->pt();
00659       oniaTracketa[ic] = tk->eta();
00660       oniaTrackphi[ic] = tk->phi();
00661       oniaTrackdr[ic] = tk->dxy(BSPosition);
00662       oniaTrackdz[ic] = tk->dz(BSPosition);
00663       oniaTrackchg[ic] = tk->charge();
00664       oniaTrackHits[ic] = tk->numberOfValidHits();
00665       oniaTrackNormChi2[ic] = tk->normalizedChi2();
00666 
00667       ic++;
00668     }
00669   }
00670   else {nOniaTrackCand = 0;}
00671 
00672 
00673   // Dealing with trackerMuons
00674   if(trkmucands.isValid()) {
00675     int itrackermuc=0;
00676     for ( unsigned int i=0; i<trkmucands->size(); ++i ){
00677       const reco::Muon& muon(trkmucands->at(i));
00678       if (muon.isTrackerMuon()) {
00679         trackermuonpt[itrackermuc] = muon.pt();
00680         trackermuoneta[itrackermuc] = muon.eta();
00681         trackermuonphi[itrackermuc] = muon.phi();
00682         trackermuonchg[itrackermuc] = muon.charge();
00683         if ( !muon.innerTrack().isNull() ){
00684           trackermuonnhits[itrackermuc] = muon.innerTrack()->numberOfValidHits();
00685         }
00686         itrackermuc++;
00687       }
00688     }
00689     ntrackermuoncand=itrackermuc;
00690   }
00691   else {ntrackermuoncand = 0;}
00692 
00694 }
00695 
00696 int HLTMuon::validChambers(const reco::TrackRef & track)
00697 {
00698   // count hits in chambers using std::maps
00699   std::map<uint32_t,int> DTchambers;
00700   std::map<uint32_t,int> CSCchambers;
00701 
00702   for (trackingRecHit_iterator hit = track->recHitsBegin();  hit != track->recHitsEnd();  ++hit) {
00703     if( !((*hit)->isValid()) ) continue;
00704 
00705     DetId id = (*hit)->geographicalId();
00706 
00707     if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::DT) {
00708       // get the DT chamber index, not the layer index, by using DTChamberId
00709       uint32_t index = DTChamberId(id).rawId();
00710 
00711       if (DTchambers.find(index) == DTchambers.end()) {
00712         DTchambers[index] = 0;
00713       }
00714       DTchambers[index]++;
00715     }
00716 
00717     else if (id.det() == DetId::Muon  &&  id.subdetId() == MuonSubdetId::CSC) {
00718       // get the CSC chamber index, not the layer index, by explicitly setting the layer id to 0
00719       CSCDetId id2(id);
00720       uint32_t index = CSCDetId(id2.endcap(), id2.station(), id2.ring(), id2.chamber(), 0);
00721 
00722       if (CSCchambers.find(index) == CSCchambers.end()) {
00723         CSCchambers[index] = 0;
00724       }
00725       CSCchambers[index]++;
00726     }
00727   }
00728 
00729   // count chambers that satisfy minimal numbers of hits per chamber
00730   int validChambers = 0;
00731 
00732   int minDThits = 1;
00733   int minCSChits = 1;
00734 
00735   for (std::map<uint32_t,int>::const_iterator iter = DTchambers.begin();  iter != DTchambers.end();  ++iter) {
00736     if (iter->second >= minDThits) {
00737       validChambers++;
00738     }
00739   }
00740   for (std::map<uint32_t,int>::const_iterator iter = CSCchambers.begin();  iter != CSCchambers.end();  ++iter) {
00741     if (iter->second >= minCSChits) {
00742       validChambers++;
00743     }
00744   }
00745   return validChambers;
00746 }
/data/refman/pasoursint/CMSSW_4_4_5_patch3/src/HLTrigger/HLTanalyzers/src/HLTMuon.cc
