CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/HLTrigger/HLTanalyzers/src/HLTMuon.cc

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   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   //pf
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   // Muon-specific branches of the tree 
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 /* **Analyze the event** */
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   //std::cout << " Beginning HLTMuon " << std::endl;
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   //  reco::RecoChargedCandidateCollection myMucands1;
00338   std::sort(myMucands1.begin(),myMucands1.end(),PtGreater()); 
00339 
00341 
00342   // Dealing with L2 muons
00343   reco::RecoChargedCandidateCollection myMucands2;
00344   if (MuCands2.isValid()) {
00345     //     reco::RecoChargedCandidateCollection myMucands2;
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       // eta (we require |eta|<2.5 in all filters
00356       muonl2eta[imu2c] = tk->eta();
00357       muonl2phi[imu2c] = tk->phi();
00358 
00359       // Dr (transverse distance to (0,0,0))
00360       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00361       // However, we use |dr|<200 microns at L3, which it probably too tough for LHC startup
00362       muonl2dr[imu2c] = fabs(tk->dxy(BSPosition));
00363 
00364       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00365       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
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       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00373       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00374       // Baseline cuts (HLT exercise):
00375       //                Relaxed Single muon:  ptLx>16 GeV
00376       //                Isolated Single muon: ptLx>11 GeV
00377       //                Relaxed Double muon: ptLx>3 GeV
00378       double l2_err0 = tk->error(0); // error on q/p
00379       double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00380       //       double ptLx = tk->pt();
00381       // convert 50% efficiency threshold to 90% efficiency threshold
00382       // For L2 muons: nsigma_Pt_ = 3.9
00383       //       double nsigma_Pt_ = 3.9;
00384       // For L3 muons: nsigma_Pt_ = 2.2
00385       // these are the old TDR values for nsigma_Pt_
00386       // We know that these values are slightly smaller for CMSSW
00387       // But as quoted above, we want to get rid of this gymnastics in the future
00388       //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00389 
00390       // Charge
00391       // We use the charge in some dimuon paths
00392       muonl2pterr[imu2c] = l2_err0/l2_abspar0;
00393       muonl2chg[imu2c] = tk->charge();
00394 
00395       if (isoMap2.isValid()){
00396         // Isolation flag (this is a bool value: true => isolated)
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       //find the corresponding L1 
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           //      std::cout << << std::endl;
00417           //          if ( tkl1 == l1 ) {break;} 
00418           imu1idx++; 
00419         } 
00420       } 
00421       else {imu1idx = -999;} 
00422       muonl21idx[imu2c] = imu1idx; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00423 
00424       imu2c++;
00425     }
00426   }
00427   else {nmu2cand = 0;}
00428 
00429   // Dealing with L3 muons
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       //find the corresponding L2 track
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; // Index of the L2 muon having matched with the L3 muon with index imu3c
00461       muonl32idx[imu3c] = imu2idx;
00462 
00463       muonl3globalpt[imu3c] = tk->pt();
00464       muonl3pt[imu3c] = candref->pt();
00465       // eta (we require |eta|<2.5 in all filters
00466       muonl3globaleta[imu3c] = tk->eta();
00467       muonl3globalphi[imu3c] = tk->phi();
00468       muonl3eta[imu3c] = candref->eta();
00469       muonl3phi[imu3c] = candref->phi();
00470 
00471       //       // Dr (transverse distance to (0,0,0))
00472       //       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00473       //       // However, we use |dr|<300 microns at L3, which it probably too tough for LHC startup
00474       muonl3dr[imu3c] = fabs(tk->dxy(BSPosition));
00475 
00476       //       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00477       //       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00478       muonl3dz[imu3c] = tk->dz(BSPosition);
00479       muonl3vtxz[imu3c] = tk->dz();
00480       muonl3nhits[imu3c] = tk->numberOfValidHits();  
00481 
00482       //       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00483       //       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00484       //       // Baseline cuts (HLT exercise):
00485       //       //                Relaxed Single muon:  ptLx>16 GeV
00486       //       //                Isolated Single muon: ptLx>11 GeV
00487       //       //                Relaxed Double muon: ptLx>3 GeV
00488       double l3_err0 = tk->error(0); // error on q/p
00489       double l3_abspar0 = fabs(tk->parameter(0)); // |q/p|
00490       // //       double ptLx = tk->pt();
00491       //       // convert 50% efficiency threshold to 90% efficiency threshold
00492       //       // For L2 muons: nsigma_Pt_ = 3.9
00493       //       // For L3 muons: nsigma_Pt_ = 2.2
00494       // //       double nsigma_Pt_ = 2.2;
00495       //       // these are the old TDR values for nsigma_Pt_
00496       //       // We know that these values are slightly smaller for CMSSW
00497       //       // But as quoted above, we want to get rid of this gymnastics in the future
00498       // //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00499 
00500       // Charge
00501       // We use the charge in some dimuon paths
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         // Isolation flag (this is a bool value: true => isolated)
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         // Isolation flag (this is a bool value: true => isolated) 
00519         edm::ValueMap<bool> ::value_type muon1IsTrk10Isolated = (*isoTrk10Map3)[tk];
00520         muonl3trk10iso[imu3c] = muon1IsTrk10Isolated;
00521       }
00522       else {muonl3trk10iso[imu3c] =-999;}
00523 
00524       //Check DCA for muon combinations
00525       int imu3c2nd = imu3c + 1;// This will be the index in the hltTree for the 2nd muon of the dimuon combination
00526 
00527       for (cand j=i; j!=myMucands3.end(); j++) if (i!=j) {//Loop over all L3 muons from the one we are already treating
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();//Save the DCA
00538             dimu1st[idimuc] = imu3c;//Save which is the index in the hltTree for the 1st muon
00539             dimu2nd[idimuc] = imu3c2nd;//Save which is the index in the hltTree for the 2nd muon
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   // Dealing with dimu vertices
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   // Dealing with L2 no-Vertex muons
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); // error on q/p
00632       double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00633 
00634       muonl2novtxpterr[imu2c] = l2_err0/l2_abspar0;
00635       muonl2novtxchg[imu2c] = tk->charge();
00636 
00637       l1extra::L1MuonParticleRef l1; 
00638       int il2 = 0; 
00639       //find the corresponding L1 
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; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00656 
00657       imu2c++;
00658     }
00659   }
00660   else {nmu2cand = 0;}
00661 
00662 
00663 
00664   // Dealing with Onia Pixel tracks
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   // Dealing with Onia Tracks
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   // Dealing with trackerMuons
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   // count hits in chambers using std::maps
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       // get the DT chamber index, not the layer index, by using DTChamberId
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       // get the CSC chamber index, not the layer index, by explicitly setting the layer id to 0
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   // count chambers that satisfy minimal numbers of hits per chamber
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 }