86 {
87
88 MsgStream log(
msgSvc(), name() );
89 log << MSG::INFO << "in execute()" << endmsg;
90
91
92
93
94 setFilterPassed( false );
95
96 SmartDataPtr<Event::EventHeader> eventHeader( eventSvc(), "/Event/EventHeader" );
97 m_runNo = eventHeader->runNumber();
98 m_event = eventHeader->eventNumber();
99 log << MSG::DEBUG << "run, evtnum = " << m_runNo << " , " << m_event << endmsg;
100
102 log << MSG::DEBUG << "ncharg, nneu, tottks = " << evtRecEvent->totalCharged() << " , "
103 << evtRecEvent->totalNeutral() << " , " << evtRecEvent->totalTracks() << endmsg;
104
106
107 if ( evtRecEvent->totalNeutral() > 100 ) { return StatusCode::SUCCESS; }
108
109 Vint iGood, ipip, ipim;
110 iGood.clear();
111 ipip.clear();
112 ipim.clear();
114 ppip.clear();
115 ppim.clear();
116
117 Hep3Vector xorigin( 0, 0, 0 );
118
119 IVertexDbSvc* vtxsvc;
120 Gaudi::svcLocator()->service( "VertexDbSvc", vtxsvc ).ignore();
122 {
125 xorigin.setX( dbv[0] );
126 xorigin.setY( dbv[1] );
127 xorigin.setZ( dbv[2] );
128 }
129
130 int nCharge = 0;
131 for ( int i = 0; i < evtRecEvent->totalCharged(); i++ )
132 {
134 if ( !( *itTrk )->isMdcTrackValid() ) continue;
135 if ( !( *itTrk )->isMdcKalTrackValid() ) continue;
136 RecMdcTrack* mdcTrk = ( *itTrk )->mdcTrack();
137
138 double pch = mdcTrk->
p();
139 double x0 = mdcTrk->
x();
140 double y0 = mdcTrk->
y();
141 double z0 = mdcTrk->
z();
142 double phi0 = mdcTrk->
helix( 1 );
143 double xv = xorigin.x();
144 double yv = xorigin.y();
145 double Rxy = fabs( ( x0 - xv ) *
cos( phi0 ) + ( y0 - yv ) *
sin( phi0 ) );
146
147 if ( fabs( z0 ) >= m_vz0cut ) continue;
148 if ( Rxy >= m_vr0cut ) continue;
149
150 iGood.push_back( ( *itTrk )->trackId() );
151 nCharge += mdcTrk->
charge();
152 if ( mdcTrk->
charge() > 0 ) { ipip.push_back( ( *itTrk )->trackId() ); }
153 else { ipim.push_back( ( *itTrk )->trackId() ); }
154 }
155
156
157
158
159 int nGood = iGood.size();
160 log << MSG::DEBUG << "ngood, totcharge = " << nGood << " , " << nCharge << endmsg;
161 if ( ( nGood != 2 ) || ( nCharge != 0 ) ) { return StatusCode::SUCCESS; }
162
164 iGam.clear();
165 for ( int i = evtRecEvent->totalCharged(); i < evtRecEvent->totalTracks(); i++ )
166 {
168 if ( !( *itTrk )->isEmcShowerValid() ) continue;
169 RecEmcShower* emcTrk = ( *itTrk )->emcShower();
170 Hep3Vector emcpos( emcTrk->
x(), emcTrk->
y(), emcTrk->
z() );
171
172 double dthe = 200.;
173 double dphi = 200.;
174 double dang = 200.;
175 for ( int j = 0; j < evtRecEvent->totalCharged(); j++ )
176 {
178 if ( !( *jtTrk )->isExtTrackValid() ) continue;
179 RecExtTrack* extTrk = ( *jtTrk )->extTrack();
182
183 double angd = extpos.angle( emcpos );
184 double thed = extpos.theta() - emcpos.theta();
185 double phid = extpos.deltaPhi( emcpos );
186 thed = fmod( thed + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
187 phid = fmod( phid + CLHEP::twopi + CLHEP::twopi +
pi, CLHEP::twopi ) - CLHEP::pi;
188
189
190
191 if ( angd < dang )
192 {
193 dang = angd;
194 dthe = thed;
195 dphi = phid;
196 }
197 }
198 if ( dang >= 200 ) continue;
199 double eraw = emcTrk->
energy();
200 dthe = dthe * 180 / ( CLHEP::pi );
201 dphi = dphi * 180 / ( CLHEP::pi );
202 dang = dang * 180 / ( CLHEP::pi );
203 if ( eraw < m_energyThreshold ) continue;
204 if ( dang < m_gammaAngCut ) continue;
205
206
207
208
209 iGam.push_back( ( *itTrk )->trackId() );
210 }
211
212
213
214
215 int nGam = iGam.size();
216
217 log << MSG::DEBUG << "num Good Photon " << nGam << " , " << evtRecEvent->totalNeutral()
218 << endmsg;
219 if ( nGam < 2 ) { return StatusCode::SUCCESS; }
220
221
222
223
224
225
226
227 m_tuple->write().ignore();
228
229
230
231
232
233
234
235 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setPartId( 2 );
236 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setPartId( 2 );
237
238
239
240
241
242 ( *( evtRecTrkCol->begin() + ipip[0] ) )->setQuality( 0 );
243 ( *( evtRecTrkCol->begin() + ipim[0] ) )->setQuality( 0 );
244
245
246
247
248 setFilterPassed( true );
249
250
251 return StatusCode::SUCCESS;
252}
EvtRecTrackCol::iterator EvtRecTrackIterator
std::vector< HepLorentzVector > Vp4
double sin(const BesAngle a)
double cos(const BesAngle a)
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