00001 #include "caf_util/MuonSelector.hpp"
00002 #include "cafe/Config.hpp"
00003
00004 #include <string>
00005 #include <sstream>
00006 #include <fstream>
00007 #include <stdexcept>
00008
00009 using namespace std ;
00010 using namespace cafe ;
00011
00012 namespace {
00013 bool moreThan(const TMBMuon& a, const TMBMuon& b) {
00014 return a.Pt() > b.Pt();
00015 }
00016 }
00017
00018 namespace caf_util {
00019
00020 MuonSelector::MuonSelector(const char *name) : SelectUserObjects<TMBMuon>(name),
00021 _nselected(0), _event(0),
00022 _jetvars("fX", "fY", "fZ"),_runvar("_runno"),
00023 _chunkvar("_nchunks"), _mc(false)
00024 {
00025
00026 cafe::Config config(name);
00027
00028
00029 _nmuons = config.get("nMuons", 1);
00030 _nmuonsmax = config.get("nMuonsMax", -1);
00031
00032 _pTcut = config.get("pT", 0.0);
00033 _eta = config.get("eta", 2.0) ;
00034
00035 _MinDeltaR = config.get("MinDeltaR", -1.0);
00036 _MaxDeltaR = config.get("MaxDeltaR", -1.0);
00037 _etHalo = config.get("etHalo", -1.0) ;
00038 _etHaloMin = config.get("etHaloMin", -1.0) ;
00039 _etHaloScaled = config.get("etHaloScaled", -1.0) ;
00040 _etHaloScaledMin = config.get("etHaloScaledMin", -1.0) ;
00041 _etTrkCone = config.get("etTrkCone", -1.0) ;
00042 _etTrkConeMin = config.get("etTrkConeMin", -1.0) ;
00043 _etTrkConeScaled = config.get("etTrkConeScaled", -1.0) ;
00044 _etTrkConeScaledMin = config.get("etTrkConeScaledMin", -1.0) ;
00045 _excludeHole = config.get("excludeHole",false);
00046
00047 _use_muid = config.get("UseMUID",0);
00048
00049 if(!_use_muid) {
00050
00051 err() << "WARNING: Standard muon ID not used - MUID group efficiencies not available for these cuts" << endl;
00052 _use_beam_spot = config.get("UseBeamSpot",false);
00053
00054 _nseg = config.get("Nseg", 3) ;
00055 _useCentralMatched = config.get("UseCentralMatched", true) ;
00056 _chisq = config.get("chisq", 4.0) ;
00057 _cosmicVeto = config.get("cosmicVeto", false) ;
00058 _dcaSignif = config.get("dcaSignif", -1.0) ;
00059 _dcaSMT = config.get("dcaSMT",-1.0);
00060 _dcaNoSMT = config.get("dcaNoSMT",-1.0);
00061 _nSMT = config.get("nSMT",-1);
00062
00063 string q = config.get("Quality", "medium") ;
00064 if ( q == "LOOSE" || q == "Loose" || q == "loose")
00065 _qual = Loose ;
00066 else if ( q == "TIGHT" || q == "Tight" || q == "tight")
00067 _qual = Tight ;
00068 else if( q == "MEDIUM" || q == "Medium" || q == "medium")
00069 _qual = Medium ;
00070 else {
00071 err() << "MUON SELECTOR WARNING! "
00072 << "Parameter Quality = \"" << q
00073 << "\" is not recognized. "
00074 << "Must be loose, medium or tight. "
00075 << "Will use default quality MEDIUM."
00076 << endl ;
00077 _qual = Medium ;
00078 }
00079 }
00080 else {
00081 _path = config.get("Path", "caf_util/configs/");
00082
00083 _muqual_file = config.get("MuonQualityFile","mucuts.txt");
00084 _muqual_name = config.get("MuonQualityName","MediumNseg3");
00085 _muqual_version = config.get("MuonQualityVersion",1);
00086
00087 _trkqual_file = config.get("TrkQualityFile","mutrkcuts.txt");
00088 _trkqual_name = config.get("TrkQualityName","None");
00089 _trkqual_version = config.get("TrkQualityVersion",1);
00090
00091 _isoqual_file = config.get("IsoQualityFile","muisocuts.txt");
00092 _isoqual_name = config.get("IsoQualityName","None");
00093 _isoqual_version = config.get("IsoQualityVersion",1);
00094
00095 _cuts_set = SetStdCuts(_path, _muqual_file,_muqual_name,_muqual_version,_trkqual_file,_trkqual_name,_trkqual_version, _isoqual_file, _isoqual_name, _isoqual_version);
00096 if(_cuts_set) {
00097 DisplayMuonQualityCuts();
00098 DisplayMuonTrackQualityCuts();
00099 DisplayMuonIsolationQualityCuts();
00100 out() << endl;
00101 out() << "Muon Cuts were set up successfully" << endl;
00102 out() << "
00103 out() << endl;
00104 }
00105 else throw runtime_error("Error setting the std muon id cuts");
00106 } //std id
00107
00108 out() << endl << "Muon cuts:" << endl;
00109 // Cut controlled by std ID only need printing if
00110 // std ID is turned off
00111 if(!_use_muid) {
00112 out() << "Muon quality = ";
00113 if(_qual==Loose) out() << "Loose" << endl;
00114 if(_qual==Medium) out() << "Medium" << endl;
00115 if(_qual==Tight) out() << "Tight" << endl;
00116 out() << "N Seg >= " << _nseg << endl;
00117 if(_useCentralMatched) out() << "Getting central matched muons" << endl;
00118 if(_chisq > 0) out() << "Muon Track Chi2 / NDof < " << _chisq << endl;
00119 if(_cosmicVeto) out() << "Vetoing cosmics" << endl;
00120 if(_dcaSignif > 0) out() << "DCA Significance < " << _dcaSignif << endl;
00121 if(_dcaSMT > 0) out() << "DCA for muon tracks with >0 SMT hits < " << _dcaSMT << endl;
00122 if(_dcaNoSMT > 0) out() << "DCA for muon tracks with 0 SMT hits < " << _dcaNoSMT << endl;;
00123 if(_nSMT > 0) out() << "nSMT hits >= " << _nSMT << endl;
00124 if(_etHalo > 0) out() << "Muon etHalo < " << _etHalo << endl;
00125 if(_etHaloMin > 0) out() << "Muon etHalo >= " << _etHaloMin << endl;
00126 if(_etHaloScaled > 0) out() << "Muon etHaloScaled < " << _etHaloScaled << endl;
00127 if(_etHaloScaledMin > 0) out() << "Muon etHaloScaled >= " << _etHaloScaledMin << endl;
00128 if(_etTrkCone > 0) out() << "Muon etTrkCone < " << _etTrkCone << endl;
00129 if(_etTrkConeMin > 0) out() << "Muon etTrkCone >= " << _etTrkConeMin << endl;
00130 if(_etTrkConeScaled > 0) out() << "Muon etTrkConeScaled < " << _etTrkConeScaled << endl;
00131 if(_etTrkConeScaledMin > 0) out() << "Muon etTrkConeScaled >= " << _etTrkConeScaledMin << endl;
00132 if(_MinDeltaR > 0 || _MaxDeltaR >0) {
00133 out() << "Removing muons near jets from branch " << _JetBranch.c_str();
00134 if(_MinDeltaR > 0) out() << ", dR > " << _MinDeltaR;
00135 if(_MaxDeltaR > 0) out() << ", dR < " << _MaxDeltaR;
00136 out() << endl;
00137 }
00138 }
00139 if(_pTcut > 0) out() << "Muon pT > " << _pTcut << endl;
00140 if(_eta > 0) out() << "Muon deteta < " << _eta << endl;
00141 if(_reject_no_beam_spot) out() << "Recjecting runs not in beam spot file" << endl;
00142 if(_corr_pT_SMT) out() << "Correcting pT of Muons with >0 SMT hits for beam spot posn." << endl;
00143 if(_corr_pT_noSMT) out() << "Correcting pT of Muons with 0 SMT hits for beam spot posn." << endl;
00144
00145
00146 out() << "Remove muons in the bottom hole: " ;
00147 if(_excludeHole) out() << "true" << endl ;
00148 else out() << "false" << endl ;
00149
00150 //If needed setup beam spot file
00151 if(_use_beam_spot) {
00152 _beam = config.get("BeamSpotFile","caf_util/configs/beamSpot-2.09");
00153 _reject_no_beam_spot = config.get("RejectNoBeamSpot",0);
00154 char dataname[50];
00155 sprintf(dataname,"%s",_beam.c_str());
00156 if((_beamspot_file = fopen(dataname,"r")) == NULL) {
00157 cerr << "Can't open beam spot file: " << dataname << endl;
00158 cerr << "Beam spot posn set to 0." << endl;
00159 cerr << "Pt corrections will not be applied." << endl;
00160 _bx=_by=_bsx=_bsy=0;
00161 _open_beam = false;
00162 }
00163 else {
00164 out() << "Opened beamspot file: " << dataname << endl;
00165 _open_beam = true;
00166 _bx=_by=_bsx=_bsy=0; //set to 0 in case using mc
00167 }
00168 }
00169 else _open_beam = false;
00170
00171 // If using beam spot have option to correct track pT
00172 if(_open_beam && _use_beam_spot) {
00173 _corr_pT_SMT = config.get("CorrectPtSMT",0);
00174 _corr_pT_noSMT = config.get("CorrectPtNoSMT",1);
00175 }
00176 else {
00177 _corr_pT_SMT = false;
00178 _corr_pT_noSMT = false;
00179 }
00180
00181 _sort = config.get("Sort", 0);
00182 _JetBranch = config.get("JetBranch", "");
00183
00184 if( (_MinDeltaR > 0 || _MaxDeltaR > 0) &&
00185 _JetBranch == "" )
00186 throw runtime_error("caf_util::MuonSelector: DeltaR cut requested, but no JetBranch specified");
00187
00188 _finder.SetVariablesTrack(Variables("fUniqueID", "_chi2_ndof", "_rdca", "_trerrs")) ;
00189 _finder.addTrackVariable("_z");
00190 _finder.addTrackVariable("_q");
00191 _finder.addTrackVariable("fX");
00192 _finder.addTrackVariable("fY");
00193 _finder.addTrackVariable("fZ");
00194 _finder.addTrackVariable("_hitmask");
00195
00196 _finder.SetVariablesVertex(Variables("fUniqueID","_VertexX","_VertexY","_VertexZ"));
00197 _finder.addVertexVariable("_CovXX");
00198 _finder.addVertexVariable("_CovYY");
00199 _finder.addVertexVariable("_CovZZ");
00200
00201 _debuglevel = config.get("DebugLevel",0);
00202 };
00203
00204 bool MuonSelector::processEvent(cafe::Event &event)
00205 {
00206
00207 //get pointer to statistics collector
00208 event.get("StatPointer", _stat) ;
00209 _event = &event ;
00210 //
00211 // get the jet branch if the DeltaR cut is wanted
00212 //
00213 if (_JetBranch != "" && (_MinDeltaR > 0 || _MaxDeltaR > 0 )) {
00214 _jets = event.getJets(_JetBranch.c_str(),_jetvars);
00215 }
00216 //
00217 // Get Beam Spot if needed
00218 //
00219 if(_use_beam_spot && _open_beam && !_mc) {
00220 int n_mc_chunks = (event.getMCEventInfo(_chunkvar))->nchunks();
00221 // If event is MC don't use beam spot file
00222 if(n_mc_chunks>0) {
00223 out() << "Sample is MC - beam spot file ignored" << endl;
00224 _mc=true;
00225 }
00226 else _mc=false;
00227 if(_use_beam_spot && !_mc) {
00228 int run = (event.getGlobal(_runvar))->runno();
00229 if(run!=_lastrun) { //if same run don't need to load new beamspot
00230 bool use=false;
00231 use = GetBeamSpot(_beamspot_file, run, _bx, _by, _bsx, _bsy);
00232 if(use) _lastrun=run;
00233 //reject event if no beam spot found for that run
00234 if(!use && _reject_no_beam_spot) return false;
00235 }
00236 }
00237 //put beam spot info into event for later processors
00238 event.put("Bx",_bx);
00239 event.put("By",_by);
00240 event.put("Bsx",_bsx);
00241 event.put("Bsy",_bsy);
00242 }
00243 //
00244 // apply selections to the muons
00245 //
00246 // initialize
00247 //
00248 _nselected = 0 ;
00249 SelectUserObjects<TMBMuon>::processEvent(event) ;
00250
00251 for (int i = 1; i <= _nmuons; i++) {
00252 if (_nselected < i) return false ;
00253 ostringstream st ;
00254 st << "Number of muons >= " << i ;
00255 _stat.EventSelected(st.str()) ;
00256 }
00257
00258 if (_nmuonsmax >= 0) {
00259 if (_nselected > _nmuonsmax) return false ;
00260 ostringstream st ;
00261 st << "Number of muons <= " << _nmuonsmax ;
00262 _stat.EventSelected(st.str()) ;
00263 }
00264
00265 return true ;
00266 };
00267
00268 bool MuonSelector::selectObject(const TMBMuon &muon)
00269 {
00270
00271 if (_excludeHole) {
00272 if (muon.isInHole()) return false;
00273 {
00274 ostringstream st ;
00275 if (_nmuons>0) _stat.EventSelected("Muon out of bottom hole") ;
00276 }
00277
00278 }
00279
00280 if(_eta > 0) {
00281 if (fabs(muon.detectorEta()) >= _eta) return false ;
00282 {
00283 ostringstream st ;
00284 st << "Muon detector eta < " << (float) _eta ;
00285 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00286 }
00287 }
00288
00289 string qual ;
00290 switch (_qual) {
00291 case Loose : {
00292 if (!muon.isLoose()) return false ;
00293 qual = "loose" ;
00294 } ; break ;
00295 case Medium : {
00296 if (!muon.isMedium()) return false ;
00297 qual = "medium" ;
00298 } ; break ;
00299 case Tight : {
00300 if (!muon.isTight()) return false ;
00301 qual = "tight" ;
00302 } ; break ;
00303 default : return false ; break ;
00304 }
00305 if (_nmuons>0) _stat.EventSelected("Muon quality is " + qual) ;
00306
00307 if (abs(muon.nseg()) < abs(_nseg)) return false ;
00308 {
00309 ostringstream st ;
00310 st << "Number of layers >= " << _nseg ;
00311 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00312 }
00313
00314 if (_cosmicVeto) {
00315 if (muon.isCosmic()) return false ;
00316 if (_nmuons>0)_stat.EventSelected("Veto on cosmic muon") ;
00317 }
00318
00319 if (_useCentralMatched) {
00320 if (!muon.hasCentral()) return false ;
00321 if (_nmuons>0)_stat.EventSelected("Matched with central track") ;
00322 }
00323
00324 const TMBTrack* trk = 0 ;
00325 if (muon.hasCentral() && (_chisq > 0. || _dcaSignif > 0. || _dcaSMT >0. || _dcaNoSMT > 0. || _nSMT >= 0 || _corr_pT_SMT || _corr_pT_noSMT)) {
00326 trk =_finder.FindTrack(*_event, muon.GetChargedTrackRef()) ;
00327
00328 if (!trk) {
00329 err() << "MUON SELECTOR ERROR!"
00330 << " Track pointer from the track matched muon is 0!"
00331 << endl ;
00332 return false ;
00333 }
00334 }
00335
00336 if (muon.hasCentral() && _nSMT > 0) {
00337 if(trk->nsmt() < _nSMT) return false;
00338 ostringstream st ;
00339 st << "N SMT hits on Muon Track >= " << _nSMT;
00340 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00341 }
00342
00343 if (muon.hasCentral() && _chisq > 0 ) {
00344 if (trk->getChi2Ndf() >= _chisq) return false ;
00345 ostringstream st ;
00346 st << "Muon global fit chi square < " << (float) _chisq ;
00347 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00348 }
00349
00350 double dca = 0;
00351 if (muon.hasCentral() && (_dcaSignif > 0. || _dcaSMT > 0. || _dcaNoSMT > 0. || _corr_pT_SMT || _corr_pT_noSMT) ) {
00352
00353 if(!_use_beam_spot) {
00354 const TMBVertex* vtx = _finder.FindVertex(*_event,
00355 muon.GetVertexRef());
00356
00357 if (!vtx) {
00358 err() << "MUON SELECTOR ERROR!"
00359 << " Vertex pointer from the track matched muon is 0!"
00360 << endl ;
00361 if(debug()>0) {
00362 Collection<TMBPrimaryVertex> all_vtxs = (*_event).getPrimaryVertices(_jetvars);
00363 int n_vtx = all_vtxs.size();
00364 err() << "Event has " << n_vtx << " vertices." << endl;
00365 }
00366 return false ;
00367 }
00368
00369 double ip[2];
00370 double iperr[3];
00371 trk->impact(vtx,ip,iperr);
00372
00373 if(_dcaSignif >= 0.) {
00374 if (fabs(ip[0])/sqrt(iperr[0]) >= _dcaSignif)
00375 return false ;
00376 ostringstream st ;
00377 st << "DCA significance < " << (float) _dcaSignif ;
00378 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00379 }
00380 dca = ip[0];
00381 if(_dcaSMT >= 0. || _dcaNoSMT >= 0.) {
00382 if(trk->nsmt() > 0 && fabs(dca) >= _dcaSMT && _dcaSMT > 0.) return false;
00383 if (trk->nsmt() == 0 && fabs(dca) >= _dcaNoSMT && _dcaNoSMT > 0.) return false;
00384 ostringstream st ;
00385 if(_dcaNoSMT > 0. && _dcaSMT > 0. )
00386 st << "DCA < " << (float) _dcaSMT << " for nSMT > 0 and DCA < " << (float) _dcaNoSMT << " for nSMT = 0";
00387 else {
00388 if(_dcaSMT > 0.) st << "DCA < " << (float) _dcaSMT << " for nSMT > 0";
00389 if(_dcaNoSMT > 0.) st << "DCA < " << (float) _dcaNoSMT << " for nSMT = 0";
00390 }
00391 if (_nmuons>0) _stat.EventSelected(st.str());
00392 }
00393 } // !_use_beam_spot
00394 else {
00395 dca = DCAfromBeamSpot(trk);
00396 if(_dcaSignif >0. ) throw runtime_error("DCA Significance cut not cuttently supported for beam spot DCA");
00397 if(_dcaSMT > 0. || _dcaNoSMT > 0.) {
00398 if(trk->nsmt() > 0 && fabs(dca) >= _dcaSMT && _dcaSMT > 0.) return false;
00399 if (trk->nsmt() == 0 && fabs(dca) >= _dcaNoSMT && _dcaNoSMT > 0.) return false;
00400 ostringstream st ;
00401 if(_dcaNoSMT > 0. && _dcaSMT > 0. )
00402 st << "DCA < " << (float) _dcaSMT << " for nSMT > 0 and DCA < " << (float) _dcaNoSMT << " for nSMT = 0";
00403 else {
00404 if(_dcaSMT > 0.) st << "DCA < " << (float) _dcaSMT << " for nSMT > 0";
00405 if(_dcaNoSMT > 0.) st << "DCA < " << (float) _dcaNoSMT << " for nSMT = 0";
00406 }
00407 if (_nmuons>0) _stat.EventSelected(st.str());
00408 }
00409 }
00410
00411 } //If dca cut
00412
00413 //Since pT correction done in before
00414 //Can now just cut on muon.Pt()
00415 if (muon.Pt() < _pTcut) return false ;
00416
00417 {
00418 ostringstream st ;
00419 if(_corr_pT_noSMT || _corr_pT_SMT) st << "Corrected ";
00420 st << "Muon pT >= " << _pTcut << " GeV" ;
00421 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00422 }
00423
00424 if (_etHalo > 0) {
00425 if (muon.etHalo() >= _etHalo) return false ;
00426 ostringstream st ;
00427 st << "Calorimeter isolation < " << _etHalo << " GeV" ;
00428 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00429 }
00430 if (_etHaloMin > 0) {
00431 if (muon.etHalo() < _etHaloMin) return false ;
00432 ostringstream st ;
00433 st << "Calorimeter isolation > " << _etHaloMin << " GeV" ;
00434 if (_nmuons>0) _stat.EventSelected(st.str()) ;
00435 }
00436
00437 if (_etHaloScaled > 0) {
00438 if (muon.etHalo() / muon.Pt() >= _etHaloScaled) return false ;
00439 ostringstream st ;
00440 st << "Calorimeter scaled isolation < " << (float) _etHaloScaled ;
00441 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00442 }
00443
00444 if (_etHaloScaledMin > 0) {
00445 if (muon.etHalo() / muon.Pt() < _etHaloScaledMin) return false ;
00446 ostringstream st ;
00447 st << "Calorimeter scaled isolation >= " << (float) _etHaloScaledMin ;
00448 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00449 }
00450
00451 if (_etTrkCone > 0) {
00452 if (muon.etTrkCone5() >= _etTrkCone) return false ;
00453 ostringstream st ;
00454 st << "Track isolation < " << _etTrkCone << " GeV" ;
00455 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00456 }
00457 if (_etTrkConeMin > 0) {
00458 if (muon.etTrkCone5() < _etTrkConeMin) return false ;
00459 ostringstream st ;
00460 st << "Track isolation > " << _etTrkConeMin << " GeV" ;
00461 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00462 }
00463
00464 if (_etTrkConeScaled > 0) {
00465 if (muon.etTrkCone5() / muon.Pt() >= _etTrkConeScaled) return false ;
00466 ostringstream st ;
00467 st << "Track scaled isolation < " << (float) _etTrkConeScaled ;
00468 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00469 }
00470
00471 if (_etTrkConeScaledMin > 0) {
00472 if (muon.etTrkCone5() / muon.Pt() < _etTrkConeScaledMin) return false ;
00473 ostringstream st ;
00474 st << "Track scaled isolation >= " << (float) _etTrkConeScaledMin ;
00475 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00476 }
00477
00478 //
00479 // apply the DeltaR cut
00480 //
00481
00482 if (_jets.size() > 0){
00483 double deltaR;
00484 double deltaR_nearestJet = 100.0;
00485 for (unsigned int i=0; i < _jets.size(); i++) {
00486 deltaR = kinem::delta_R(_jets[i].Eta(), _jets[i].Phi(), muon.Eta(), muon.Phi());
00487 if (deltaR < deltaR_nearestJet)
00488 deltaR_nearestJet = deltaR;
00489 } // Loop over jets
00490 // Check if user wants to apply cut on DeltaR
00491 if (_MinDeltaR > 0){
00492 if (deltaR_nearestJet <= _MinDeltaR)
00493 return false;
00494 ostringstream st ;
00495 st << "DeltaR(mu, jet) > " << (float) _MinDeltaR ;
00496 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00497 } // if (_MinDeltaR > 0)
00498 if (_MaxDeltaR > 0){
00499 if (deltaR_nearestJet > _MaxDeltaR)
00500 return false;
00501 ostringstream st ;
00502 st << "DeltaR(mu, jet) <= " << (float) _MaxDeltaR ;
00503 if (_nmuons>0)_stat.EventSelected(st.str()) ;
00504 } // if (_MaxDeltaR > 0)
00505 } // if (_jets.size() > 0)
00506
00507 _nselected++ ;
00508
00509 if (debug()>1)
00510 out() << "PROCESSOR \"" << name()
00511 << "\" SELECTED OBJECT: pT = " << muon.Pt()
00512 << " eta = " << muon.Eta()
00513 << " phi = " << muon.Phi()
00514 << endl ;
00515 return true ;
00516
00517 };
00518
00519 void MuonSelector::before(cafe::Collection<TMBMuon>& from) {
00520
00521
00522 if(_corr_pT_SMT || _corr_pT_noSMT) {
00523 for(Collection<TMBMuon>::iterator it_muon = from.begin(); it_muon != from.end(); ++it_muon) {
00524
00525
00526 if(it_muon->hasCentral()) {
00527 const TMBTrack* trk = it_muon->GetChargedTrack();
00528 if(trk) {
00529 int smt = trk->nsmt();
00530 if(_corr_pT_SMT && smt>0) {
00531 double dca = DCAfromBeamSpot(trk);
00532 if(_debuglevel>1) out() << "Muon pT before correction = " << it_muon->Pt() << endl;
00533 it_muon->CorrectPt(dca);
00534 if(_debuglevel>1) out() << "Muon pT after correction = " << it_muon->Pt() << endl;
00535 if(_debuglevel>1) out() << "Correcting track with SMT hits" << endl;
00536 }
00537 else if(_debuglevel>1 && smt>0) out() << "Not Correcting track since nSMT > 0 & _corr_pT_SMT is false" << endl;
00538 if(_corr_pT_noSMT && smt==0) {
00539 double dca = DCAfromBeamSpot(trk);
00540 if(_debuglevel>1) out() << "Muon pT before correction = " << it_muon->Pt() << endl;
00541 it_muon->CorrectPt(dca);
00542 if(_debuglevel>1) out() << "Correcting track with no SMT hits" << endl;
00543 if(_debuglevel>1) out() << "Muon pT after correction = " << it_muon->Pt() << endl;
00544 }
00545 }
00546 else cerr << "A muon with hasCentral() == 1, has no track, it will not be corrected" << endl;
00547 }
00548 }
00549 }
00550
00551
00552 if(_sort) {
00553 from.sort(from.begin(),from.end(),::moreThan) ;
00554 }
00555 }
00556
00557 bool MuonSelector::SetStdCuts(std::string path, std::string mu_file_name, std::string mu_name, int mu_ver, std::string trk_file_name, std::string trk_name, int trk_ver, std::string iso_file_name, std::string iso_name, int iso_ver) {
00558
00559 out() << endl;
00560 out() << "
00561 out() << "
00562 out() << "
00563
00564 ifstream mu_file;
00565 ifstream trk_file;
00566 ifstream iso_file;
00567
00568 string mu_path_file = path + mu_file_name;
00569 mu_file.open(mu_path_file.c_str());
00570
00571 string trk_path_file = path + trk_file_name;
00572 trk_file.open(trk_path_file.c_str());
00573
00574 string iso_path_file = path + iso_file_name;
00575 iso_file.open(iso_path_file.c_str());
00576
00577 if(!mu_file) {
00578 cerr << "In Muon Selector - cannot find the muon quality file " << mu_file_name.c_str() << endl;
00579 return false;
00580 }
00581 if(!trk_file) {
00582 cerr << "In Muon Selector - cannot find the muon trk quality file " << trk_file_name.c_str() << endl;
00583 return false;
00584 }
00585 if(!iso_file) {
00586 cerr << "In Muon Selector - cannot find the muon iso quality file " << iso_file_name.c_str() << endl;
00587 return false;
00588 }
00589
00590 //MuID uses beam spot for DCA correction
00591 _use_beam_spot = true;
00592
00593 string token = "";
00594 string qual = "";
00595 int ver;
00596
00597 while(!mu_file.eof()) {
00598 mu_file >> token;
00599
00600 if(token=="Version:") {
00601
00602 mu_file >> ver;
00603
00604 mu_file >> token;
00605 if(token!="Name:") cerr << "ERROR in setting muon cuts: name did not follow version in muon quality file" << endl;
00606 else {
00607 mu_file >> qual;
00608 if((ver==mu_ver)&&(qual==mu_name)) {
00609 //out() << "Using:" << endl;
00610 out() << "Muon ID Name: " << qual << endl;
00611 out() << "Version: " << ver << endl;
00612 break;
00613 }
00614 }
00615 }
00616 } //mu_file.eof()
00617
00618 if((ver!=mu_ver)||(qual!=mu_name)) {
00619 cerr << "Error no muon quality criteria found for " << mu_name << " and version " << mu_ver << endl;
00620 return false;
00621 }
00622 mu_file >> token;
00623 if(token=="END") {
00624 cerr << "No muon criteria in the muon quality file" << endl;
00625 return false;
00626 }
00627 while((token!="END")&&(!mu_file.eof())) {
00628 if(token=="MuonQuality:") {
00629 string q;
00630 mu_file >> q;
00631 if ( q == "LOOSE" || q == "Loose" || q == "loose")
00632 _qual = Loose ;
00633 else if ( q == "TIGHT" || q == "Tight" || q == "tight")
00634 _qual = Tight ;
00635 else if( q == "MEDIUM" || q == "Medium" || q == "medium")
00636 _qual = Medium ;
00637 else {
00638 err() << "MUON SELECTOR WARNING! "
00639 << "Parameter Quality = \"" << q
00640 << "\" is not recognized. "
00641 << "Must be loose, medium or tight. "
00642 << "Will use default quality MEDIUM."
00643 << endl ;
00644 _qual = Medium ;
00645 }
00646 }
00647 else {
00648 cerr << "Muon Quality not found" << endl;
00649 cerr << "ID file: " << mu_path_file << endl;
00650 return false;
00651 }
00652 mu_file >> token;
00653 if(token=="Nseg:") mu_file >> _nseg;
00654 else {
00655 cerr << "Nseg not found" << endl;
00656 cerr << "ID file: " << mu_path_file << endl;
00657 return false;
00658 }
00659 _useCentralMatched = true;
00660 mu_file >> token;
00661 if(token=="Cosmic:") mu_file >> _cosmicVeto;
00662 else {
00663 cerr << "Cosmic not found" << endl;
00664 cerr << "ID file: " << mu_path_file << endl;
00665 return false;
00666 }
00667 mu_file >> token;
00668 }
00669
00670 mu_file.close();
00671
00672 token = "";
00673 qual = "";
00674 ver = -1;
00675
00676 while(!trk_file.eof()) {
00677 trk_file >> token;
00678
00679 if(token=="Version:") {
00680
00681 trk_file >> ver;
00682
00683 trk_file >> token;
00684 if(token!="Name:") cerr << "ERROR in setting muon trk cuts: name did not follow version in muon trk quality file" << endl;
00685 else {
00686 trk_file >> qual;
00687 if((ver==trk_ver)&&(qual==trk_name)) {
00688
00689 out() << "Muon Track ID Name: " << qual << endl;
00690 out() << "Version: " << ver << endl;
00691 break;
00692 }
00693 }
00694 }
00695 }
00696
00697 if((ver!=trk_ver)||(qual!=trk_name)) {
00698 cerr << "Error no muon track criteria found for " << trk_name << " and version " << mu_ver << endl;
00699 return false;
00700 }
00701
00702 trk_file >> token;
00703 while( (!trk_file.eof()) && (token!="END") ) {
00704 if(token=="ChiSq:") trk_file >> _chisq;
00705 else {
00706 cerr << "ChiSq not found in Muon track file" << endl;
00707 cerr << "Track file: " << trk_path_file << endl;
00708 return false;
00709 }
00710 trk_file >> token;
00711 if(token=="DCASMT:") trk_file >> _dcaSMT;
00712 else {
00713 cerr << "DCASMT not found in Muon Track File" << endl;
00714 cerr << "Track file: " << trk_path_file << endl;
00715 return false;
00716 }
00717 trk_file >> token;
00718 if(token=="DCANoSMT:") trk_file >> _dcaNoSMT;
00719 else {
00720 cerr << "DCANoSMT not found in Muon Track File" << endl;
00721 cerr << "Track file: " << trk_path_file << endl;
00722 return false;
00723 }
00724 trk_file >> token;
00725 if(token=="nSMT:") trk_file >> _nSMT;
00726 else {
00727 cerr << "nSMT not found in Muon Track File" << endl;
00728 cerr << "Track file: " << trk_path_file << endl;
00729 return false;
00730 }
00731 trk_file >> token;
00732 }
00733
00734 trk_file.close();
00735
00736
00737
00738
00739 token = "";
00740 qual = "";
00741 ver = -1;
00742
00743 while(!iso_file.eof()) {
00744 iso_file >> token;
00745
00746 if(token=="Version:") {
00747
00748 iso_file >> ver;
00749
00750 iso_file >> token;
00751 if(token!="Name:") cerr << "ERROR in setting muon iso cuts: name did not follow version in muon iso quality file" << endl;
00752 else {
00753 iso_file >> qual;
00754 if((ver==iso_ver)&&(qual==iso_name)) {
00755
00756 out() << "Muon Isolation ID Name: " << qual << endl;
00757 out() << "Version: " << ver << endl;
00758 break;
00759 }
00760 }
00761 }
00762 }
00763
00764 if((ver!=iso_ver)||(qual!=iso_name)) {
00765 cerr << "Error no muon isolation criteria found for " << iso_name << " and version " << mu_ver << endl;
00766 return false;
00767 }
00768
00769 iso_file >> token;
00770 while( (!iso_file.eof()) && (token!="END") ) {
00771
00772 if(token=="etHalo:") iso_file >> _etHalo;
00773 else {
00774 cerr << "etHalo not found Muon Isolation file" << endl;
00775 cerr << "Isolation file: " << iso_path_file << endl;
00776 return false;
00777 }
00778 iso_file >> token;
00779
00780 if(token=="etHaloMin:") iso_file >> _etHaloMin;
00781 else {
00782 cerr << "etHaloMin not found Muon Isolation file" << endl;
00783 cerr << "Isolation file: " << iso_path_file << endl;
00784 return false;
00785 }
00786 iso_file >> token;
00787
00788 if(token=="etTrkCone:") iso_file >> _etTrkCone;
00789 else {
00790 cerr << "etTrkCone not found Muon Isolation file" << endl;
00791 cerr << "Isolation file: " << iso_path_file << endl;
00792 return false;
00793 }
00794 iso_file >> token;
00795
00796 if(token=="etTrkConeMin:") iso_file >> _etTrkConeMin;
00797 else {
00798 cerr << "etTrkConeMin not found Muon Isolation file" << endl;
00799 cerr << "Isolation file: " << iso_path_file << endl;
00800 return false;
00801 }
00802 iso_file >> token;
00803
00804 if(token=="etHaloScaled:") iso_file >> _etHaloScaled;
00805 else {
00806 cerr << "etHaloScaled not found Muon Isolation file" << endl;
00807 cerr << "Isolation file: " << iso_path_file << endl;
00808 return false;
00809 }
00810 iso_file >> token;
00811
00812 if(token=="etHaloScaledMin:") iso_file >> _etHaloScaledMin;
00813 else {
00814 cerr << "etHaloScaledMin not found Muon Isolation file" << endl;
00815 cerr << "Isolation file: " << iso_path_file << endl;
00816 return false;
00817 }
00818 iso_file >> token;
00819
00820 if(token=="etTrkConeScaled:") iso_file >> _etTrkConeScaled;
00821 else {
00822 cerr << "etTrkConeScaled not found Muon Isolation file" << endl;
00823 cerr << "Isolation file: " << iso_path_file << endl;
00824 return false;
00825 }
00826 iso_file >> token;
00827
00828 if(token=="etTrkConeScaledMin:") iso_file >> _etTrkConeScaledMin;
00829 else {
00830 cerr << "etTrkConeScaledMin not found Muon Isolation file" << endl;
00831 cerr << "Isolation file: " << iso_path_file << endl;
00832 return false;
00833 }
00834 iso_file >> token;
00835
00836 if(token=="MinDeltaR:") iso_file >> _MinDeltaR;
00837 else {
00838 cerr << "MinDeltaR not found Muon Isolation file" << endl;
00839 cerr << "Isolation file: " << iso_path_file << endl;
00840 return false;
00841 }
00842 iso_file >> token;
00843
00844 if(token=="MaxDeltaR:") iso_file >> _MaxDeltaR;
00845 else {
00846 cerr << "MaxDeltaR not found Muon Isolation file" << endl;
00847 cerr << "Isolation file: " << iso_path_file << endl;
00848 return false;
00849 }
00850 iso_file >> token;
00851 }
00852
00853 iso_file.close();
00854
00855
00856 _dcaSignif = -1.0;
00857
00858 return true;
00859 }
00860
00861 void MuonSelector::DisplayMuonQualityCuts() {
00862 out() << endl;
00863 out() << "Muon Defintion corresponds to cuts:" << endl;
00864 string quality = "";
00865 if(_qual==Tight) quality = "Tight";
00866 if(_qual==Medium) quality = "Medium";
00867 if(_qual==Loose) quality = "Loose";
00868 out() << "\tMuon Quality: " << quality << endl;
00869 out() << "\tNseg >= " << _nseg << endl;
00870 out() << "\tHasCentral: " << _useCentralMatched << endl;
00871 out() << "\tCosmicVeto: " << _cosmicVeto << endl;
00872 }
00873
00874 void MuonSelector::DisplayMuonTrackQualityCuts() {
00875 out() << endl;
00876 out() << "Muon Track Definition corresponds to cuts:" << endl;
00877 if(_chisq>0) out() << "\tChiSq < " << _chisq << endl;
00878 else out() << "\tNo ChiSq cut" << endl;
00879 if(_dcaSMT>0) out() << "\tDCA < " << _dcaSMT << " for nSMT > 0" << endl;
00880 else out() << "\tNo DCA cut for SMT > 0" << endl;
00881 if(_dcaNoSMT>0) out() << "\tDCA < " << _dcaNoSMT << " for nSMT = 0" << endl;
00882 else out() << "\tNo DCA cut for SMT = 0" << endl;
00883 if(_nSMT>0) out() << "\tN SMT hits >= " << _nSMT << endl;
00884 else out() << "\tNo SMT hit cut" << endl;
00885 }
00886
00887 void MuonSelector::DisplayMuonIsolationQualityCuts() {
00888 out() << endl;
00889 out() << "Muon Isolation Definition corresponds to cuts:" << endl;
00890 if(_etHaloScaled>0) out() << "\tCalorimeter Isolation < " << _etHaloScaled << endl;
00891 else out() << "\tNo calorimeter isolation cut" << endl;
00892 if(_etHaloScaledMin>0) out() << "\tCalorimeter Isolation > " << _etHaloScaled << endl;
00893 else out() << "\tNo calorimeter isolation min cut" << endl;
00894 if(_etTrkConeScaled>0) out() << "\tTrack Isolation < " << _etTrkConeScaled << endl;
00895 else out() << "\tNo track isolation cut" << endl;
00896 if(_etTrkConeScaledMin>0) out() << "\tTrack Isolation > " << _etTrkConeScaled << endl;
00897 else out() << "\tNo track isolation min cut" << endl;
00898 if(_MinDeltaR>0) out() << "\tDeltaR >= " << _MinDeltaR << endl;
00899 else out() << "\tNo DeltaR min cut" << endl;
00900 if(_MaxDeltaR>0) out() << "\tDeltaR < " << _MaxDeltaR << endl;
00901 else out() << "\tNo DeltaR max cut" << endl;
00902 }
00903
00904 bool MuonSelector::GetBeamSpot(FILE *fp, int run, float &bx, float &by, float &bsx, float &bsy) {
00905 rewind(fp);
00906 int cont=1, rn, en, pass=0;
00907 float xe, ye, xye;
00908 while(cont) {
00909 bx=by=bsx=bsy=0;
00910 fscanf(fp,"%i %i %f %f %f %f %f %f %f",&rn,&en,&bx,&by,&xe,&ye,&xye, &bsx, &bsy);
00911 if(rn==run) {
00912 cont=0;
00913 pass=1;
00914 }
00915 if(rn==0) cont=0;
00916 }
00917 if(!pass){
00918 bx = by = 0;
00919 bsx = bsy = 0;
00920
00921 std::cerr<<"Could not find run: "<<run<< " in beam spot file" << std::endl;
00922 return false;
00923
00924 }
00925
00926 return true;
00927
00928 }
00929
00930 double MuonSelector::DCAfromBeamSpot(const TMBTrack* muon_trk) {
00931 double z0 = muon_trk->z();
00932 double phi = muon_trk->phid();
00933 double beam_x = _bx + z0*_bsx;
00934 double beam_y = _by + z0*_bsy;
00935 double dca = muon_trk->rdca() - beam_x*sin(phi) + beam_y*cos(phi);
00936 if(_debuglevel>1) out() << "DCA = " << dca << endl;
00937 return dca;
00938 }
00939
00940 double MuonSelector::GetCorrectedPt(const TMBMuon& muon, double dca) {
00941 if(!muon.hasCentral()) {
00942 if(_debuglevel>0) cerr << "Cannot correct muon pT if muon is not central" << endl;
00943 return muon.Pt();
00944 }
00945 const TMBTrack* trk = muon.GetChargedTrack();
00946 if(trk) {
00947 double err_rqpt = trk->trerrs(4, 0);
00948 double err_rr = trk->trerrs(0, 0);
00949 float qopt = trk->qpt();
00950 qopt -= dca * err_rqpt / err_rr;
00951 double pTcorr = 1 / qopt;
00952 int q = 1;
00953 if(pTcorr < 0) {
00954 pTcorr *= -1;
00955 q = -1;
00956 }
00957 return pTcorr;
00958 }
00959 else cerr << "Cannot get corrected pT for muon with no matched track" << endl;
00960
00961 return muon.Pt();
00962 }
00963
00964
00965
00966 void MuonSelector::after(cafe::Collection<TMBMuon>& acc, cafe::Collection<TMBMuon>& rej) {
00967
00968 }
00969
00970
00971 }
00972 ClassImp(caf_util::MuonSelector) ;
00973
