CMSSW Reference Manual

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   vector<std::string> parameterNames = myEmParams.getParameterNames() ;
00027   
00028   for ( 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   const int kMaxMuonL2 = 500;
00041   muonl2pt = new float[kMaxMuonL2];
00042   muonl2phi = new float[kMaxMuonL2];
00043   muonl2eta = new float[kMaxMuonL2];
00044   muonl2dr = new float[kMaxMuonL2];
00045   muonl2dz = new float[kMaxMuonL2];
00046   muonl2chg = new int[kMaxMuonL2];
00047   muonl2pterr = new float[kMaxMuonL2];
00048   muonl2iso = new int[kMaxMuonL2];
00049   const int kMaxMuonL3 = 500;
00050   muonl3pt = new float[kMaxMuonL3];
00051   muonl3phi = new float[kMaxMuonL3];
00052   muonl3eta = new float[kMaxMuonL3];
00053   muonl3dr = new float[kMaxMuonL3];
00054   muonl3dz = new float[kMaxMuonL3];
00055   muonl3chg = new int[kMaxMuonL3];
00056   muonl3pterr = new float[kMaxMuonL3];
00057   muonl3iso = new int[kMaxMuonL3];
00058   muonl32idx = new int[kMaxMuonL3];
00059 
00060   // Muon-specific branches of the tree 
00061   HltTree->Branch("NrecoMuon",&nmuon,"NrecoMuon/I");
00062   HltTree->Branch("recoMuonPt",muonpt,"recoMuonPt[NrecoMuon]/F");
00063   HltTree->Branch("recoMuonPhi",muonphi,"recoMuonPhi[NrecoMuon]/F");
00064   HltTree->Branch("recoMuonEta",muoneta,"recoMuonEta[NrecoMuon]/F");
00065   HltTree->Branch("recoMuonEt",muonet,"recoMuonEt[NrecoMuon]/F");
00066   HltTree->Branch("recoMuonE",muone,"recoMuonE[NrecoMuon]/F");
00067   HltTree->Branch("NohMuL2",&nmu2cand,"NohMuL2/I");
00068   HltTree->Branch("ohMuL2Pt",muonl2pt,"ohMuL2Pt[NohMuL2]/F");
00069   HltTree->Branch("ohMuL2Phi",muonl2phi,"ohMuL2Phi[NohMuL2]/F");
00070   HltTree->Branch("ohMuL2Eta",muonl2eta,"ohMuL2Eta[NohMuL2]/F");
00071   HltTree->Branch("ohMuL2Chg",muonl2chg,"ohMuL2Chg[NohMuL2]/I");
00072   HltTree->Branch("ohMuL2PtErr",muonl2pterr,"ohMuL2PtErr[NohMuL2]/F");
00073   HltTree->Branch("ohMuL2Iso",muonl2iso,"ohMuL2Iso[NohMuL2]/I");
00074   HltTree->Branch("ohMuL2Dr",muonl2dr,"ohMuL2Dr[NohMuL2]/F");
00075   HltTree->Branch("ohMuL2Dz",muonl2dz,"ohMuL2Dz[NohMuL2]/F");
00076   HltTree->Branch("NohMuL3",&nmu3cand,"NohMuL3/I");
00077   HltTree->Branch("ohMuL3Pt",muonl3pt,"ohMuL3Pt[NohMuL3]/F");
00078   HltTree->Branch("ohMuL3Phi",muonl3phi,"ohMuL3Phi[NohMuL3]/F");
00079   HltTree->Branch("ohMuL3Eta",muonl3eta,"ohMuL3Eta[NohMuL3]/F");
00080   HltTree->Branch("ohMuL3Chg",muonl3chg,"ohMuL3Chg[NohMuL3]/I");
00081   HltTree->Branch("ohMuL3PtErr",muonl3pterr,"ohMuL3PtErr[NohMuL3]/F");
00082   HltTree->Branch("ohMuL3Iso",muonl3iso,"ohMuL3Iso[NohMuL3]/I");
00083   HltTree->Branch("ohMuL3Dr",muonl3dr,"ohMuL3Dr[NohMuL3]/F");
00084   HltTree->Branch("ohMuL3Dz",muonl3dz,"ohMuL3Dz[NohMuL3]/F");
00085   HltTree->Branch("ohMuL3L2idx",muonl32idx,"ohMuL3L2idx[NohMuL3]/I");
00086 
00087 }
00088 
00089 /* **Analyze the event** */
00090 void HLTMuon::analyze(const edm::Handle<MuonCollection>                 & Muon,
00091                       const edm::Handle<RecoChargedCandidateCollection> & MuCands2,
00092                       const edm::Handle<edm::ValueMap<bool> >           & isoMap2,
00093                       const edm::Handle<RecoChargedCandidateCollection> & MuCands3,
00094                       const edm::Handle<edm::ValueMap<bool> >           & isoMap3,
00095                       const edm::Handle<MuonTrackLinksCollection>       & mulinks,
00096                       TTree* HltTree) {
00097 
00098   //std::cout << " Beginning HLTMuon " << std::endl;
00099 
00100   if (Muon.isValid()) {
00101     MuonCollection mymuons;
00102     mymuons = * Muon;
00103     std::sort(mymuons.begin(),mymuons.end(),PtGreater());
00104     nmuon = mymuons.size();
00105     typedef MuonCollection::const_iterator muiter;
00106     int imu=0;
00107     for (muiter i=mymuons.begin(); i!=mymuons.end(); i++) {
00108       muonpt[imu] = i->pt();
00109       muonphi[imu] = i->phi();
00110       muoneta[imu] = i->eta();
00111       muonet[imu] = i->et();
00112       muone[imu] = i->energy();
00113       imu++;
00114     }
00115   }
00116   else {nmuon = 0;}
00117 
00119 
00120   // Dealing with L2 muons
00121   RecoChargedCandidateCollection myMucands2;
00122   if (MuCands2.isValid()) {
00123 //     RecoChargedCandidateCollection myMucands2;
00124     myMucands2 = * MuCands2;
00125     std::sort(myMucands2.begin(),myMucands2.end(),PtGreater());
00126     nmu2cand = myMucands2.size();
00127     typedef RecoChargedCandidateCollection::const_iterator cand;
00128     int imu2c=0;
00129     for (cand i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00130       TrackRef tk = i->get<TrackRef>();
00131 
00132       muonl2pt[imu2c] = tk->pt();
00133       // eta (we require |eta|<2.5 in all filters
00134       muonl2eta[imu2c] = tk->eta();
00135       muonl2phi[imu2c] = tk->phi();
00136 
00137       // Dr (transverse distance to (0,0,0))
00138       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00139       // However, we use |dr|<200 microns at L3, which it probably too tough for LHC startup
00140       muonl2dr[imu2c] = fabs(tk->d0());
00141 
00142       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00143       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00144       muonl2dz[imu2c] = fabs(tk->dz());
00145 
00146       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00147       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00148       // Baseline cuts (HLT exercise):
00149       //                Relaxed Single muon:  ptLx>16 GeV
00150       //                Isolated Single muon: ptLx>11 GeV
00151       //                Relaxed Double muon: ptLx>3 GeV
00152        double l2_err0 = tk->error(0); // error on q/p
00153        double l2_abspar0 = fabs(tk->parameter(0)); // |q/p|
00154 //       double ptLx = tk->pt();
00155       // convert 50% efficiency threshold to 90% efficiency threshold
00156       // For L2 muons: nsigma_Pt_ = 3.9
00157 //       double nsigma_Pt_ = 3.9;
00158       // For L3 muons: nsigma_Pt_ = 2.2
00159       // these are the old TDR values for nsigma_Pt_
00160       // We know that these values are slightly smaller for CMSSW
00161       // But as quoted above, we want to get rid of this gymnastics in the future
00162 //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00163 
00164       // Charge
00165       // We use the charge in some dimuon paths
00166                         muonl2pterr[imu2c] = l2_err0/l2_abspar0;
00167       muonl2chg[imu2c] = tk->charge();
00168 
00169       if (isoMap2.isValid()){
00170         // Isolation flag (this is a bool value: true => isolated)
00171         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap2)[tk];
00172         muonl2iso[imu2c] = muon1IsIsolated;
00173       }
00174       else {muonl2iso[imu2c] = -999;}
00175 
00176       imu2c++;
00177     }
00178   }
00179   else {nmu2cand = 0;}
00180 
00181   // Dealing with L3 muons
00182   RecoChargedCandidateCollection myMucands3;
00183   if (MuCands3.isValid()) {
00184     myMucands3 = * MuCands3;
00185     std::sort(myMucands3.begin(),myMucands3.end(),PtGreater());
00186     nmu3cand = myMucands3.size();
00187     typedef RecoChargedCandidateCollection::const_iterator cand;
00188     int imu3c=0;
00189     for (cand i=myMucands3.begin(); i!=myMucands3.end(); i++) {
00190       TrackRef tk = i->get<TrackRef>();
00191 
00192       TrackRef staTrack;
00193       typedef MuonTrackLinksCollection::const_iterator l3muon;
00194       int il3 = 0;
00195       for (l3muon j=mulinks->begin(); j!=mulinks->end(); j++){
00196         if (j->globalTrack() == tk) {
00197           staTrack = j->standAloneTrack();
00198           break;
00199         }
00200         il3++;
00201       }
00202       int imu2idx = 0;
00203       if (MuCands2.isValid()) {
00204         typedef RecoChargedCandidateCollection::const_iterator candl2;
00205         for (candl2 i=myMucands2.begin(); i!=myMucands2.end(); i++) {
00206           TrackRef tkl2 = i->get<TrackRef>();
00207           if ( tkl2 == staTrack ) {break;}
00208           imu2idx++;
00209         }
00210       }
00211       else {imu2idx = -999;}
00212       muonl32idx[imu3c] = imu2idx; // Index of the L2 muon having matched with the L3 muon with index imu3c
00213       
00214       muonl3pt[imu3c] = tk->pt();
00215       // eta (we require |eta|<2.5 in all filters
00216       muonl3eta[imu3c] = tk->eta();
00217       muonl3phi[imu3c] = tk->phi();
00218 
00219 //       // Dr (transverse distance to (0,0,0))
00220 //       // For baseline triggers, we do no cut at L2 (|dr|<9999 cm)
00221 //       // However, we use |dr|<300 microns at L3, which it probably too tough for LHC startup
00222       muonl3dr[imu3c] = fabs(tk->d0());
00223 
00224 //       // Dz (longitudinal distance to z=0 when at minimum transverse distance)
00225 //       // For baseline triggers, we do no cut (|dz|<9999 cm), neither at L2 nor at L3
00226       muonl3dz[imu3c] = fabs(tk->dz());
00227 
00228 //       // At present we do not cut on this, but on a 90% CL value "ptLx" defined here below
00229 //       // We should change this in the future and cut directly on "pt", to avoid unnecessary complications and risks
00230 //       // Baseline cuts (HLT exercise):
00231 //       //                Relaxed Single muon:  ptLx>16 GeV
00232 //       //                Isolated Single muon: ptLx>11 GeV
00233 //       //                Relaxed Double muon: ptLx>3 GeV
00234         double l3_err0 = tk->error(0); // error on q/p
00235         double l3_abspar0 = fabs(tk->parameter(0)); // |q/p|
00236 // //       double ptLx = tk->pt();
00237 //       // convert 50% efficiency threshold to 90% efficiency threshold
00238 //       // For L2 muons: nsigma_Pt_ = 3.9
00239 //       // For L3 muons: nsigma_Pt_ = 2.2
00240 // //       double nsigma_Pt_ = 2.2;
00241 //       // these are the old TDR values for nsigma_Pt_
00242 //       // We know that these values are slightly smaller for CMSSW
00243 //       // But as quoted above, we want to get rid of this gymnastics in the future
00244 // //       if (abspar0>0) ptLx += nsigma_Pt_*err0/abspar0*tk->pt();
00245 
00246       // Charge
00247       // We use the charge in some dimuon paths
00248       muonl3pterr[imu3c] = l3_err0/l3_abspar0;
00249       muonl3chg[imu3c] = tk->charge();
00250 
00251       if (isoMap3.isValid()){
00252         // Isolation flag (this is a bool value: true => isolated)
00253         edm::ValueMap<bool> ::value_type muon1IsIsolated = (*isoMap3)[tk];
00254         muonl3iso[imu3c] = muon1IsIsolated;
00255       }
00256       else {muonl3iso[imu3c] = -999;}
00257 
00258       imu3c++;
00259     }
00260   }
00261   else {nmu3cand = 0;}
00262 
00264 
00265 
00266 
00267 }
HLTMuon.cc
