CAF: JetSelector.cpp Source File

00001 #include "caf_util/JetSelector.hpp"
00002 #include "cafe/Config.hpp"
00003 #include "jetcorr/CalTool.hpp"
00004 #include <sstream>
00005 
00006 using namespace std ;
00007 using namespace cafe ;
00008  
00009 namespace{
00010     bool moreThan(const TMBJet& a, const TMBJet& b) {
00011         return a.Pt() > b.Pt();    
00012     }
00013 }
00014   
00015   
00016 namespace caf_util {
00017 
00018     JetSelector::JetSelector(const char *name) : cafe::SelectUserObjects<TMBJet>(name),
00019                                                  _nselected(0), _vertexVars ("fUniqueID"), 
00020                                                  _vars("_nchunks"),_varsglob("_runno") 
00021 
00022     {
00023         // get parameters from the configuration file
00024         cafe::Config config(name);
00025     
00026         _useStandardCuts = config.get("UseStandardID",true) ;    
00027         _useSmearedJets  = config.get("UseSmearedJets",false);
00028         _applyStandardIDCuts = config.get("UseStandardQualityCuts",false) ;    
00029         _applyStandardEMJetsRemoval = config.get("UseStandardEMJetsRemoval",false) ;    
00030         _applyStandardL1Conf = config.get("UseStandardL1Conf",false) ;    
00031         _requireVertexConf = config.get("RequireVertexConf",false) ;    
00032         _vertexBranch = config.get("VertexBranch","PrimaryVertex") ;    
00033         _jes = config.get("JES","JES") ;    
00034         _pTmin = config.get("pT",0.0) ;
00035         _eta = config.get("eta", 10.0);
00036         _detEta = config.get("detEta", 10.0);
00037         _sort  = config.get("Sort", 0);
00038 
00039         _EMFmin = config.get("EMFmin", 0.05);
00040         _EMFmax = config.get("EMFmax", 0.95);
00041         _CHFmin = config.get("CHFmin", -1.);
00042         _CHFmax = config.get("CHFmax", 0.4);
00043 
00044     
00045         _doL1cut = config.get("L1cut", false);
00046         _L1boundaryOffset = config.get("L1BoundOffs", 0.2);
00047         _f90noL1 = config.get("f90noL1", 0.5);
00048         _chfnoL1 = config.get("CHFnoL1", 0.5);
00049         _l1cutCC = config.get("l1cutCC", 0.4);
00050         _l1cutICD = config.get("l1cutICD", 0.4);
00051         _l1cutEC = config.get("l1cutEC", 0.4);
00052 
00053         _CPFmin = config.get ("CPFmin", -999.);
00054 
00055         _njets = config.get("nJets", 0);
00056         _njetsmax = config.get("nJetsMax", -1);
00057     };
00058   
00059     bool JetSelector::processEvent(cafe::Event &event)
00060     {
00061         _isMC=false;
00062         if(event.getMCEventInfo(_vars)!=NULL) _isMC=true;
00063         
00064         if (_useStandardCuts) _runnum = 0 ;
00065         else {
00066           const TMBGlobal* global = event.getGlobal (_varsglob) ;
00067           assert (global);
00068           _runnum = global->runno();
00069         }
00070 
00071         _nselected = 0 ;
00072 
00073         //get pointer to statistics collector
00074         event.get("StatPointer", _stat) ;
00075 
00076         //find PV0 if needed
00077         if (_requireVertexConf) {
00078           Collection<TMBPrimaryVertex> PVs (event.getCollection<TMBPrimaryVertex>(_vertexBranch,_vertexVars)) ;
00079           _PV0id = PVs.size() > 0 ? 0xffffff & PVs[0].GetUniqueID() : 0 ;
00080           if (debug()>5) // note the vz() will be junk unless partialReads are disabled
00081             for (unsigned int i=0; i<PVs.size(); ++i) 
00082               out() << i<<") uid: "<<PVs[i].GetUniqueID()<<" z: "<<PVs[i].vz()<<endl ;
00083         }
00084         
00085         SelectUserObjects<TMBJet>::processEvent(event) ;
00086 
00087         if (debug() >=2 ) 
00088           out() << " JET SELECTED: " << _nselected 
00089                 << std::endl ;
00090         
00091         for (int i = 1; i <= _njets; i++) {
00092             if (_nselected < i) return false ; 
00093             ostringstream st ;
00094             st << "Number of jets >= " << i ;
00095             _stat.EventSelected(st.str()) ;         
00096         }
00097 
00098         if (_njetsmax >= 0) {
00099           if (_nselected > _njetsmax) return false ; 
00100           ostringstream st ;
00101           st << "Number of jets <= " << _njetsmax ;
00102           _stat.EventSelected(st.str()) ;      
00103         }
00104 
00105         if (_nselected ==1 ) 
00106           event.tag("1jet") ;
00107         else if (_nselected > 1 ) 
00108           event.tag("2jet") ;
00109 
00110         return true  ;  
00111     };
00112   
00113     bool JetSelector::selectObject(const TMBJet &jet)
00114     {     
00115 
00116         // if JES was not set-up before in the sort metod,
00117         // set-up it here
00118         if (!_sort) { 
00119           SetJESType(jet) ;
00120           if (_useSmearedJets && _isMC) Smear(jet);
00121         }
00122 
00123         // jet ID cuts
00124         if( _useStandardCuts){
00125           if (!jet.isGood()) return false ;
00126           if (_njets>0) _stat.EventSelected("Jet ID") ;
00127     
00128           if (jet.isEM()) return false ;
00129           if (_njets>0) _stat.EventSelected("EM jets removal") ;
00130 
00131         } else {
00132 
00133           // non-standard cuts
00134 
00135           if (_applyStandardIDCuts) {
00136             if (!jet.isGood()) return false ;
00137             if (_njets>0) _stat.EventSelected("Jet ID") ;
00138           } else {
00139             if( jet.emf()<=_EMFmin ) return false;
00140             { ostringstream st ;
00141               st << "EM Fraction > " << _EMFmin ;
00142             if (_njets>0) _stat.EventSelected(st.str()) ; }
00143             
00144             if( jet.emf()>=_EMFmax ) return false ;
00145             { ostringstream st ;
00146               st << "EM Fraction < " << _EMFmax ;
00147               if (_njets>0) _stat.EventSelected(st.str()) ; }
00148             
00149             if( jet.chf()<=_CHFmin ) return false ;
00150             { ostringstream st ;
00151               st << "CH Fraction > " << _CHFmin ;
00152               if (_njets>0) _stat.EventSelected(st.str()) ; }
00153             
00154             if( jet.chf()>=_CHFmax ) return false ;
00155             { ostringstream st ;
00156               st << "CH Fraction < " << _CHFmax ;
00157               if (_njets>0) _stat.EventSelected(st.str()) ; }
00158           }
00159           
00160           if (_applyStandardEMJetsRemoval) {
00161             if (jet.isEM()) return false ;
00162             if (_njets>0) _stat.EventSelected("EM jets removal") ;
00163           }
00164 
00165           //L1 confirmation
00166           if (_applyStandardL1Conf) {
00167             if (!jet.isL1Conf()) return false ;
00168             if (_njets>0) _stat.EventSelected("Jet L1 confirmation") ;
00169           } else {
00170             if(!_isMC && _doL1cut && jet.L1_energy()<55.) {
00171               if(fabs(jet.detEta()*.1)>MaxAbsEtaOfL1CalCoverage(_runnum) - _L1boundaryOffset) {
00172                 float f90=0.;
00173                 if(jet.nItems()!=0) {
00174                 f90=(float)jet.n90()/(float)jet.nItems();
00175                 if( !((f90<_f90noL1)||(jet.chf()<_chfnoL1)) ) 
00176                   return false ;
00177                 }
00178               } else { 
00179                 float ratio=0.;
00180                 bool CC=false, ICD=false, EC=false;
00181                 if(fabs(jet.detEta()*0.1)<=0.8) CC=true;
00182                 if(fabs(jet.detEta()*0.1)>0.8) ICD=true;
00183                 if(fabs(jet.detEta()*0.1)>1.5) {ICD=false; EC=true;};
00184                 
00185                 TMBJet _jet = jet;
00186                 _jet.ActAsUnCorrected();
00187                 ratio = _jet.L1_energy()/(jet.Pt()*(1.-jet.chf()));
00188                 
00189                 if(CC&&(ratio<_l1cutCC)) return false ;
00190                 if(ICD&&(ratio<_l1cutICD)) return false ;
00191                 if(EC&&(ratio<_l1cutEC)) return false ;
00192               }
00193             }
00194             if(!_isMC && _doL1cut && _njets>0)
00195               _stat.EventSelected("Jet L1 confirmation") ;
00196           }
00197         } // non-standard ID
00198 
00199         float etad = CalTool::EtaDToEtaDet(jet.detEta());
00200         if (fabs(etad) >= _detEta)  return false ;
00201         { ostringstream st ;
00202           st << "Jet |detEta| < " << (float) _detEta ;
00203           if (_njets>0 && _detEta < 5.0) _stat.EventSelected(st.str()) ; }
00204 
00205         if( fabs(jet.Eta()) >= _eta) return false ;
00206         { ostringstream st ;
00207           st << "Jet eta < " << (float) _eta ;
00208           if (_njets>0 && _eta < 5.0) _stat.EventSelected(st.str()) ; }
00209         
00210         if((jet.Pt()<=0)||(jet.E()<=0)) {
00211           err() << "caf_util/JetSelector: pT or energy <= 0!"
00212                 << " Probably your jet is not good or is an EM jet."
00213                 << endl ;
00214           return false ;
00215         }
00216 
00217         if(jet.Pt() < _pTmin) return false ;
00218         {
00219           ostringstream st ;
00220           st << "Jet pT >= " << _pTmin << " GeV" ;
00221           if (_njets>0) _stat.EventSelected(st.str()) ;
00222         }
00223 
00224         if (_requireVertexConf) {
00225           const TMBVertex *matchingPV = jet.cpfVertex () ;
00226           if (matchingPV == 0) return false;
00227 
00228           UInt_t mPVuid = 0xffffff & matchingPV->GetUniqueID() ;
00229           if (debug()>9) out() << "z: "<<matchingPV->vz()<<" mPVuid: "<<mPVuid<<" =? "<<_PV0id<<endl ;
00230           if (mPVuid != _PV0id) return false ;
00231 
00232           if (jet.cpf0() < _CPFmin) return false ;
00233 
00234           ostringstream st ;
00235           st << "Jet matchs P.V.[0] with CPF >= " << _CPFmin ;
00236           if (_njets>0) _stat.EventSelected(st.str()) ;
00237         }
00238         
00239         
00240         _nselected++ ;
00241         
00242     if (debug()>2) 
00243       cout << "PROCESSOR \"" << name() 
00244            << "\" SELECTED OBJECT: pT = " << jet.Pt() 
00245            << " eta = " << jet.Eta() 
00246            << " phi = " << jet.Phi() 
00247            << endl ;
00248 
00249         return true ;
00250         
00251     };
00252   
00253     float JetSelector::MaxAbsEtaOfL1CalCoverage(int runnum)
00254     {
00255         if(runnum<=157712) return 0.8;
00256         if(runnum<=172358) return 2.4;
00257         return 1000.;
00258     };
00259     
00260     void JetSelector::before(Collection<TMBJet>& from) 
00261     {
00262         if(_sort) {
00263             for(Collection<TMBJet>::const_iterator it = from.begin() ;
00264                 it != from.end() ; 
00265                 it++ ) {
00266                 SetJESType(*it) ;
00267                 if (_useSmearedJets && _isMC) Smear(*it);
00268             }
00269             from.sort(from.begin(),from.end(),::moreThan) ;
00270         }
00271     }
00272 
00274     void JetSelector::Smear(const TMBJet& jet) {
00275        if (_jes == "Uncorrected") {
00276          throw std::runtime_error(std::string("JetSelector[" + name() + "] : You want to use the jet smearing," +
00277                                   " you need to specify after which JES correction " +
00278                                   " using the field JES in your config file"));
00279        }
00280        if (jet.CurrentlyActAs() == TMBJet::kSmeared || jet.CurrentlyActAs() == TMBJet::kSmearedMU) return;
00281        if (_jes == "JES" || _jes == "JESShiftedPlus" || _jes == "JESShiftedMinus") jet.ActAsSmeared();
00282        if (_jes == "JESMU" || _jes == "JESMUShiftedPlus" || _jes == "JESMUShiftedMinus") jet.ActAsSmearedMU();
00283     }
00284 
00286     void JetSelector::SetJESType(const TMBJet& jet) {
00287         if (_jes == "UnCorrected") {
00288             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected)
00289                 err() << "JET SELECTION WARNING! "
00290                       << "Jet was asked to act as uncorrected with JES, "
00291                       << " but type for JES was setup before with value "
00292                       << jet.CurrentlyActAs() << "." << endl ;
00293             jet.ActAsUnCorrected() ;
00294       
00295         } else if (_jes == "JESMU") {
00296             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00297                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrected &&
00298                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedPlus &&
00299                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedMinus &&
00300                 jet.CurrentlyActAs() != TMBJet::kSmearedMU)
00301                 err() << "JET SELECTION WARNING! "
00302                       << "Jet was asked to act as JESMU corrected, "
00303                       << " but type for JES was setup before with value "
00304                       << jet.CurrentlyActAs() << "." << endl ;
00305             jet.ActAsJESMUCorrected() ;
00306       
00307         } else if (_jes == "JESMUShiftedPlus") {
00308             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00309                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrected &&
00310                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrectedShiftedPlus &&
00311                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrectedShiftedMinus &&
00312                 jet.CurrentlyActAs() != TMBJet::kSmearedMU)
00313                 err() << "JET SELECTION WARNING! "
00314                       << "Jet was asked to act as JESMU corrected, "
00315                       << " but type for JES was setup before with value "
00316                       << jet.CurrentlyActAs() << "." << endl ;
00317             jet.ActAsJESMUCorrectedShiftedPlus() ;
00318         } else if (_jes == "JESMUShiftedMinus") {
00319             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00320                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrected &&
00321                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrectedShiftedPlus &&
00322                 jet.CurrentlyActAs() != TMBJet::kJESMUCorrectedShiftedMinus)            
00323                 err() << "JET SELECTION WARNING! "
00324                       << "Jet was asked to act as JESMU corrected, "
00325                       << " but type for JES was setup before with value "
00326                       << jet.CurrentlyActAs() << "." << endl ;
00327             jet.ActAsJESMUCorrectedShiftedMinus() ;
00328         } else if (_jes == "JES") {
00329             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00330                 jet.CurrentlyActAs() != TMBJet::kJESCorrected && 
00331                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedPlus &&
00332                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedMinus && 
00333                 jet.CurrentlyActAs() != TMBJet::kSmeared)
00334               err() << "JET SELECTION WARNING! "
00335                       << "Jet was asked to act as JES corrected, "
00336                       << " but type for JES was setup before with value "
00337                       << jet.CurrentlyActAs() << "." << endl ;  
00338             jet.ActAsJESCorrected() ;
00339         } else if (_jes == "JESShiftedPlus") {
00340             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00341                 jet.CurrentlyActAs() != TMBJet::kJESCorrected &&
00342                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedPlus &&
00343                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedMinus &&
00344                 jet.CurrentlyActAs() != TMBJet::kSmeared)
00345                 err() << "JET SELECTION WARNING! "
00346                       << "Jet was asked to act as JES corrected, "
00347                       << " but type for JES was setup before with value "
00348                       << jet.CurrentlyActAs() << "." << endl ;
00349             jet.ActAsJESCorrectedShiftedPlus() ;
00350         } else if (_jes == "JESShiftedMinus") {
00351             if (jet.CurrentlyActAs() != TMBJet::kUnCorrected &&
00352                 jet.CurrentlyActAs() != TMBJet::kJESCorrected &&
00353                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedPlus &&
00354                 jet.CurrentlyActAs() != TMBJet::kJESCorrectedShiftedMinus &&
00355                 jet.CurrentlyActAs() != TMBJet::kSmeared)
00356                 err() << "JET SELECTION WARNING! "
00357                       << "Jet was asked to act as JES corrected, "
00358                       << " but type for JES was setup before with value "
00359                       << jet.CurrentlyActAs() << "." << endl ;
00360             jet.ActAsJESCorrectedShiftedMinus() ;
00361         } else {
00362             err() << "JET SELECTION WARNING! "
00363                   << "JES type \"" << _jes << "\" is not defined! "
00364                   << "Default JES will be applied!" << endl ;
00365             jet.ActAsJESCorrected() ;
00366         }
00367         return ;
00368     }
00369 }
00370 ClassImp(caf_util::JetSelector) ;
00371