129 {
130 StatusCode sc = StatusCode::SUCCESS;
131
132 MsgStream log(
msgSvc(), name());
133 log << MSG::INFO << "in execute()" << endreq;
134
135
136
137
138 setFilterPassed(false);
139
140 m_pass[0] += 1;
141
142 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
145
146 Vint iGood, ikm, ikp, ipip, ipim, iGam;
147 iGood.clear();
148 ikm.clear();
149 ikp.clear();
150 ipip.clear();
151 ipim.clear();
152 iGam.clear();
153
154 Vp4 ppip, ppim, pkm, pkp;
155 ppip.clear();
156 ppim.clear();
157 pkm.clear();
158 pkp.clear();
159
160 int TotCharge = 0;
161 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
163 if(!(*itTrk)->isMdcTrackValid()) continue;
165 if(fabs(mdcTrk->
z()) >= m_vz0cut)
continue;
166 if(mdcTrk->
r() >= m_vr0cut)
continue;
167 iGood.push_back(i);
168 TotCharge += mdcTrk->
charge();
169 }
170
171
172
173 int nGood = iGood.size();
174
175
176
177
178
179 if((nGood < 2) || (TotCharge!=0)) return sc;
180
181 m_pass[1] += 1;
182
183
184
185
187 for(int i = 0; i < nGood; i++) {
189
196
197
200
201 if(!(*itTrk)->isMdcKalTrackValid()) continue ;
205 HepLorentzVector ptrk;
206 ptrk.setPx(mdcKalTrk->
px());
207 ptrk.setPy(mdcKalTrk->
py());
208 ptrk.setPz(mdcKalTrk->
pz());
209 double p3 = ptrk.mag();
210 ptrk.setE(sqrt(p3*p3+
mk*
mk));
211 if(mdcKalTrk->
charge() >0) {
212 ikp.push_back(iGood[i]);
213 pkp.push_back(ptrk);
214 } else {
215 ikm.push_back(iGood[i]);
216 pkm.push_back(ptrk);
217 }
218 }
219
222 HepLorentzVector ptrk;
223 ptrk.setPx(mdcKalTrk->
px());
224 ptrk.setPy(mdcKalTrk->
py());
225 ptrk.setPz(mdcKalTrk->
pz());
226 double p3 = ptrk.mag();
227 ptrk.setE(sqrt(p3*p3+
mpi*
mpi));
228 if(mdcKalTrk->
charge() >0) {
229 ipip.push_back(iGood[i]);
230 ppip.push_back(ptrk);
231 } else {
232 ipim.push_back(iGood[i]);
233 ppim.push_back(ptrk);
234 }
235 }
236 }
237 int npip = ipip.size();
238 int npim = ipim.size();
239 int nkm = ikm.size();
240 int nkp = ikp.size();
241
242 m_nkm=nkm;
243 m_nkp=nkp;
244 m_npip=npip;
245 m_npim=npim;
246 m_ncharge=nGood;
247 m_pass[3] += 1;
248
249 if(npip < 1 && npim < 1) return sc;
250 if(nkp < 1 && nkm < 1) return sc;
251
252 m_pass[4] += 1;
253
254
255
256
257
258 HepLorentzVector pkstar, pkstar1, ppi1, ppi2, pTot;
260 ikstar.clear();
261
262 double difchi0=99999.0;
263 int ixpim = -1;
264 int ixkp = -1;
265
266 for (int i = 0; i < npim; i++) {
267 for (int j = 0; j < nkp; j++) {
268 pkstar = ppim[i] + pkp[j];
269 double difchi =fabs(pkstar.m()-
mkstar0);
270 if(difchi < difchi0){
271 difchi0 = difchi;
272 ixpim = i;
273 ixkp = j;
274 }
275 }
276 }
277
278 m_difchikp = difchi0;
279
280 if(ixpim != -1) m_kstarkp=(ppim[ixpim] + pkp[ixkp]).m();
281
282 double difchi1=99999.0;
283 int ixpip = -1;
284 int ixkm = -1;
285
286 for (int ii = 0; ii < npip; ii++) {
287 for (int jj = 0; jj < nkm; jj++) {
288 pkstar1 = ppip[ii] + pkm[jj];
289 double difchi2 =fabs(pkstar1.m()-
mkstar0);
290 if(difchi2 < difchi1){
291 difchi1 = difchi2;
292 ixpip = ii;
293 ixkm = jj;
294 }
295 }
296 }
297
298 m_difchikm = difchi1;
299
300 if(ixpip != -1) m_kstarkm=(ppip[ixpip] + pkm[ixkm]).m();
301
302
303 if(ixpip == -1 && ixpim == -1) return sc;
304
305 if(difchi0 < difchi1){
306 pTot = ppim[ixpim] + pkp[ixkp];
307
308 } else{
309 pTot = ppip[ixpip] + pkm[ixkm];
310
311 }
312 m_mkstar = pTot.m();
313 m_pkstar = pTot.rho();
314
315
316 TH1 *h(0);
317 if (m_thsvc->getHist("/DQAHist/InclKstar/InclKstar_mass", h).isSuccess()) {
318 h->Fill(m_mkstar);
319 } else {
320 log << MSG::ERROR << "Couldn't retrieve InclKstar_mass" << endreq;
321 }
322
323 m_tuple2->write();
324
325
326
327 if(ixpim != -1){
328
329
330 (*(evtRecTrkCol->begin()+ipim[ixpim]))->tagPion();
331 (*(evtRecTrkCol->begin()+ikp[ixkp]))->tagKaon();
332 (*(evtRecTrkCol->begin()+ipim[ixpim]))->setQuality(3);
333 (*(evtRecTrkCol->begin()+ikp[ixkp]))->setQuality(3);
334 } else {
335
336
337 (*(evtRecTrkCol->begin()+ipip[ixpip]))->tagPion();
338 (*(evtRecTrkCol->begin()+ikm[ixkm]))->tagKaon();
339 (*(evtRecTrkCol->begin()+ipip[ixpip]))->setQuality(3);
340 (*(evtRecTrkCol->begin()+ikm[ixkm]))->setQuality(3);
341 }
342
343
344
345 setFilterPassed(true);
346
347
348 return StatusCode::SUCCESS;
349}
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
static void setPidType(PidType pidType)
int methodProbability() const
void setChiMinCut(const double chi=4)
void setRecTrack(EvtRecTrack *trk)
void setMethod(const int method)
void identify(const int pidcase)
void usePidSys(const int pidsys)
static ParticleID * instance()
bool IsPidInfoValid() const
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol