BOSS 7.0.7
BESIII Offline Software System
Loading...
Searching...
No Matches
CosmicGenerator Class Reference

#include <CosmicGenerator.h>

+ Inheritance diagram for CosmicGenerator:

Public Member Functions

 CosmicGenerator (const std::string &name, ISvcLocator *pSvcLocator)
 
 ~CosmicGenerator ()
 
StatusCode initialize ()
 
StatusCode execute ()
 
StatusCode finalize ()
 
HepLorentzVector generateVertex (void)
 

Static Public Attributes

static IBesRndmGenSvcp_BesRndmGenSvc = 0
 

Detailed Description

Definition at line 51 of file CosmicGenerator.h.

Constructor & Destructor Documentation

◆ CosmicGenerator()

CosmicGenerator::CosmicGenerator ( const std::string &  name,
ISvcLocator *  pSvcLocator 
)

Definition at line 153 of file CosmicGenerator.cxx.

154 : Algorithm(name,pSvcLocator)
155//--------------------------------------------------------------------------
156{
157 //
158 // Migration to MeV and mm units: all conversions are done in this interface
159 // to the CosmicGun. The CosmicGun itself uses GeV units internally - to call
160 // the fortran code.
161 //
162
163 m_GeV = 1000;
164 m_mm = 10;
165 m_readfile = false;
166 m_exzCut=false ;
167 m_tuple=0;
168 m_tuple1=0;
169 declareProperty("eventfile", m_infile = "NONE" );
170 declareProperty("emin", m_emin =0.1 );
171 declareProperty("emax", m_emax =10. );
172 declareProperty("xvert_low", m_xlow =-110.7); //EMC BSCRmax2
173 declareProperty("xvert_hig", m_xhig =110.7 );
174 declareProperty("zvert_low", m_zlow =-171.2);//EMC WorldZposition+0.5*CrystalLength
175 declareProperty("zvert_hig", m_zhig = 171.2 );
176 declareProperty("yvert_val", m_yval = 0 );
177 declareProperty("tmin", m_tmin =0. );
178 declareProperty("tmax", m_tmax =24. );
179 declareProperty("tcor", m_tcor =0. );
180 declareProperty("IPx", m_IPx =0. );
181 declareProperty("IPy", m_IPy =0. );
182 declareProperty("IPz", m_IPz =0. );
183 declareProperty("Radius", m_radius =0. );
184 declareProperty("ExzCut", m_exzCut );
185 declareProperty("CubeProjection", m_cubeProj = true );
186 declareProperty("OptimizeSphere", m_sphereOpt = false );
187
188
189 declareProperty("SwapYZAxis", m_swapYZAxis = false);
190 declareProperty("ctcut", m_ctcut =0.0 );// according to sin(acos(0.93)) =0.368
191 // declareProperty("ReadTTR", m_readTTR =false );
192 // declareProperty("ReadTR", m_readTTR =false );
193 declareProperty("PrintEvent", m_printEvent=10);
194 declareProperty("PrintMod", m_printMod=100);
195 declareProperty("RMax", m_rmax = 10000000. );
196 declareProperty("ThetaMin", m_thetamin = 0.);
197 declareProperty("ThetaMax", m_thetamax = 1.);
198 declareProperty("PhiMin", m_phimin = -1*PI);
199 declareProperty("PhiMax", m_phimax = PI);
200 declareProperty("Xposition", m_xpos = 263.5-0.0001); //263.5*cm,263.5*cm,287.5*cm
201 declareProperty("Yposition", m_ypos = 263.5-0.0001);
202 declareProperty("Zposition", m_zpos = 287.5-0.0001);
203
204 declareProperty("DumpMC", m_dumpMC = 0);
205}
const double PI
Definition: PipiJpsi.cxx:55

◆ ~CosmicGenerator()

CosmicGenerator::~CosmicGenerator ( )

Definition at line 208 of file CosmicGenerator.cxx.

210{}

Member Function Documentation

◆ execute()

StatusCode CosmicGenerator::execute ( )

Definition at line 363 of file CosmicGenerator.cxx.

