CMS 3D CMS Logo

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