001 #include "JetTagTools/EmulTag.h"
002
003 #include "GaudiKernel/IToolSvc.h"
004 #include "GaudiKernel/ITHistSvc.h"
005 #include <CLHEP/Random/Randomize.h>
006
007 #include "JetEvent/Jet.h"
008 #include "JetTagInfo/BaseTagInfo.h"
009 #include "JetTagInfo/AtlfInfo.h"
010
011 #include "TH1F.h"
012 #include "TH2F.h"
013
014 namespace Analysis {
015
016 EmulTag::EmulTag(const std::string& t, const std::string& n, const IInterface* p)
017 : AthAlgTool(t,n,p) {
018 declareInterface<ITagTool>(this);
019 declareProperty("Mode", m_mode = 2);
020 declareProperty("ParamSet", m_atlfBNSet = 10);
021 declareProperty("Tagger", m_algoType);
022 m_algoType.push_back("SV1IP3D");
023 m_algoType.push_back("6d");
024 }
025
026 EmulTag::~EmulTag() { }
027
028 StatusCode EmulTag::initialize() {
029
030 m_imod = m_mode;
031 if (m_imod != 1 && m_imod != 2) {
032 ATH_MSG_ERROR("Unknown mode, use emulation instead !");
033 m_imod = 2;
034 }
035
036 if (m_imod == 2) {
037
038 for (unsigned int i = 0; i < m_algoType.size(); i++) m_ialgoType.push_back(m_algoType[i]);
039 if (m_ialgoType.size() != 2) {
040 ATH_MSG_ERROR("Tagger should be a 2-vector of strings, use SV1+IP3D emulation !");
041 m_ialgoType.clear();
042 m_ialgoType.push_back("SV1IP3D");
043 m_ialgoType.push_back("6d");
044 }
045
046 ATH_MSG_INFO("Initializing " << name() << " Tagger = " << m_ialgoType[0]);
047 std::string type = m_ialgoType[1];
048 StatusCode sc = initHistos(type);
049 if (sc.isFailure()) {
050 ATH_MSG_ERROR("Unable to initialize weight histograms");
051 return StatusCode::FAILURE;
052 }
053 } else {
054 ATH_MSG_INFO("Initializing " << name() << " with parameterisation " << m_atlfBNSet);
055 StatusCode sc = initParam(m_atlfBNSet);
056 if (sc.isFailure()) {
057 ATH_MSG_ERROR("Unable to initialize parameterisation");
058 return StatusCode::FAILURE;
059 }
060 }
061
062 return StatusCode::SUCCESS;
063 }
064
065
066 StatusCode EmulTag::finalize() {
067 return StatusCode::SUCCESS;
068 }
069
070 void EmulTag::tagJet(Jet& jetToTag) {
071
072
073 double ptjet = jetToTag.et()/1.e3;
074 double etajet = jetToTag.eta();
075 int ip = findEt(ptjet);
076 int ie = findEta(etajet);
077 int pdg = abs(jetToTag.pdgId());
078 double dRb = 999., dRc = 999., dRt = 999.;
079 const Analysis::AtlfInfo* ainfo = jetToTag.tagInfo<Analysis::AtlfInfo>("AtlfInfo");
080 if (ainfo) {
081 dRb = ainfo->deltaRMinTo("B");
082 dRc = ainfo->deltaRMinTo("C");
083 dRt = ainfo->deltaRMinTo("T");
084 }
085
086 bool isB = false;
087 double epsilon = 0.6;
088 double randnum = 0.; if (m_imod == 1) randnum = RandFlat::shoot();
089 double w = -20.;
090 bool firstBin = true;
091
092
093 if (fabs(etajet) <= 2.5) {
094 if (pdg == 5) {
095 if (m_imod == 2) {
096 if (m_hWb[ie][ip]) {
097 w = m_Wb[ie][ip]->GetRandom();
098 if (m_Wb[ie][ip]->GetXaxis()->FindBin(w) > 1) firstBin = false;
099 }
100 } else {
101 epsilon = m_EpsB->GetBinContent(ip+1,ie+1);
102 }
103 } else {
104 int ifl = 0;
105 int iso = 1;
106 if (pdg == 0) {
107 if (dRb < 0.8 || dRc < 0.8 || dRt < 0.8) iso = 0;
108 }
109 if (pdg == 4) {
110 ifl = 1;
111 if (dRb < 0.8) iso = 0;
112 }
113 if (pdg == 15) {
114 ifl = 2;
115 if (dRb < 0.8 || dRc < 0.8) iso = 0;
116 }
117 if (m_imod == 2) {
118 if (m_Wl[iso][ifl][ie][ip]) {
119 w = m_Wl[iso][ifl][ie][ip]->GetRandom();
120 if (m_Wl[iso][ifl][ie][ip]->GetXaxis()->FindBin(w) > 1) firstBin = false;
121 }
122 } else {
123 if (ifl == 0) {
124 if (iso == 1) epsilon = 1./m_RejLightNP->GetBinContent(ip+1,ie+1);
125 if (iso == 0) epsilon = 1./m_RejLightP->GetBinContent(ip+1,ie+1);
126 } else if (ifl == 4) {
127 epsilon = 1./m_RejC->GetBinContent(ip+1,ie+1);
128 } else if (ifl == 15) {
129 epsilon = 1./m_RejTau->GetBinContent(ip+1,ie+1);
130 }
131 }
132 }
133 }
134
135 if (m_imod == 2) {
136
137 BaseTagInfo* info = new BaseTagInfo(m_ialgoType[0]);
138 std::vector<double> like;
139 like.push_back(1.);
140 double pu = firstBin ? 1.e9 : exp(-w);
141 like.push_back(pu);
142 info->setTagLikelihood(like);
143 info->makeValid();
144 jetToTag.addInfo(info);
145 } else {
146 if (randnum < epsilon) isB = true;
147 jetToTag.setCombinedLikelihood( std::vector<double>( 1, isB
148 ? 1.
149 : 0. ) );
150
151 }
152 }
153
154 void EmulTag::finalizeHistos() { }
155
156 StatusCode EmulTag::initHistos(std::string algN) {
157
158 StatusCode sc = StatusCode::SUCCESS;
159
160 ITHistSvc* histoSvc;
161 sc = service("THistSvc", histoSvc);
162 if (sc.isFailure()) {
163 ATH_MSG_FATAL("THistSvc not found!");
164 return StatusCode::FAILURE;
165 }
166
167 std::string fullPath = "/WeightFile/";
168 std::string flav[4] = {"bw","uw","cw","tw"};
169 for (int fl = 0; fl < 4; fl++) {
170 for (int ie = 0; ie < m_etabin-1; ie++) {
171 for (int ip = 0; ip < m_ptbin-1; ip++) {
172 std::ostringstream oe, op;
173 oe << ie;
174 op << ip;
175 std::string hName = fullPath + flav[fl] + algN + "e" + oe.str() + "p" + op.str();
176 if (fl == 0) {
177 sc = histoSvc->regHist(hName);
178 sc = histoSvc->getHist(hName,m_Wb[ie][ip]);
179 } else {
180 sc = histoSvc->regHist(hName+"i");
181 sc = histoSvc->getHist(hName+"i",m_Wl[1][fl-1][ie][ip]);
182 sc = histoSvc->regHist(hName+"ni");
183 sc = histoSvc->getHist(hName+"ni",m_Wl[0][fl-1][ie][ip]);
184 }
185 if (sc.isFailure()) {
186 ATH_MSG_FATAL("Missing weight histograms " << hName);
187 return StatusCode::FAILURE;
188 }
189 }
190 }
191 }
192
193 for (int fl = 0; fl < 4; fl++) {
194 for (int ie = 0; ie < m_etabin-1; ie++) {
195 for (int ip = 0; ip < m_ptbin-1; ip++) {
196 if (fl == 0) {
197 m_hWb[ie][ip] = true;
198 if (m_Wb[ie][ip]->Integral() <= 0.) {
199 ATH_MSG_WARNING("histo " << flav[fl] << " " << algN << " " << ie << " " << ip << " empty");
200 m_hWb[ie][ip] = false;
201 }
202 } else {
203 m_hWl[0][fl-1][ie][ip] = true;
204 m_hWl[1][fl-1][ie][ip] = true;
205 if (m_Wl[0][fl-1][ie][ip]->Integral() <= 0.) {
206 ATH_MSG_WARNING("histo " << flav[fl] << " " << algN << " " << ie << " " << ip << " empty");
207 m_hWl[0][fl-1][ie][ip] = false;
208 }
209 if (m_Wl[1][fl-1][ie][ip]->Integral() <= 0.) {
210 ATH_MSG_WARNING("histo " << flav[fl] << " " << algN << " " << ie << " " << ip << " empty");
211 m_hWl[1][fl-1][ie][ip] = false;
212 }
213 }
214 }
215 }
216 }
217 return sc;
218 }
219
220 StatusCode EmulTag::initParam(unsigned int iset) {
221
222 StatusCode sc = StatusCode::SUCCESS;
223
224 double m_epsilonBjet = 0.6;
225 unsigned int ialgo = iset;
226 std::string algName = "IP2D";
227 if ((ialgo > 7 && ialgo) < 15 || (ialgo > 21 && ialgo < 29)) algName = "SV1+IP3D";
228 if (ialgo == 100) {
229 algName = "Canonical";
230 }
231
232 if (m_atlfBNSet <= 14 && ialgo > 7) ialgo -= 7;
233 m_epsilonBjet = 0.5 + (ialgo - 1)*0.05;
234 msg(MSG::INFO) << "Running with " << algName;
235 if (m_atlfBNSet < 15 || m_atlfBNSet == 100) msg(MSG::INFO) << " at eps_b = " << m_epsilonBjet << endreq;
236 if (m_atlfBNSet > 14) msg(MSG::INFO) << " with a fixed cut" << endreq;
237
238
239 ITHistSvc* histoSvc;
240 sc = service("THistSvc", histoSvc);
241 if (sc.isFailure()) {
242 ATH_MSG_FATAL("THistSvc not found!");
243 return StatusCode::FAILURE;
244 }
245
246 std::string fullPath = "/ParamFile/";
247 std::ostringstream o;
248 o << iset;
249 std::string hisname = fullPath+"BEfficiency"+o.str();
250 sc = histoSvc->regHist(hisname);
251 sc = histoSvc->getHist(hisname,m_EpsB);
252 hisname = fullPath+"RejectionFactorPU"+o.str();
253 sc = histoSvc->regHist(hisname);
254 sc = histoSvc->getHist(hisname,m_RejLightP);
255 hisname = fullPath+"RejectionFactorNPU"+o.str();
256 sc = histoSvc->regHist(hisname);
257 sc = histoSvc->getHist(hisname,m_RejLightNP);
258 hisname = fullPath+"RejectionFactorC"+o.str();
259 sc = histoSvc->regHist(hisname);
260 sc = histoSvc->getHist(hisname,m_RejC);
261 hisname = fullPath+"RejectionFactorTau"+o.str();
262 sc = histoSvc->regHist(hisname);
263 sc = histoSvc->getHist(hisname,m_RejTau);
264 if (sc.isFailure()) {
265 ATH_MSG_FATAL("Missing parameterisation histograms ");
266 return StatusCode::FAILURE;
267 }
268
269 return sc;
270 }
271
272 int EmulTag::findEta(double eta) {
273 if (fabs(eta) < m_ebin[1] && fabs(eta) >= m_ebin[0]) return 0;
274 if (fabs(eta) < m_ebin[2] && fabs(eta) >= m_ebin[1]) return 1;
275 if (fabs(eta) < m_ebin[3] && fabs(eta) >= m_ebin[2]) return 2;
276 if (fabs(eta) < m_ebin[4] && fabs(eta) >= m_ebin[3]) return 3;
277 if (fabs(eta) <= m_ebin[5] && fabs(eta) >= m_ebin[4]) return 4;
278 return 0;
279 }
280 int EmulTag::findEt(double et) {
281 if (fabs(et) < m_pbin[1] && fabs(et) >= m_pbin[0]) return 0;
282 if (fabs(et) < m_pbin[2] && fabs(et) >= m_pbin[1]) return 1;
283 if (fabs(et) < m_pbin[3] && fabs(et) >= m_pbin[2]) return 2;
284 if (fabs(et) < m_pbin[4] && fabs(et) >= m_pbin[3]) return 3;
285 if (fabs(et) < m_pbin[5] && fabs(et) >= m_pbin[4]) return 4;
286 if (fabs(et) < m_pbin[6] && fabs(et) >= m_pbin[5]) return 5;
287 if (fabs(et) < m_pbin[7] && fabs(et) >= m_pbin[6]) return 6;
288 if (fabs(et) < m_pbin[8] && fabs(et) >= m_pbin[7]) return 7;
289 if (fabs(et) <= m_pbin[9] && fabs(et) >= m_pbin[8]) return 8;
290 return 0;
291 }
292
293 }
294
| Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
|
This page was automatically generated by the
LXR engine.
|
|