363 {
364//---------------------------------------------------------------------------
365 MsgStream log(messageService(), name());
366 log << MSG::INFO << "CosmicGenerator executing" << endreq;
367
368
369 ++m_events;
370 // MsgStream log(messageService(), name());
371 log << MSG::DEBUG << "Event #" << m_events << endreq;
372
373
374 // clear up the vectors
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 // units are already converted to MeV's and mm.
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++;
424 vert = generateVertex();
425 Hep3Vector vert3(vert.x(),vert.y(),vert.z());
426
427 pp = gun->GenerateEvent();
428 if(m_dumpMC==1) {
429 m_cosmicE=pp.e()*m_GeV;
430 m_cosmicTheta=pp.theta();
431 m_cosmicPhi=pp.phi();
432 m_cosmicCharge=gun->GetMuonCharge();
433 m_tuple->write();
434 }
435 double theta1=pp.theta();
436 double phi1=pp.phi();
437 double mag1=pp.vect().mag();
438
439 Hep3Vector pp_corr(mag1*sin(theta1)*cos(phi1),
440 -mag1*cos(theta1),
441 mag1*sin(theta1)*sin(phi1));
442 Hep3Vector direction(pp_corr.x(),pp_corr.y(), pp_corr.z());
443 Hep3Vector center_dir=m_center-vert3;
444 // if optimization activated, check for the direction of the generated muon
445
446
447 if (m_cubeProj) {
448 double alpha,beta;
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;
467 HepPoint3D vert_pro0;
468 HepPoint3D vert_pro1;
469 HepPoint3D vert_pro2;
470 HepPoint3D vert_org(vert3);
471
472 coor_y0 = m_ypos; // defined in jobOpt.
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; // defined in jobOpt.
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; // defined in jobOpt.
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 //reset vert(x,y,z,t), after projection
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; // if no opt required accept the first muon
534 }
535
536 // pp.setX(pp.x()*m_GeV);
537 // pp.setY(pp.y()*m_GeV);
538 // pp.setZ(pp.z()*m_GeV);
539 // pp.setT(pp.t()*m_GeV);
540
541 pp.setX(pp.x());
542 pp.setY(pp.y());
543 pp.setZ(pp.z());
544 pp.setT(pp.t());
545 // Get the mass of the particle to be generated
546 int charge = gun->GetMuonCharge();
547 // m_pdgCode.push_back(charge*13);
548 m_pdgCode.push_back(charge*-13);
549
550 // const ParticleData* pdtparticle = m_particleTable->particle(ParticleID(abs(m_pdgCode.back())));
551 // double mass = pdtparticle->mass().value();
552 // double mass =105.5658;
553
554 // Compute the kinematic values. First, the vertex 4-vector:
555
556
557 // Set the polarization. Realistically, this is going to be zero
558 // for most studies, but you never know...
559 double polx = 0;
560 double poly = 0;
561 double polz = 0;
562 //m_polarization.set_normal3d(HepNormal3D(polx,poly,polz));
563 Polarization thePolarization;
564
565 // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ?
566 // if not...do nothing...if so, invert position of y- and z- coordinate
567 //
568 // well and don't forget about the direction of the incoming cosmic muon(s) either
569 // that means: y -> -y
570 //
571 if(!m_swapYZAxis){
572 // thePolarization.set_normal3d(HepNormal3D(polx,-poly,polz));
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 // The method of calculating e, theta, and phi depends on the user's
582 // commands. Let the KinematicManager handle it.
583 double e = pp.e();
584 double theta = pp.theta();
585 double phi = pp.phi();
586
587 // At this point, we have e, theta, and phi. Put them together to
588 // get the four-momentum.
589
590 double p2 = e*e - mass*mass;
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 // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ?
609 // if not...do nothing...if so, invert position of y- and z- coordinate
610 //
611 // well and don't forget about the direction of the incoming cosmic muon(s) either
612 // that means: y -> -y
613 //
614 if(!m_swapYZAxis) {
615 // Line below corrupted py sign and forces muons to be upwards, not downwards.
616 // m_fourMom.push_back(HepLorentzVector(px,-py,pz,pp.e()));
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()));
621 double x = vert.x();
622 double y = vert.y();
623 double z = vert.z();
624 double t = vert.t();
625 if(projected){
626
627 x = vert_proj.x();
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 // log << MSG::DEBUG
637 // << "proj:"<<m_cubeProj<<" ("<<x<<", "<<y<<", "<<z<<", "<<t<<")"
638 // <<endreq;
639 // Do we need to swap Y- and Z-axis for the PixelEndCap A Cosmic test ?
640 // if not...do nothing...if so, invert position of y- and z- coordinate
641 //
642 // but not only that...change also the direction of the incoming cosmic muon(s),
643 // they must go towards the pixel endcap A, i.e. y -> -y
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// m_cosmicE=pp.e()*m_GeV;
672// m_cosmicTheta=pp.theta();
673// m_cosmicPhi=pp.phi();
674// m_cosmicCharge=gun->GetMuonCharge();
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 // m_cosmicE=pz;
681 m_tuple1->write();
682 }
683 }
684
685// fill evt
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 // Note: The vertex and particle are owned by the event, so the
694 // event is responsible for those pointers.
695
696 // Create the particle, and specify its polarization.
697
698 GenParticle* particle = new GenParticle( m_fourMom[v], m_pdgCode[v], 1);
699 particle->set_polarization( m_polarization[v] );
700
701 // Create the vertex, and add the particle to the vertex.
702 GenVertex* vertex = new GenVertex(m_fourPos[v]);
703 vertex->add_particle_out( particle );
704
705 // Add the vertex to the event.
706 event->add_vertex( vertex );
707
708 }
709
710 event->set_event_number(m_events); // Set the event number
711
712
713
714 } else {
715
716 log<<MSG::ERROR<<"Wrong different number of vertexes/momenta/polaritazions!"<<endreq;
717 return StatusCode::FAILURE;
718 }
719 // Check if the McCollection already exists
720 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen");
721 if (anMcCol!=0)
722 {
723 // Add event to existing collection
724
725 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq;
726 McGenEvent* mcEvent = new McGenEvent(event);
727 anMcCol->push_back(mcEvent);
728 }
729 else
730 {
731 // Create Collection and add to the transient store
732 McGenEventCol *mcColl = new McGenEventCol;
733 McGenEvent* mcEvent = new McGenEvent(event);
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}
double sin(const BesAngle a)
Definition: BesAngle.h:210
double cos(const BesAngle a)
Definition: BesAngle.h:213
double mass
HepGeom::Point3D< double > HepPoint3D
HepGeom::Normal3D< float > HepNormal3D
const double alpha
**********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
Definition: KarLud.h:35
ObjectVector< McGenEvent > McGenEventCol
Definition: McGenEvent.h:39
@ theta1
Definition: TrkKalDeriv.h:24
TTree * t
Definition: binning.cxx:23
const HepLorentzVector & Vertex(void)
const HepLorentzVector & Momentum(void)
HepLorentzVector generateVertex(void)
int GetMuonCharge(void)
Definition: CosmicGun.cxx:134
HepLorentzVector GenerateEvent(void)
Definition: CosmicGun.cxx:109
static CosmicGun * GetCosmicGun(void)
Definition: CosmicGun.cxx:49
double y[1000]
double x[1000]
float charge

