363 {
364
365 MsgStream log(messageService(), name());
366 log << MSG::INFO << "CosmicGenerator executing" << endreq;
367
368
369 ++m_events;
370
371 log << MSG::DEBUG << "Event #" << m_events << endreq;
372
373
374
375 m_fourPos.clear();
376 m_fourMom.clear();
377 m_polarization.clear();
378 m_pdgCode.clear();
379
380
381 if(m_readfile)
382 {
383 if(!m_ffile.eof())
384 {
386 m_ffile >> evt;
387
388 std::cout << evt << std::endl;
389
390 double polx = 0;
391 double poly = 0;
392 double polz = 0;
393 Polarization thePolarization;
394 thePolarization.set_normal3d(
HepNormal3D(polx,poly,polz));
395 m_polarization.push_back(thePolarization);
396
397
398
399
400 m_fourPos.push_back(evt.
Vertex());
401 m_fourMom.push_back(evt.
Momentum());
402 m_pdgCode.push_back(evt.
pdgID());
403
404 }
405 else
406 {
407 log << MSG::FATAL << "End of file reached - stop " << endreq;
408 exit(1);
409 return StatusCode::FAILURE;
410 }
411 }
412
413 else
414 {
415
416 bool accepted=false;
417 bool projected=false;
418 HepLorentzVector pp;
420 HepLorentzVector vert;
421 HepLorentzVector vert_proj;
422 while(!accepted){
423 m_tried++;
425 Hep3Vector vert3(vert.x(),vert.y(),vert.z());
426
428 if(m_dumpMC==1) {
429 m_cosmicE=pp.e()*m_GeV;
430 m_cosmicTheta=pp.theta();
431 m_cosmicPhi=pp.phi();
433 m_tuple->write();
434 }
436 double phi1=pp.phi();
437 double mag1=pp.vect().mag();
438
442 Hep3Vector direction(pp_corr.x(),pp_corr.y(), pp_corr.z());
443 Hep3Vector center_dir=m_center-vert3;
444
445
446
447 if (m_cubeProj) {
449 if(m_sphereOpt){
450
451 beta=direction.angle(center_dir);
452 alpha=asin(m_radius/center_dir.mag());
453 if(fabs(beta)>
alpha){
454 log<<MSG::DEBUG<<
alpha<<
", "<<beta<<
" rejected muon due to sphere cut! "<<endreq;
455 m_rejected++;
456 continue;
457 }
458 }
459
460 double l_fight0,l_fight1, l_fight2;
461 double coor_x0, coor_y0, coor_z0;
462 double coor_x1, coor_y1, coor_z1;
463 double coor_x2, coor_y2, coor_z2;
464 bool isIn0=false;
465 bool isIn1=false;
466 bool isIn2=false;
471
472 coor_y0 = m_ypos;
473 coor_x0 = direction.x()*(coor_y0 - vert.y())/direction.y() +vert.x();
474 coor_z0 = direction.z()*(coor_y0 - vert.y())/direction.y() +vert.z();
475
476
477
478 if( fabs(coor_x0) <= m_xpos && fabs(coor_z0) <= m_zpos){
479 vert_pro0=
HepPoint3D (coor_x0, coor_y0, coor_z0);
480 l_fight0=vert_org.distance(vert_pro0);
481 isIn0=true;
482 }
483 else{
484
485 coor_z1 = sign(coor_z0-vert.z())*m_zpos;
486 coor_x1 = direction.x()*(coor_z1 - vert.z())/direction.z() +vert.x();
487 coor_y1 = direction.y()*(coor_z1 - vert.z())/direction.z() +vert.y();
488
489 vert_pro1=
HepPoint3D (coor_x1, coor_y1, coor_z1);
490 l_fight1=vert_org.distance(vert_pro1);
491
492
493 coor_x2 = sign(coor_x0-vert.x())*m_xpos;
494 coor_z2 = direction.z()*(coor_x2 - vert.x())/direction.x() +vert.z();
495 coor_y2 = direction.y()*(coor_x2 - vert.x())/direction.x() +vert.y();
496 vert_pro2=
HepPoint3D (coor_x2, coor_y2, coor_z2);
497 l_fight2=vert_org.distance(vert_pro2);
498 if(l_fight1<l_fight2)
499 isIn1=true;
500 else
501 isIn2=true;
502 }
503
504 if(isIn0){
505 vert_proj=HepLorentzVector (vert_pro0.x(),vert_pro0.y(),vert_pro0.z() , l_fight0);
506
507 projected =true;
508 accepted = true;
509 m_accepted++;
510 }
511 else if(isIn1){
512 vert_proj=HepLorentzVector (vert_pro1.x(),vert_pro1.y(),vert_pro1.z() , l_fight1);
513
514 projected =true;
515 accepted = true;
516 m_accepted++;
517 }
518 else if(isIn2){
519 vert_proj=HepLorentzVector (vert_pro2.x(),vert_pro2.y(),vert_pro2.z() , l_fight2);
520
521 projected =true;
522 accepted = true;
523 m_accepted++;
524 } else {
525 log<<MSG::DEBUG<<" rejected muon due to cube projection request! "<<endreq;
526 m_rejected++;
527 }
528
529
530
531
532 }
533 else accepted=true;
534 }
535
536
537
538
539
540
541 pp.setX(pp.x());
542 pp.setY(pp.y());
543 pp.setZ(pp.z());
544 pp.setT(pp.t());
545
547
548 m_pdgCode.push_back(charge*-13);
549
550
551
552
553
554
555
556
557
558
559 double polx = 0;
560 double poly = 0;
561 double polz = 0;
562
563 Polarization thePolarization;
564
565
566
567
568
569
570
571 if(!m_swapYZAxis){
572
573 thePolarization.set_normal3d(
HepNormal3D(polx,poly,polz));
574 }
575 else
576 thePolarization.set_normal3d(
HepNormal3D(polx,polz,-poly));
577
578 m_polarization.push_back(thePolarization);
579
580
581
582
583 double e = pp.e();
584 double theta = pp.theta();
585 double phi = pp.phi();
586
587
588
589
591 if ( p2 < 0 )
592 {
593 log << MSG::ERROR
594 << "Event #" << m_events
595 <<
" E=" << e <<
", mass=" <<
mass
596 << " -- you have generated a tachyon! Increase energy or change particle ID."
597 << endmsg;
598 return StatusCode::FAILURE;
599 }
600
601 double p = sqrt(p2);
602 double px = p*
sin(theta)*
cos(phi);
603 double pz = p*
sin(theta)*
sin(phi);
604 double py = -p*
cos(theta);
605
606
607
608
609
610
611
612
613
614 if(!m_swapYZAxis) {
615
616
617 m_fourMom.push_back(HepLorentzVector(px,py,pz,pp.e()));
618 }
619 else
620 m_fourMom.push_back(HepLorentzVector(px,pz,-py,pp.e()));
622 double y = vert.y();
623 double z = vert.z();
625 if(projected){
626
628 y = vert_proj.y();
629 z = vert_proj.z();
630 t = vert.t()-vert_proj.t()/HepLorentzVector(px,py,pz,pp.e()).beta()
631 +m_tcor;
632
633
634
635 }
636
637
638
639
640
641
642
643
644
645
646
647 if(!m_swapYZAxis)
648 m_fourPos.push_back(HepLorentzVector(x,y,z,
t));
649 else
650 m_fourPos.push_back(HepLorentzVector(x,z,y,
t));
651
652 log << MSG::DEBUG
653 << " (x,y,z,t) = ("
654 << m_fourPos.back().x() << ","
655 << m_fourPos.back().y() << ","
656 << m_fourPos.back().z() << ","
657 << m_fourPos.back().t() << "), (Px,Py,Pz,E) = ("
658 << m_fourMom.back().px() << ","
659 << m_fourMom.back().py() << ","
660 << m_fourMom.back().pz() << ","
661 << m_fourMom.back().e() << ")"
662 << endreq;
663 log << MSG::DEBUG
664 << " (theta,phi) = (" << theta << "," << phi << "), "
665 << "polarization(x,y,z) = ("
666 << m_polarization.back().normal3d().x() << ","
667 << m_polarization.back().normal3d().y() << ","
668 << m_polarization.back().normal3d().z() << ")"
669 << endreq;
670 if(m_dumpMC==1) {
671
672
673
674
675 mc_theta=m_fourMom.back().theta();
676 mc_phi=m_fourMom.back().phi();
677 mc_px=m_fourMom.back().px();
678 mc_py=m_fourMom.back().py();
679 mc_pz=m_fourMom.back().pz();
680
681 m_tuple1->write();
682 }
683 }
684
685
686
687 GenEvent* event = new GenEvent(1,1);
688
689 if(m_fourMom.size()==m_fourPos.size()&&m_fourMom.size()==m_polarization.size()){
690
691 for(std::size_t
v=0;
v<m_fourMom.size();++
v){
692
693
694
695
696
697
698 GenParticle* particle =
new GenParticle( m_fourMom[
v], m_pdgCode[
v], 1);
699 particle->set_polarization( m_polarization[
v] );
700
701
702 GenVertex* vertex =
new GenVertex(m_fourPos[
v]);
703 vertex->add_particle_out( particle );
704
705
706 event->add_vertex( vertex );
707
708 }
709
710 event->set_event_number(m_events);
711
712
713
714 } else {
715
716 log<<MSG::ERROR<<"Wrong different number of vertexes/momenta/polaritazions!"<<endreq;
717 return StatusCode::FAILURE;
718 }
719
720 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
721 if (anMcCol!=0)
722 {
723
724
725 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
727 anMcCol->push_back(mcEvent);
728 }
729 else
730 {
731
734 mcColl->push_back(mcEvent);
735 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl);
736 if (sc != StatusCode::SUCCESS)
737 {
738 log << MSG::ERROR << "Could not register McGenEvent" << endreq;
739 delete mcColl;
740 delete event;
741 delete mcEvent;
742 return StatusCode::FAILURE;
743 }
744
745 }
746
747
748
749
750
751 return StatusCode::SUCCESS;
752
753}
ObjectVector< McGenEvent > McGenEventCol
HepGeom::Point3D< double > HepPoint3D
HepGeom::Normal3D< float > HepNormal3D
double sin(const BesAngle a)
double cos(const BesAngle a)
**********Class see also m_nmax DOUBLE PRECISION m_amel DOUBLE PRECISION m_x2 DOUBLE PRECISION m_alfinv DOUBLE PRECISION m_Xenph INTEGER m_KeyWtm INTEGER m_idyfs DOUBLE PRECISION m_zini DOUBLE PRECISION m_q2 DOUBLE PRECISION m_Wt_KF DOUBLE PRECISION m_WtCut INTEGER m_KFfin *COMMON c_KarLud $ !Input CMS energy[GeV] $ !CMS energy after beam spread beam strahlung[GeV] $ !Beam energy spread[GeV] $ !z boost due to beam spread $ !electron beam mass *ff pair spectrum $ !minimum v
const HepLorentzVector & Vertex(void)
const HepLorentzVector & Momentum(void)
HepLorentzVector generateVertex(void)
HepLorentzVector GenerateEvent(void)
static CosmicGun * GetCosmicGun(void)