99 {
100
101 MsgStream log(
msgSvc(), name());
102 log << MSG::INFO << "in execute()" << endreq;
103
104
105
106
107 setFilterPassed(false);
108
109 SmartDataPtr<Event::EventHeader> eventHeader(eventSvc(),"/Event/EventHeader");
110 m_runNo=eventHeader->runNumber();
111 m_event=eventHeader->eventNumber();
112 log << MSG::DEBUG <<"run, evtnum = "
113 << m_runNo << " , "
114 << m_event <<endreq;
115
116
118 log << MSG::DEBUG <<"ncharg, nneu, tottks = "
119 << evtRecEvent->totalCharged() << " , "
120 << evtRecEvent->totalNeutral() << " , "
121 << evtRecEvent->totalTracks() <<endreq;
122
124
125 if(evtRecEvent->totalNeutral()>100) {
126 return StatusCode::SUCCESS;
127 }
128
129 Vint iGood, ipip, ipim;
130 iGood.clear();
131 ipip.clear();
132 ipim.clear();
134 ppip.clear();
135 ppim.clear();
136
137 Hep3Vector xorigin(0,0,0);
138
140 Gaudi::svcLocator()->service("VertexDbSvc", vtxsvc);
144 xorigin.setX(dbv[0]);
145 xorigin.setY(dbv[1]);
146 xorigin.setZ(dbv[2]);
147 }
148
149 int nCharge = 0;
150 for(int i = 0; i < evtRecEvent->totalCharged(); i++){
152 if(!(*itTrk)->isMdcTrackValid()) continue;
153 if (!(*itTrk)->isMdcKalTrackValid()) continue;
155
156 double pch =mdcTrk->
p();
157 double x0 =mdcTrk->
x();
158 double y0 =mdcTrk->
y();
159 double z0 =mdcTrk->
z();
160 double phi0=mdcTrk->
helix(1);
161 double xv=xorigin.x();
162 double yv=xorigin.y();
163 double Rxy=fabs((x0-xv)*
cos(phi0)+(y0-yv)*
sin(phi0));
164
165 if(fabs(z0) >= m_vz0cut) continue;
166 if(Rxy >= m_vr0cut) continue;
167
168 iGood.push_back((*itTrk)->trackId());
169 nCharge += mdcTrk->
charge();
170 if (mdcTrk->
charge() > 0) {
171 ipip.push_back((*itTrk)->trackId());
172 } else {
173 ipim.push_back((*itTrk)->trackId());
174 }
175 }
176
177
178
179
180 int nGood = iGood.size();
181 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endreq;
182 if((nGood != 2)||(nCharge!=0)){
183 return StatusCode::SUCCESS;
184 }
185
187 iGam.clear();
188 for(int i = evtRecEvent->totalCharged(); i< evtRecEvent->totalTracks(); i++) {
190 if(!(*itTrk)->isEmcShowerValid()) continue;
192 Hep3Vector emcpos(emcTrk->
x(), emcTrk->
y(), emcTrk->
z());
193
194 double dthe = 200.;
195 double dphi = 200.;
196 double dang = 200.;
197 for(int j = 0; j < evtRecEvent->totalCharged(); j++) {
199 if(!(*jtTrk)->isExtTrackValid()) continue;
203
204 double angd = extpos.angle(emcpos);
205 double thed = extpos.theta() - emcpos.theta();
206 double phid = extpos.deltaPhi(emcpos);
207 thed = fmod(thed+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
208 phid = fmod(phid+CLHEP::twopi+CLHEP::twopi+
pi, CLHEP::twopi) - CLHEP::pi;
209
210
211
212 if(angd < dang) {
213 dang = angd;
214 dthe = thed;
215 dphi = phid;
216 }
217 }
218 if(dang>=200) continue;
219 double eraw = emcTrk->
energy();
220 dthe = dthe * 180 / (CLHEP::pi);
221 dphi = dphi * 180 / (CLHEP::pi);
222 dang = dang * 180 / (CLHEP::pi);
223 if(eraw < m_energyThreshold) continue;
224 if(dang < m_gammaAngCut) continue;
225
226
227
228
229 iGam.push_back((*itTrk)->trackId());
230 }
231
232
233
234
235 int nGam = iGam.size();
236
237 log << MSG::DEBUG << "num Good Photon " << nGam << " , " <<evtRecEvent->totalNeutral()<<endreq;
238 if(nGam<2){
239 return StatusCode::SUCCESS;
240 }
241
242
243
244
245
246
247
248 m_tuple->write();
249
250
251
252
253
254
255
256 (*(evtRecTrkCol->begin()+ipip[0]))->setPartId(2);
257 (*(evtRecTrkCol->begin()+ipim[0]))->setPartId(2);
258
259
260
261
262
263 (*(evtRecTrkCol->begin()+ipip[0]))->setQuality(0);
264 (*(evtRecTrkCol->begin()+ipim[0]))->setQuality(0);
265
266
267
268
269 setFilterPassed(true);
270
271
272 return StatusCode::SUCCESS;
273
274}
double sin(const BesAngle a)
double cos(const BesAngle a)
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
const Hep3Vector emcPosition() const
const int emcVolumeNumber() const
const HepVector helix() const
......
virtual bool isVertexValid()=0
virtual double * SigmaPrimaryVertex()=0
virtual double * PrimaryVertex()=0
_EXTERN_ std::string EvtRecEvent
_EXTERN_ std::string EvtRecTrackCol