◆ finalize()

StatusCode CosmicGenerator::finalize ( )

Definition at line 757 of file CosmicGenerator.cxx.

757 {
758//---------------------------------------------------------------------------
759 // Get the KinematicManager.
760
761 if(m_cubeProj) {
762
763 MsgStream log(messageService(), name());
764
765 log << MSG::INFO<<"***************************************************"<<endreq;
766 log << MSG::INFO <<"** you have run CosmicGenerator with some " <<endreq;
767 log << MSG::INFO <<"** filters for cosmic muon simulation" <<endreq;
768 log << MSG::INFO <<"** "<<m_tried<<" muons were tried" <<endreq;
769 log << MSG::INFO <<"** "<<m_accepted<<" muons were accepted" <<endreq;
770 log << MSG::INFO <<"** "<<m_rejected<<" muons were rejected" <<endreq;
771 log << MSG::INFO<<"***************************************************"<<endreq;
772 std::cout<<"***************************************************"<<std::endl;
773 std::cout <<"** you have run CosmicGenerator with some " <<std::endl;
774 std::cout <<"** filters for cosmic muon simulation" <<std::endl;
775 std::cout <<"** "<<m_tried<<" muons were tried" <<std::endl;
776 std::cout <<"** "<<m_accepted<<" muons were accepted" <<std::endl;
777 std::cout <<"** "<<m_rejected<<" muons were rejected" <<std::endl;
778 std::cout<<"***************************************************"<<std::endl;
779
780 }
781
782 return StatusCode::SUCCESS;
783}

