199 {
200 MsgStream log(
msgSvc(), name());
201 log << MSG::INFO << "Ekhara in execute()" << endreq;
202
204
205
206
207
208
209 double p1[4],p2[4],q1[4],q2[4];
210 double pi1[4],pi2[4],qpion[4];
211
212 for (int i=0;i<4;i++) {
213 p1[i]=0.0;
214 p2[i]=0.0;
215 q1[i]=0.0;
216 q2[i]=0.0;
217 pi1[i]=0.0;
218 pi2[i]=0.0;
219 qpion[i]=0.0;
220 }
221
223
224 double PosiMom = sqrt(q1[1]*q1[1]+q1[2]*q1[2]+q1[3]*q1[3]);
225 double PosiThe = acos(q1[3]/PosiMom);
226 double PosiPhi = atan2(q1[2],q1[1]);
227
228 double ElecMom = sqrt(q2[1]*q2[1]+q2[2]*q2[2]+q2[3]*q2[3]);
229 double ElecThe = acos(q2[3]/ElecMom);
230 double ElecPhi = atan2(q2[2],q2[1]);
231
232 double EtaPMom = sqrt(qpion[1]*qpion[1]+qpion[2]*qpion[2]+qpion[3]*qpion[3]);
233 double EtaPThe = acos(qpion[3]/EtaPMom);
234 double EtaPPhi = atan2(qpion[2],qpion[1]);
235
236 hMCPosiMom->fill(PosiMom);
237 hMCElecMom->fill(ElecMom);
238 hMCEtaPMom->fill(EtaPMom);
239
240 hMCPosiPhi->fill(PosiPhi*TMath::RadToDeg());
241 hMCElecPhi->fill(ElecPhi*TMath::RadToDeg());
242 hMCEtaPPhi->fill(EtaPPhi*TMath::RadToDeg());
243
244 hMCPosiThe->fill(PosiThe*TMath::RadToDeg());
245 hMCElecThe->fill(ElecThe*TMath::RadToDeg());
246 hMCEtaPThe->fill(EtaPThe*TMath::RadToDeg());
247
248
249
250
251 GenEvent* evt = new GenEvent(1,1);
252
253 GenVertex* prod_vtx = new GenVertex();
254 evt->add_vertex( prod_vtx );
255
256
257 GenParticle* p = new GenParticle( HepLorentzVector(p1[1],p1[2],p1[3],p1[0]), -11, 3);
258 p->suggest_barcode(1 );
259 prod_vtx->add_particle_in(p);
260
261
262 p = new GenParticle(HepLorentzVector(p2[1],p2[2],p2[3],p2[0]), 11, 3);
263 p->suggest_barcode( 2 );
264 prod_vtx->add_particle_in(p);
265
266
267
268 p = new GenParticle(HepLorentzVector(q1[1],q1[2],q1[3],q1[0]), -11, 1);
269 p->suggest_barcode(3 );
270 prod_vtx->add_particle_out(p);
271
272
273 p = new GenParticle( HepLorentzVector(q2[1],q2[2],q2[3],q2[0] ), 11, 1);
274 p->suggest_barcode( 4 );
275 prod_vtx->add_particle_out(p);
276
277 if (m_channel_id == 2){
278 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 111, 1);
279 p->suggest_barcode( 5 );
280 prod_vtx->add_particle_out(p);
281 }
282 else if (m_channel_id == 3){
283 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 221, 1);
284 p->suggest_barcode( 5 );
285 prod_vtx->add_particle_out(p);
286 }
287 else if (m_channel_id == 4){
288 p = new GenParticle( HepLorentzVector( qpion[1],qpion[2],qpion[3],qpion[0]), 331, 1);
289 p->suggest_barcode( 5 );
290 prod_vtx->add_particle_out(p);
291 }
292 else if (m_channel_id == 1) {
293
294 p = new GenParticle( HepLorentzVector( pi1[1],pi1[2],pi1[3],pi1[0]),211, 1);
295 p->suggest_barcode( 5 );
296 prod_vtx->add_particle_out(p);
297
298 p = new GenParticle( HepLorentzVector(pi2[1],pi2[2],pi2[3],pi2[0]), -211, 1);
299 p->suggest_barcode( 6 );
300 prod_vtx->add_particle_out(p);
301 }
302
303 if( log.level() <= MSG::INFO )
304 {
305
306 evt->print();
307
308
309
310 }
311
312
313 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
314 if (anMcCol!=0)
315 {
316
317 MsgStream log(messageService(), name());
318 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
320 anMcCol->push_back(mcEvent);
321 }
322 else
323 {
324
327 mcColl->push_back(mcEvent);
328 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
329 if (sc != StatusCode::SUCCESS)
330 {
331 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
332 delete mcColl;
333 delete evt;
334 delete mcEvent;
335 return StatusCode::FAILURE;
336 }
337 else
338 {
339 log << MSG::INFO << "McGenEventCol created " << endreq;
340 }
341
342 }
343
344 log<<MSG::INFO<< "before execute() return"<<endreq;
345 return StatusCode::SUCCESS;
346
347}
#define GET_MOMENTUM(p1, p2, q1, q2, pi1, pi2, qpion)
ObjectVector< McGenEvent > McGenEventCol