CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/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   muonl2chg = new int[kMaxMuonL2];
00056   muonl2pterr = new float[kMaxMuonL2];
00057   muonl2iso = new int[kMaxMuonL2];
00058   muonl21idx = new int[kMaxMuonL2];
00059   const int kMaxMuonL3 = 500;
00060   muonl3pt = new float[kMaxMuonL3];
00061   muonl3phi = new float[kMaxMuonL3];
00062   muonl3eta = new float[kMaxMuonL3];
00063   muonl3dr = new float[kMaxMuonL3];
00064   muonl3dz = new float[kMaxMuonL3];
00065   muonl3chg = new int[kMaxMuonL3];
00066   muonl3pterr = new float[kMaxMuonL3];
00067   muonl3iso = new int[kMaxMuonL3];
00068   muonl32idx = new int[kMaxMuonL3];
00069   const int kMaxOniaPixel = 500;
00070   oniaPixelpt = new float[kMaxOniaPixel];
00071   oniaPixelphi = new float[kMaxOniaPixel];
00072   oniaPixeleta = new float[kMaxOniaPixel];
00073   oniaPixeldr = new float[kMaxOniaPixel];
00074   oniaPixeldz = new float[kMaxOniaPixel];
00075   oniaPixelchg = new int[kMaxOniaPixel];
00076   oniaPixelHits = new int[kMaxOniaPixel];
00077   oniaPixelNormChi2 = new float[kMaxOniaPixel];
00078   const int kMaxTrackPixel = 500;
00079   oniaTrackpt = new float[kMaxTrackPixel];
00080   oniaTrackphi = new float[kMaxTrackPixel];
00081   oniaTracketa = new float[kMaxTrackPixel];
00082   oniaTrackdr = new float[kMaxTrackPixel];
00083   oniaTrackdz = new float[kMaxTrackPixel];
00084   oniaTrackchg = new int[kMaxTrackPixel];
00085   oniaTrackHits = new int[kMaxTrackPixel];
00086   oniaTrackNormChi2 = new float[kMaxTrackPixel];
00087   const int kMaxMuonL2NoVtx = 500; 
00088   muonl2novtxpt = new float[kMaxMuonL2NoVtx]; 
00089   muonl2novtxphi = new float[kMaxMuonL2NoVtx]; 
00090   muonl2novtxeta = new float[kMaxMuonL2NoVtx]; 
00091   muonl2novtxdr = new float[kMaxMuonL2NoVtx]; 
00092   muonl2novtxdz = new float[kMaxMuonL2NoVtx]; 
00093   muonl2novtxchg = new int[kMaxMuonL2NoVtx]; 
00094   muonl2novtxpterr = new float[kMaxMuonL2NoVtx]; 
00095   muonl2novtx1idx = new int[kMaxMuonL2NoVtx]; 
00096 
00097   // Muon-specific branches of the tree 
00098   HltTree->Branch("NrecoMuon",&nmuon,"NrecoMuon/I");
00099   HltTree->Branch("recoMuonPt",muonpt,"recoMuonPt[NrecoMuon]/F");
00100   HltTree->Branch("recoMuonPhi",muonphi,"recoMuonPhi[NrecoMuon]/F");
00101   HltTree->Branch("recoMuonEta",muoneta,"recoMuonEta[NrecoMuon]/F");
00102   HltTree->Branch("recoMuonEt",muonet,"recoMuonEt[NrecoMuon]/F");
00103   HltTree->Branch("recoMuonE",muone,"recoMuonE[NrecoMuon]/F");
00104   HltTree->Branch("recoMuonChi2NDF",        muonchi2NDF,       "recoMuonChi2NDF[NrecoMuon]/F");
00105   HltTree->Branch("recoMuonCharge",         muoncharge  ,      "recoMuonCharge[NrecoMuon]/F");
00106   HltTree->Branch("recoMuonTrkIsoR03",      muonTrkIsoR03 ,    "recoMuonTrkIsoR03[NrecoMuon]/F");
00107   HltTree->Branch("recoMuonECalIsoR03",     muonECalIsoR03 ,   "recoMuonECalIsoR03[NrecoMuon]/F");
00108   HltTree->Branch("recoMuonHCalIsoR03",     muonHCalIsoR03 ,   "recoMuonHCalIsoR03[NrecoMuon]/F");
00109   HltTree->Branch("recoMuonD0",             muonD0 ,           "recoMuonD0[NrecoMuon]/F");
00110   HltTree->Branch("recoMuonType",           muontype       ,   "recoMuonType[NrecoMuon]/I");
00111   HltTree->Branch("recoMuonNValidTrkHits",  muonNValidTrkHits, "recoMuonNValidTrkHits[NrecoMuon]/I");
00112   HltTree->Branch("recoMuonNValidMuonHits", muonNValidMuonHits,"recoMuonNValidMuonHits[NrecoMuon]/I");
00113   
00114   HltTree->Branch("NohMuL2",&nmu2cand,"NohMuL2/I");
00115   HltTree->Branch("ohMuL2Pt",muonl2pt,"ohMuL2Pt[NohMuL2]/F");
00116   HltTree->Branch("ohMuL2Phi",muonl2phi,"ohMuL2Phi[NohMuL2]/F");
00117   HltTree->Branch("ohMuL2Eta",muonl2eta,"ohMuL2Eta[NohMuL2]/F");
00118   HltTree->Branch("ohMuL2Chg",muonl2chg,"ohMuL2Chg[NohMuL2]/I");
00119   HltTree->Branch("ohMuL2PtErr",muonl2pterr,"ohMuL2PtErr[NohMuL2]/F");
00120   HltTree->Branch("ohMuL2Iso",muonl2iso,"ohMuL2Iso[NohMuL2]/I");
00121   HltTree->Branch("ohMuL2Dr",muonl2dr,"ohMuL2Dr[NohMuL2]/F");
00122   HltTree->Branch("ohMuL2Dz",muonl2dz,"ohMuL2Dz[NohMuL2]/F");
00123   HltTree->Branch("ohMuL2L1idx",muonl21idx,"ohMuL2L1idx[NohMuL2]/I");   
00124   HltTree->Branch("NohMuL3",&nmu3cand,"NohMuL3/I");
00125   HltTree->Branch("ohMuL3Pt",muonl3pt,"ohMuL3Pt[NohMuL3]/F");
00126   HltTree->Branch("ohMuL3Phi",muonl3phi,"ohMuL3Phi[NohMuL3]/F");
00127   HltTree->Branch("ohMuL3Eta",muonl3eta,"ohMuL3Eta[NohMuL3]/F");
00128   HltTree->Branch("ohMuL3Chg",muonl3chg,"ohMuL3Chg[NohMuL3]/I");
00129   HltTree->Branch("ohMuL3PtErr",muonl3pterr,"ohMuL3PtErr[NohMuL3]/F");
00130   HltTree->Branch("ohMuL3Iso",muonl3iso,"ohMuL3Iso[NohMuL3]/I");
00131   HltTree->Branch("ohMuL3Dr",muonl3dr,"ohMuL3Dr[NohMuL3]/F");
00132   HltTree->Branch("ohMuL3Dz",muonl3dz,"ohMuL3Dz[NohMuL3]/F");
00133   HltTree->Branch("ohMuL3L2idx",muonl32idx,"ohMuL3L2idx[NohMuL3]/I");
00134   HltTree->Branch("NohOniaPixel",&nOniaPixelCand,"NohOniaPixel/I");
00135   HltTree->Branch("ohOniaPixelPt",oniaPixelpt,"ohOniaPixelPt[NohOniaPixel]/F");
00136   HltTree->Branch("ohOniaPixelPhi",oniaPixelphi,"ohOniaPixelPhi[NohOniaPixel]/F");
00137   HltTree->Branch("ohOniaPixelEta",oniaPixeleta,"ohOniaPixelEta[NohOniaPixel]/F");
00138   HltTree->Branch("ohOniaPixelChg",oniaPixelchg,"ohOniaPixelChg[NohOniaPixel]/I");
00139   HltTree->Branch("ohOniaPixelDr",oniaPixeldr,"ohOniaPixelDr[NohOniaPixel]/F");
00140   HltTree->Branch("ohOniaPixelDz",oniaPixeldz,"ohOniaPixelDz[NohOniaPixel]/F");
00141   HltTree->Branch("ohOniaPixelHits",oniaPixelHits,"ohOniaPixelHits[NohOniaPixel]/I");
00142   HltTree->Branch("ohOniaPixelNormChi2",oniaPixelNormChi2,"ohOniaPixelNormChi2[NohOniaPixel]/F");
00143   HltTree->Branch("NohOniaTrack",&nOniaTrackCand,"NohOniaTrack/I");
00144   HltTree->Branch("ohOniaTrackPt",oniaTrackpt,"ohOniaTrackPt[NohOniaTrack]/F");
00145   HltTree->Branch("ohOniaTrackPhi",oniaTrackphi,"ohOniaTrackPhi[NohOniaTrack]/F");
00146   HltTree->Branch("ohOniaTrackEta",oniaTracketa,"ohOniaTrackEta[NohOniaTrack]/F");
00147   HltTree->Branch("ohOniaTrackChg",oniaTrackchg,"ohOniaTrackChg[NohOniaTrack]/I");
00148   HltTree->Branch("ohOniaTrackDr",oniaTrackdr,"ohOniaTrackDr[NohOniaTrack]/F");
00149   HltTree->Branch("ohOniaTrackDz",oniaTrackdz,"ohOniaTrackDz[NohOniaTrack]/F");
00150   HltTree->Branch("ohOniaTrackHits",oniaTrackHits,"ohOniaTrackHits[NohOniaTrack]/I");
00151   HltTree->Branch("ohOniaTrackNormChi2",oniaTrackNormChi2,"ohOniaTrackNormChi2[NohOniaTrack]/F");
00152   HltTree->Branch("NohMuL2NoVtx",&nmu2cand,"NohMuL2NoVtx/I"); 
00153   HltTree->Branch("ohMuL2NoVtxPt",muonl2novtxpt,"ohMuL2NoVtxPt[NohMuL2NoVtx]/F"); 
00154   HltTree->Branch("ohMuL2NoVtxPhi",muonl2novtxphi,"ohMuL2NoVtxPhi[NohMuL2NoVtx]/F"); 
00155   HltTree->Branch("ohMuL2NoVtxEta",muonl2novtxeta,"ohMuL2NoVtxEta[NohMuL2NoVtx]/F"); 
00156   HltTree->Branch("ohMuL2NoVtxChg",muonl2novtxchg,"ohMuL2NoVtxChg[NohMuL2NoVtx]/I"); 
00157   HltTree->Branch("ohMuL2NoVtxPtErr",muonl2novtxpterr,"ohMuL2NoVtxPtErr[NohMuL2NoVtx]/F"); 
00158   HltTree->Branch("ohMuL2NoVtxDr",muonl2novtxdr,"ohMuL2NoVtxDr[NohMuL2NoVtx]/F"); 
00159   HltTree->Branch("ohMuL2NoVtxDz",muonl2novtxdz,"ohMuL2NoVtxDz[NohMuL2NoVtx]/F"); 
00160   HltTree->Branch("ohMuL2NoVtxL1idx",muonl2novtx1idx,"ohMuL2NoVtxL1idx[NohMuL2NoVtx]/I");    
00161 
00162 }
00163 
00164 /* **Analyze the event** */
00165 void HLTMuon::analyze(const edm::Handle<reco::MuonCollection>                 & Muon,
00166                       const edm::Handle<l1extra::L1MuonParticleCollection>    & MuCands1, 
00167                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands2,
00168                       const edm::Handle<edm::ValueMap<bool> >                 & isoMap2,
00169                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuCands3,
00170                       const edm::Handle<edm::ValueMap<bool> >                 & isoMap3,
00171                       const edm::Handle<reco::RecoChargedCandidateCollection> & oniaPixelCands,
00172                       const edm::Handle<reco::RecoChargedCandidateCollection> & oniaTrackCands,
00173                       const edm::Handle<reco::RecoChargedCandidateCollection> & MuNoVtxCands2, 
00174                       const reco::BeamSpot::Point & BSPosition,
00175                       TTree* HltTree) {
00176 
00177   //std::cout << " Beginning HLTMuon " << std::endl;
00178 
00179   if (Muon.isValid()) {
00180     reco::MuonCollection mymuons;
00181     mymuons = * Muon;
00182     std::sort(mymuons.begin(),mymuons.end(),PtGreater());
00183     nmuon = mymuons.size();
00184     typedef reco::MuonCollection::const_iterator muiter;
00185     int imu=0;
00186     for (muiter i=mymuons.begin(); i!=mymuons.end(); i++) 
00187     {
00188       muonpt[imu]         = i->pt();
00189       muonphi[imu]        = i->phi();
00190       muoneta[imu]        = i->eta();
00191       muonet[imu]         = i->et();
00192       muone[imu]          = i->energy(); 
00193       muontype[imu]       = i->type();
00194       muoncharge[imu]     = i->charge(); 
00195       muonTrkIsoR03[imu]  = i->isolationR03().sumPt;
00196       muonECalIsoR03[imu] = i->isolationR03().emEt;
00197       muonHCalIsoR03[imu] = i->isolationR03().hadEt;
00198       
00199       
00200       if (i->globalTrack().isNonnull())
00201       {
00202         muonchi2NDF[imu] = i->globalTrack()->normalizedChi2();
00203         muonD0[imu] = i->globalTrack()->dxy(BSPosition);
00204        }
00205       else 
00206       {
00207         muonchi2NDF[imu] = -99.;
00208         muonD0[imu] = -99.;}
00209       
00210       if (i->innerTrack().isNonnull()) muonNValidTrkHits[imu] = i->innerTrack()->numberOfValidHits();
00211       else muonNValidTrkHits[imu] = -99;
00212       
00213       if (i->isGlobalMuon()!=0) muonNValidMuonHits[imu] = i->globalTrack()->hitPattern().numberOfValidMuonHits();
00214       else muonNValidMuonHits[imu] = -99;
00215       
00216       imu++;
00217     }
00218   }
00219   else {nmuon = 0;}
00220 
00221   l1extra::L1MuonParticleCollection myMucands1; 
00222   myMucands1 = * MuCands1; 
00223   //  reco::RecoChargedCandidateCollection myMucands1;
00224   std::sort(myMucands1.begin(),myMucands1.end(),PtGreater()); 
00225 
00227 
00228   // Dealing with L2 muons
00229   reco::RecoChargedCandidateCollection myMucands2;
00230   if (MuCands2.isValid()) {
00231 //     reco::RecoChargedCandidateCollection myMucands2;
00232     myMucands2 = * MuCands2;
00233     std::sort(myMucands2.begin(),myMucands2.end(),PtGreater());
00234     nmu2cand = myMucands2.size();
00235     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00236     int imu2c=0;
00237     for (cand i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00238       reco::TrackRef tk = i->get<reco::TrackRef>();
00239 
00240       muonl2pt[imu2c] = tk->pt();
00241       // eta (we require |eta|<2.5 in all filters
00242       muonl2eta[imu2c] = tk->eta();
00243       muonl2phi[imu2c] = tk->phi();
00244 
00245       // Dr (transverse distance to (0,0,0))
00246       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00247       // However, we use |dr|<200 microns at L3, which it probably too tough for LHC startup
00248       muonl2dr[imu2c] = fabs(tk->dxy(BSPosition));
00249 
00250       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00251       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00252       muonl2dz[imu2c] = tk->dz(BSPosition);
00253 
00254       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00255       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00256       // Baseline cuts (HLT exercise):
00257       //                Relaxed Single muon:  ptLx>16 GeV
00258       //                Isolated Single muon: ptLx>11 GeV
00259       //                Relaxed Double muon: ptLx>3 GeV
00260        double l2_err0 = tk->error(0); // error on q/p
00261        double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00262 //       double ptLx = tk->pt();
00263       // convert 50% efficiency threshold to 90% efficiency threshold
00264       // For L2 muons: nsigma_Pt_ = 3.9
00265 //       double nsigma_Pt_ = 3.9;
00266       // For L3 muons: nsigma_Pt_ = 2.2
00267       // these are the old TDR values for nsigma_Pt_
00268       // We know that these values are slightly smaller for CMSSW
00269       // But as quoted above, we want to get rid of this gymnastics in the future
00270 //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00271 
00272       // Charge
00273       // We use the charge in some dimuon paths
00274                         muonl2pterr[imu2c] = l2_err0/l2_abspar0;
00275       muonl2chg[imu2c] = tk->charge();
00276       
00277       if (isoMap2.isValid()){
00278         // Isolation flag (this is a bool value: true => isolated)
00279         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap2)[tk];
00280         muonl2iso[imu2c] = muon1IsIsolated;
00281       }
00282       else {muonl2iso[imu2c] = -999;}
00283 
00284       //JH
00285       l1extra::L1MuonParticleRef l1; 
00286       int il2 = 0; 
00287       //find the corresponding L1 
00288       l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00289       il2++; 
00290       int imu1idx = 0; 
00291       if (MuCands1.isValid()) { 
00292         typedef l1extra::L1MuonParticleCollection::const_iterator candl1; 
00293         for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) { 
00294           if((j->pt() == l1->pt()) &&
00295              (j->eta() == l1->eta()) &&
00296              (j->phi() == l1->phi()) &&
00297              (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00298             {break;}
00299           //      std::cout << << std::endl;
00300           //          if ( tkl1 == l1 ) {break;} 
00301           imu1idx++; 
00302         } 
00303       } 
00304       else {imu1idx = -999;} 
00305       muonl21idx[imu2c] = imu1idx; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00306       //end JH
00307 
00308       imu2c++;
00309     }
00310   }
00311   else {nmu2cand = 0;}
00312 
00313   // Dealing with L3 muons
00314   reco::RecoChargedCandidateCollection myMucands3;
00315   if (MuCands3.isValid()) {
00316     myMucands3 = * MuCands3;
00317     std::sort(myMucands3.begin(),myMucands3.end(),PtGreater());
00318     nmu3cand = myMucands3.size();
00319     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00320     int imu3c=0;
00321     for (cand i=myMucands3.begin(); i!=myMucands3.end(); i++) {
00322       reco::TrackRef tk = i->get<reco::TrackRef>();
00323 
00324       reco::TrackRef staTrack;
00325       typedef reco::MuonTrackLinksCollection::const_iterator l3muon;
00326       int il3 = 0;
00327       //find the corresponding L2 track
00328       staTrack = tk->seedRef().castTo<edm::Ref< L3MuonTrajectorySeedCollection> >()->l2Track();
00329       il3++;
00330       int imu2idx = 0;
00331       if (MuCands2.isValid()) {
00332         typedef reco::RecoChargedCandidateCollection::const_iterator candl2;
00333         for (candl2 i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00334           reco::TrackRef tkl2 = i->get<reco::TrackRef>();
00335           if ( tkl2 == staTrack ) {break;}
00336           imu2idx++;
00337         }
00338       }
00339       else {imu2idx = -999;}
00340       muonl32idx[imu3c] = imu2idx; // Index of the L2 muon having matched with the L3 muon with index imu3c
00341       
00342       muonl3pt[imu3c] = tk->pt();
00343       // eta (we require |eta|<2.5 in all filters
00344       muonl3eta[imu3c] = tk->eta();
00345       muonl3phi[imu3c] = tk->phi();
00346 
00347 //       // Dr (transverse distance to (0,0,0))
00348 //       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00349 //       // However, we use |dr|<300 microns at L3, which it probably too tough for LHC startup
00350       muonl3dr[imu3c] = fabs(tk->dxy(BSPosition));
00351 
00352 //       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00353 //       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00354       muonl3dz[imu3c] = tk->dz(BSPosition);
00355 
00356 //       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00357 //       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00358 //       // Baseline cuts (HLT exercise):
00359 //       //                Relaxed Single muon:  ptLx>16 GeV
00360 //       //                Isolated Single muon: ptLx>11 GeV
00361 //       //                Relaxed Double muon: ptLx>3 GeV
00362         double l3_err0 = tk->error(0); // error on q/p
00363         double l3_abspar0 = fabs(tk->parameter(0)); // |q/p|
00364 // //       double ptLx = tk->pt();
00365 //       // convert 50% efficiency threshold to 90% efficiency threshold
00366 //       // For L2 muons: nsigma_Pt_ = 3.9
00367 //       // For L3 muons: nsigma_Pt_ = 2.2
00368 // //       double nsigma_Pt_ = 2.2;
00369 //       // these are the old TDR values for nsigma_Pt_
00370 //       // We know that these values are slightly smaller for CMSSW
00371 //       // But as quoted above, we want to get rid of this gymnastics in the future
00372 // //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00373 
00374       // Charge
00375       // We use the charge in some dimuon paths
00376       muonl3pterr[imu3c] = l3_err0/l3_abspar0;
00377       muonl3chg[imu3c] = tk->charge();
00378 
00379       if (isoMap3.isValid()){
00380         // Isolation flag (this is a bool value: true => isolated)
00381         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap3)[tk];
00382         muonl3iso[imu3c] = muon1IsIsolated;
00383       }
00384       else {muonl3iso[imu3c] = -999;}
00385 
00386       imu3c++;
00387     }
00388   }
00389   else {nmu3cand = 0;}
00390 
00391   // Dealing with L2 no-Vertex muons
00392   reco::RecoChargedCandidateCollection muNoVtxMucands2;
00393   if (MuNoVtxCands2.isValid()) {
00394     muNoVtxMucands2 = * MuNoVtxCands2;
00395     std::sort(muNoVtxMucands2.begin(),muNoVtxMucands2.end(),PtGreater());
00396     nmu2cand = muNoVtxMucands2.size();
00397     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00398     int imu2c=0;
00399     for (cand i=muNoVtxMucands2.begin(); i!=muNoVtxMucands2.end(); i++) {
00400       reco::TrackRef tk = i->get<reco::TrackRef>();
00401 
00402       muonl2novtxpt[imu2c] = tk->pt();
00403       muonl2novtxeta[imu2c] = tk->eta();
00404       muonl2novtxphi[imu2c] = tk->phi();
00405       muonl2novtxdr[imu2c] = fabs(tk->dxy(BSPosition));
00406       muonl2novtxdz[imu2c] = tk->dz(BSPosition);
00407 
00408       double l2_err0 = tk->error(0); // error on q/p
00409       double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00410       
00411       muonl2novtxpterr[imu2c] = l2_err0/l2_abspar0;
00412       muonl2novtxchg[imu2c] = tk->charge();
00413       
00414       l1extra::L1MuonParticleRef l1; 
00415       int il2 = 0; 
00416       //find the corresponding L1 
00417       l1 = tk->seedRef().castTo<edm::Ref< L2MuonTrajectorySeedCollection> >()->l1Particle();
00418       il2++; 
00419       int imu1idx = 0; 
00420       if (MuCands1.isValid()) { 
00421         typedef l1extra::L1MuonParticleCollection::const_iterator candl1; 
00422         for (candl1 j=myMucands1.begin(); j!=myMucands1.end(); j++) { 
00423           if((j->pt() == l1->pt()) &&
00424              (j->eta() == l1->eta()) &&
00425              (j->phi() == l1->phi()) &&
00426              (j->gmtMuonCand().quality() == l1->gmtMuonCand().quality()))
00427             {break;}
00428           imu1idx++; 
00429         } 
00430       } 
00431       else {imu1idx = -999;} 
00432       muonl2novtx1idx[imu2c] = imu1idx; // Index of the L1 muon having matched with the L2 muon with index imu2c 
00433 
00434       imu2c++;
00435     }
00436   }
00437   else {nmu2cand = 0;}
00438  
00439 
00440 
00441   // Dealing with Onia Pixel tracks
00442   reco::RecoChargedCandidateCollection myOniaPixelCands;
00443   if (oniaPixelCands.isValid()) {
00444     myOniaPixelCands = * oniaPixelCands;
00445     std::sort(myOniaPixelCands.begin(),myOniaPixelCands.end(),PtGreater());
00446     nOniaPixelCand = myOniaPixelCands.size();
00447     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00448     int ic=0;
00449     for (cand i=myOniaPixelCands.begin(); i!=myOniaPixelCands.end(); i++) {
00450       reco::TrackRef tk = i->get<reco::TrackRef>();
00451 
00452       oniaPixelpt[ic] = tk->pt();
00453       oniaPixeleta[ic] = tk->eta();
00454       oniaPixelphi[ic] = tk->phi();
00455       oniaPixeldr[ic] = tk->dxy(BSPosition);
00456       oniaPixeldz[ic] = tk->dz(BSPosition);
00457       oniaPixelchg[ic] = tk->charge();
00458       oniaPixelHits[ic] = tk->numberOfValidHits();
00459       oniaPixelNormChi2[ic] = tk->normalizedChi2();
00460 
00461       ic++;
00462     }
00463   }
00464   else {nOniaPixelCand = 0;}
00465 
00466   // Dealing with Onia Tracks
00467   reco::RecoChargedCandidateCollection myOniaTrackCands;
00468   if (oniaTrackCands.isValid()) {
00469     myOniaTrackCands = * oniaTrackCands;
00470     std::sort(myOniaTrackCands.begin(),myOniaTrackCands.end(),PtGreater());
00471     nOniaTrackCand = myOniaTrackCands.size();
00472     typedef reco::RecoChargedCandidateCollection::const_iterator cand;
00473     int ic=0;
00474     for (cand i=myOniaTrackCands.begin(); i!=myOniaTrackCands.end(); i++) {
00475       reco::TrackRef tk = i->get<reco::TrackRef>();
00476 
00477       oniaTrackpt[ic] = tk->pt();
00478       oniaTracketa[ic] = tk->eta();
00479       oniaTrackphi[ic] = tk->phi();
00480       oniaTrackdr[ic] = tk->dxy(BSPosition);
00481       oniaTrackdz[ic] = tk->dz(BSPosition);
00482       oniaTrackchg[ic] = tk->charge();
00483       oniaTrackHits[ic] = tk->numberOfValidHits();
00484       oniaTrackNormChi2[ic] = tk->normalizedChi2();
00485 
00486       ic++;
00487     }
00488   }
00489   else {nOniaTrackCand = 0;}
00490 
00491 
00493 
00494 
00495 
00496 }