◆ generateVertex()

HepLorentzVector CosmicGenerator::generateVertex ( void  )

Definition at line 334 of file CosmicGenerator.cxx.

334 {
335
336 MsgStream log(messageService(), name());
337
338 // Get the pointer to the engine of the stream named SINGLE. If the
339 // stream does not exist is created automaticaly
340// HepRandomEngine* engine = CosmicGenerator::p_AtRndmGenSvc->GetEngine("COSMICS");
341 HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine("PYTHIA");
342 // Generate a random number according to the distribution.
343
344 float x_val = RandFlat::shoot(engine, m_xlow, m_xhig);
345 float z_val = RandFlat::shoot(engine, m_zlow, m_zhig);
346
347 // Generate a random number for time offset
348 float t_val=0;
349 if(m_tmin < m_tmax){
350 t_val = RandFlat::shoot(engine, m_tmin, m_tmax);
351 }
352 else if(m_tmin == m_tmax){
353 t_val = m_tmin;
354 }
355 else log << MSG::FATAL << " You specified m_tmin = " << m_tmin << " and m_tmax " << m_tmax << endreq;
356 HepLorentzVector p(x_val,m_yval,z_val, t_val*c_light);
357
358 return p;
359
360}
static IBesRndmGenSvc * p_BesRndmGenSvc
virtual CLHEP::HepRandomEngine * GetEngine(const std::string &StreamName)=0
Interface to the CLHEP engine.

Referenced by execute().

◆ initialize()

StatusCode CosmicGenerator::initialize ( )

Definition at line 213 of file CosmicGenerator.cxx.

