Go to the documentation of this file.00001 #include "ForIA/AnalysisTools/JetSelection.hh"
00002
00003
00004 #include "ForIA/Units.hh"
00005 #include <cmath>
00006 #include <stdexcept>
00007
00008 namespace ForIA {
00009
00010 JetSelection::JetSelection():
00011 FourVectorSelection(),
00012 m_jetAlgKey(Jet::ANTI_KT, 0.4, Jet::EM_TOPO),
00013 m_haveSeteMin(false),
00014 m_haveSetetMin(false),
00015 m_haveSetptMin(false),
00016 m_haveSetetaMin(false),
00017 m_haveSetetaMax(false),
00018 m_haveSetetMinJet1(false),
00019 m_haveSetetMinJet2(false),
00020 m_haveSetetaMinJet1(false),
00021 m_haveSetetaMaxJet1(false),
00022 m_haveSetetaMinJet2(false),
00023 m_haveSetetaMaxJet2(false),
00024 m_haveSetDeltaPhi(false),
00025 m_haveSetJet2ETOverJet1ET(false),
00026 m_haveSetApplyJESUncertCorr(false),
00027 m_haveSetApplyJERUncertCorr(false)
00028 {
00029 std::cout<<"JetSelector initialising JESTool with: AntiKt4TopoJetsEM"<<std::endl;
00030 m_JESTool = new JESUncertaintyProvider("AntiKt4TopoJetsEM");
00031 m_JESTool->init();
00032 std::cout<<"JetSelector initialising JERTool with: AntiKt4TopoJES"<<std::endl;
00033
00034 m_JERTool = new JERProvider("AntiKt4TopoJES","Truth", "./JERProviderPlots.root");
00035 m_JERTool->init();
00036
00037 }
00038
00039 const JetVector &JetSelection::acceptedJets() const{
00040 if(m_freshEvent) fillVectors();
00041 return m_acceptedJets;
00042 }
00043 const JetVector &JetSelection::rejectedJets() const{
00044 if(m_freshEvent) fillVectors();
00045 return m_rejectedJets;
00046 }
00047
00048 bool JetSelection::setSelection(DefinedSelection selection){
00049 if(!m_canSetCuts) return false;
00050
00051 m_haveSeteMin = false;
00052 m_haveSetetMin = false;
00053 m_haveSetptMin = false;
00054 m_haveSetetaMin= false;
00055 m_haveSetetaMax= false;
00056
00057 switch (selection) {
00058 case LOOSE:
00059 setSelection(LOOSE);
00060 setETMinJet1(10.0 * GeV);
00061 break;
00062
00063 case MEDIUM:
00064 setSelection(MEDIUM);
00065 setETMinJet1(10.0 * GeV);
00066 setETMinJet2(10.0 * GeV);
00067 break;
00068
00069 case TIGHT:
00070 setSelection(TIGHT);
00071 setEMin(20.0 * GeV);
00072 break;
00073
00074 }
00075
00076 return true;
00077 }
00078
00079 void JetSelection::fillVectors() const{
00080
00081
00082
00083 m_freshEvent = false;
00084
00085 m_acceptedJets.clear();
00086 m_passedFourVectors.clear();
00087 m_rejectedJets.clear();
00088 m_rejectedFourVectors.clear();
00089
00090 if(m_event == 0) return;
00091
00092 vector<JetPtr> input;
00093
00094 JetVector toCopy = m_event->jets(m_jetAlgKey);
00095
00096 for(JetVector::const_iterator jet = toCopy.begin();
00097 jet != toCopy.end(); ++jet){
00098 input.push_back(JetPtr(new Jet(*(*jet))));
00099 }
00100
00101 if(input.size() == 0 ) return;
00102
00103
00104 ClusterVector clusters = m_event->clusters();
00105 bool isMC = clusters.at(0)->isMC();
00106
00107 if(!isMC){
00108 bool eventContainsBadJet = false;
00109 for(vector<JetPtr>::const_iterator jet = input.begin();
00110 jet != input.end(); ++jet){
00111 if( (*jet)->isBadLoose() ){
00112 eventContainsBadJet = true;
00113 break;
00114 }
00115 }
00116 if(eventContainsBadJet) return;
00117 }
00118
00119 bool passedCuts = true;
00120
00121 if(passedCuts){
00122 for(vector<JetPtr>::const_iterator jet = input.begin();
00123 jet != input.end(); ++jet){
00124
00125 bool passedCuts = true;
00126
00127
00128
00129
00130
00131
00132 if( !isMC && (*jet)->isBadLoose() ) continue;
00133
00134
00135 if(m_haveSetApplyJESUncertCorr && m_JESTool!=0){
00136 if( (*jet)->PT()>15000. && fabs((*jet)->eta())<4.5 ){
00137 double JESCorr = m_sign * m_JESTool->getRelUncert( (*jet)->PT(), (*jet)->eta(), JESUncertaintyProvider::ALL, 1);
00138 std::cout<<"JESCorr: " << JESCorr <<std::endl;
00139 std::cout<<"Before scaling: "<< (*jet)->E() <<" et: " << (*jet)->ET() <<" eta: "<< (*jet)->eta() <<" phi: "<< (*jet)->phi() <<std::endl;
00140 (*jet)->scaleEnergy( JESCorr, true );
00141 std::cout<<"After scaling: "<< (*jet)->E() <<" et: " << (*jet)->ET() <<" eta: "<< (*jet)->eta() <<" phi: "<< (*jet)->phi() <<std::endl;
00142 std::cout<<"Done scaling: -------------------------------- "<<std::endl;
00143 }
00144 }
00145
00146 if(m_haveSetApplyJERUncertCorr && m_JERTool!=0){
00147 double sigma = m_JERTool->getSigma( (*jet)->PT()/1000., (*jet)->eta() );
00148 double uncert= m_JERTool->getUncert( (*jet)->PT()/1000., (*jet)->eta() );
00149 double smear = sqrt( pow( sigma+uncert, 2) - pow( sigma, 2) );
00150 double JERCorr = gRandom->Gaus( 1., smear );
00151 (*jet)->scaleEnergy( JERCorr, false);
00152 }
00153
00154
00155 if(m_haveSeteMin && passedCuts){
00156 passedCuts = ((*jet)->E() > m_eMin);
00157 }
00158 if (m_haveSetetMin && passedCuts){
00159 passedCuts = ((*jet)->ET() >= m_etMin);
00160 }
00161 if (m_haveSetptMin && passedCuts){
00162 passedCuts = ((*jet)->PT() >= m_ptMin);
00163 }
00164 if(m_haveSetetaMin && passedCuts){ passedCuts = (fabs((*jet)->eta()) >= m_etaMin);}
00165 if(m_haveSetetaMax && passedCuts){ passedCuts = (fabs((*jet)->eta()) <= m_etaMax);}
00166
00167 if(passedCuts){
00168 m_acceptedJets.push_back( *jet );
00169 m_passedFourVectors.push_back( *jet );
00170 } else {
00171 m_rejectedJets.push_back( *jet );
00172 m_rejectedFourVectors.push_back( *jet );
00173 }
00174 }
00175 }
00176
00177
00178 return;
00179 }
00180
00181 bool JetSelection::setETMin(double etMin){
00182 if(!m_canSetCuts) return false;
00183 m_haveSetetMin = true;
00184 m_etMin = etMin;
00185 return true;
00186 }
00187
00188 bool JetSelection::setETMinJet1(double etMinJet1){
00189 if(!m_canSetCuts) return false;
00190 m_haveSetetMinJet1 = true;
00191 m_etMinJet1 = etMinJet1;
00192 return true;
00193 }
00194
00195 bool JetSelection::setETMinJet2(double etMinJet2){
00196 if(!m_canSetCuts) return false;
00197 m_haveSetetMinJet2 = true;
00198 m_etMinJet2 = etMinJet2;
00199 return true;
00200 }
00201
00202 bool JetSelection::setEtaMinJet1(double etaMinJet1){
00203 if(!m_canSetCuts) return false;
00204 m_haveSetetaMinJet1 = true;
00205 m_etaMinJet1 = etaMinJet1;
00206 return true;
00207 }
00208
00209 bool JetSelection::setEtaMaxJet1(double etaMaxJet1){
00210 if(!m_canSetCuts) return false;
00211 m_haveSetetaMaxJet1 = true;
00212 m_etaMaxJet1 = etaMaxJet1;
00213 return true;
00214 }
00215
00216 bool JetSelection::setEtaMinJet2(double etaMinJet2){
00217 if(!m_canSetCuts) return false;
00218 m_haveSetetaMinJet2 = true;
00219 m_etaMinJet2 = etaMinJet2;
00220 return true;
00221 }
00222
00223 bool JetSelection::setEtaMaxJet2(double etaMaxJet2){
00224 if(!m_canSetCuts) return false;
00225 m_haveSetetaMaxJet2 = true;
00226 m_etaMaxJet2 = etaMaxJet2;
00227 return true;
00228 }
00229
00230 bool JetSelection::setDeltaPhi( double deltaPhi ){
00231 if(!m_canSetCuts) return false;
00232 m_haveSetDeltaPhi = true;
00233 m_deltaPhi = deltaPhi;
00234 return true;
00235 }
00236
00237 bool JetSelection::setJet2ETOverJet1ET( double jet2ETOverjet1ET ){
00238 if(!m_canSetCuts) return false;
00239 m_haveSetJet2ETOverJet1ET = true;
00240 m_jet2ETOverjet1ET = jet2ETOverjet1ET;
00241 return true;
00242 }
00243
00244 bool JetSelection::setJESUncertCorr(int sign){
00245 if(!m_canSetCuts) return false;
00246 m_haveSetApplyJESUncertCorr = true;
00247 m_sign = sign;
00248 return true;
00249 }
00250
00251 bool JetSelection::setJERUncertCorr(bool jer){
00252 if(!m_canSetCuts) return false;
00253 m_haveSetApplyJERUncertCorr = jer;
00254 return true;
00255 }
00256 bool JetSelection::setPTMin(double ptMin){
00257 if(!m_canSetCuts) return false;
00258 m_haveSetptMin = true;
00259 m_ptMin = ptMin;
00260 return true;
00261 }
00262
00263 bool JetSelection::setEMin(double eMin){
00264 if(!m_canSetCuts) return false;
00265 m_haveSeteMin = true;
00266 m_eMin = eMin;
00267 return true;
00268 }
00269
00270 bool JetSelection::setEtaMin(double etaMin){
00271 if(!m_canSetCuts) return false;
00272 m_haveSetetaMin = true;
00273 m_etaMin = etaMin;
00274 return true;
00275 }
00276
00277 bool JetSelection::setEtaMax(double etaMax){
00278 if(!m_canSetCuts) return false;
00279 m_haveSetetaMax = true;
00280 m_etaMax = etaMax;
00281 return true;
00282 }
00283
00284 }