213 {
214//---------------------------------------------------------------------------
215
216 // Initialize event count.
217
218 m_events = 0;
219 m_tried=0;
220 m_accepted=0;
221 m_rejected=0;
222 if(m_emin<0.1) {m_emin=0.1;std::cout<<"WARNING: set emin to 0.1 GeV"<<std::endl;}
223 m_emin =m_emin *m_GeV;
224 m_emax =m_emax *m_GeV;
225 m_xlow =m_xlow *m_mm;
226 m_xhig =m_xhig *m_mm;
227 m_zlow =m_zlow *m_mm;
228 m_zhig =m_zhig *m_mm;
229 m_yval =m_yval *m_mm;
230 m_xpos =m_xpos *m_mm;
231 m_ypos =m_ypos *m_mm;
232 m_zpos =m_zpos *m_mm;
233 m_radius=m_radius*m_mm;
234 m_tcor=m_tcor
235 +sqrt((m_xpos-m_xlow)*(m_xpos-m_xlow)+(m_zpos-m_zlow)*(m_zpos-m_zlow)+(m_ypos-m_yval)*(m_ypos-m_yval))
236 /(m_emin/sqrt(m_emin*m_emin+mass*mass*m_GeV*m_GeV));
237
238 ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap.h
239 p_msgSvc = 0;
240 StatusCode result = svcLocator->service("MessageSvc", p_msgSvc);
241
242 if ( !result.isSuccess() || p_msgSvc == 0 )
243 {
244 std::cerr << "VGenCommand(): could not initialize MessageSvc!" << std::endl;
245 throw GaudiException("Could not initalize MessageSvc","CosmicGenerator",StatusCode::FAILURE);
246 }
247
248 MsgStream log (p_msgSvc,"CosmicGenerator");
249 if(m_dumpMC==1) {
250 StatusCode status;
251 NTuplePtr nt(ntupleSvc(), "FILE1/comp");
252 if(nt) m_tuple = nt;
253 else {
254 m_tuple = ntupleSvc()->book ("FILE1/comp", CLID_ColumnWiseTuple, "ks N-Tuple example");
255 if ( m_tuple ) {
256 status = m_tuple->addItem ("cosmic_e", m_cosmicE);
257 status = m_tuple->addItem ("cosmic_theta", m_cosmicTheta);
258 status = m_tuple->addItem ("cosmic_phi", m_cosmicPhi);
259 status = m_tuple->addItem ("cosmic_charge", m_cosmicCharge);
260 }
261 else {
262 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg;
263 return StatusCode::FAILURE;
264 }
265 }
266 NTuplePtr nt1(ntupleSvc(), "FILE1/compp");
267 if(nt1) m_tuple1 = nt1;
268 else {
269 m_tuple1 = ntupleSvc()->book ("FILE1/compp", CLID_ColumnWiseTuple, "ks N-Tuple example");
270 if ( m_tuple1 ) {
271 status = m_tuple1->addItem ("mc_theta", mc_theta);
272 status = m_tuple1->addItem ("mc_phi", mc_phi);
273 status = m_tuple1->addItem ("mc_px", mc_px);
274 status = m_tuple1->addItem ("mc_py", mc_py);
275 status = m_tuple1->addItem ("mc_pz", mc_pz);
276 }
277 else {
278 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg;
279 return StatusCode::FAILURE;
280 }
281 }
282 }
283 if(m_infile=="NONE")
284
285 {
286 // Initialize the BES random number services.
287
288 static const bool CREATEIFNOTTHERE(true);
289 StatusCode RndmStatus = svcLocator->service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE);
290
291 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc)
292 {
293 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
294 return RndmStatus;
295 }
296
298
299 gun->SetEnergyRange(m_emin/m_GeV,m_emax/m_GeV);
300 gun->SetCosCut(m_ctcut);
301 gun->PrintLevel(m_printEvent, m_printMod);
302 float flux_withCT = gun->InitializeGenerator();
303
304 log << MSG::INFO << "Initialisation cosmic gun done." << endreq;
305 log << MSG::INFO << "Accepted diff flux after E and cos(theta) cuts = " <<
306 flux_withCT << " /cm^2/s" << endreq;
307 log << MSG::INFO << "Accepted total flux after E and cos(theta) cuts = " <<
308 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << endreq;
309
310 std::cout<< "Accepted diff flux after E and cos(theta) cuts = " <<
311 flux_withCT << " /cm^2/s" << std::endl;
312 std::cout << "Accepted total flux after E and cos(theta) cuts = " <<
313 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << std::endl;
314
315 }
316 else
317 {
318 log << MSG::INFO << "Cosmics are read from file " << m_infile << endreq;
319 m_ffile.open(m_infile.c_str());
320 if(!m_ffile)
321 {
322 log << MSG::FATAL << "Could not open input file - stop! " << endreq;
323 return StatusCode::FAILURE;
324 }
325 m_readfile = true;
326 }
327
328 m_center=Hep3Vector(m_IPx, m_IPy, m_IPz);
329 log << MSG::INFO << "Initialisation done" << endreq;
330 return StatusCode::SUCCESS;
331
332}
INTupleSvc * ntupleSvc()
void PrintLevel(int printevt, int printmod)
Definition: CosmicGun.cxx:90
void SetEnergyRange(float emin, float emax)
Definition: CosmicGun.cxx:138
float InitializeGenerator()
Definition: CosmicGun.cxx:81
void SetCosCut(float ctcut)
Definition: CosmicGun.cxx:154

Member Data Documentation

◆ p_BesRndmGenSvc

IBesRndmGenSvc * CosmicGenerator::p_BesRndmGenSvc = 0
static

Definition at line 62 of file CosmicGenerator.h.

Referenced by cosmicrndm_(), generateVertex(), and initialize().


The documentation for this class was generated from the following files: