#include <CosmicGenerator.h>
Public Member Functions | |
CosmicGenerator (const std::string &name, ISvcLocator *pSvcLocator) | |
~CosmicGenerator () | |
StatusCode | initialize () |
StatusCode | execute () |
StatusCode | finalize () |
HepLorentzVector | generateVertex (void) |
Static Public Attributes | |
static IBesRndmGenSvc * | p_BesRndmGenSvc = 0 |
Private Member Functions | |
bool | exzCut (const Hep3Vector &pos, const HepLorentzVector &p) |
Private Attributes | |
StoreGateSvc * | m_sgSvc |
ActiveStoreSvc * | m_activeStore |
NTuple::Tuple * | m_tuple |
NTuple::Item< double > | m_cosmicE |
NTuple::Item< double > | m_cosmicTheta |
NTuple::Item< double > | m_cosmicPhi |
NTuple::Item< double > | m_cosmicCharge |
NTuple::Tuple * | m_tuple1 |
NTuple::Item< double > | mc_phi |
NTuple::Item< double > | mc_theta |
NTuple::Item< double > | mc_px |
NTuple::Item< double > | mc_py |
NTuple::Item< double > | mc_pz |
long int | m_events |
long int | m_rejected |
long int | m_accepted |
long int | m_tried |
std::vector< int > | m_pdgCode |
float | m_emin |
float | m_emax |
float | m_ctcut |
float | m_xlow |
float | m_xhig |
float | m_zlow |
float | m_zhig |
float | m_yval |
float | m_IPx |
float | m_IPy |
float | m_IPz |
float | m_radius |
float | m_xpos |
float | m_ypos |
float | m_zpos |
float | m_tmin |
float | m_tmax |
float | m_tcor |
bool | m_sphereOpt |
bool | m_swapYZAxis |
bool | m_cubeProj |
long int | m_printEvent |
long int | m_printMod |
long int | m_selection |
int | m_dumpMC |
float | m_thetamin |
float | m_thetamax |
float | m_phimin |
float | m_phimax |
float | m_GeV |
float | m_mm |
bool | m_readfile |
std::string | m_infile |
std::ifstream | m_ffile |
std::vector< HepLorentzVector > | m_fourPos |
std::vector< HepLorentzVector > | m_fourMom |
Hep3Vector | m_center |
std::vector< HepMC::Polarization > | m_polarization |
IMessageSvc * | p_msgSvc |
bool | m_exzCut |
float | m_rmax |
Definition at line 51 of file CosmicGenerator.h.
CosmicGenerator::CosmicGenerator | ( | const std::string & | name, | |
ISvcLocator * | pSvcLocator | |||
) |
Definition at line 153 of file CosmicGenerator.cxx.
References false, m_ctcut, m_cubeProj, m_dumpMC, m_emax, m_emin, m_exzCut, m_GeV, m_infile, m_IPx, m_IPy, m_IPz, m_mm, m_phimax, m_phimin, m_printEvent, m_printMod, m_radius, m_readfile, m_rmax, m_sphereOpt, m_swapYZAxis, m_tcor, m_thetamax, m_thetamin, m_tmax, m_tmin, m_tuple, m_tuple1, m_xhig, m_xlow, m_xpos, m_ypos, m_yval, m_zhig, m_zlow, m_zpos, and PI.
00154 : Algorithm(name,pSvcLocator) 00155 //-------------------------------------------------------------------------- 00156 { 00157 // 00158 // Migration to MeV and mm units: all conversions are done in this interface 00159 // to the CosmicGun. The CosmicGun itself uses GeV units internally - to call 00160 // the fortran code. 00161 // 00162 00163 m_GeV = 1000; 00164 m_mm = 10; 00165 m_readfile = false; 00166 m_exzCut=false ; 00167 m_tuple=0; 00168 m_tuple1=0; 00169 declareProperty("eventfile", m_infile = "NONE" ); 00170 declareProperty("emin", m_emin =0.1 ); 00171 declareProperty("emax", m_emax =10. ); 00172 declareProperty("xvert_low", m_xlow =-110.7); //EMC BSCRmax2 00173 declareProperty("xvert_hig", m_xhig =110.7 ); 00174 declareProperty("zvert_low", m_zlow =-171.2);//EMC WorldZposition+0.5*CrystalLength 00175 declareProperty("zvert_hig", m_zhig = 171.2 ); 00176 declareProperty("yvert_val", m_yval = 0 ); 00177 declareProperty("tmin", m_tmin =0. ); 00178 declareProperty("tmax", m_tmax =24. ); 00179 declareProperty("tcor", m_tcor =0. ); 00180 declareProperty("IPx", m_IPx =0. ); 00181 declareProperty("IPy", m_IPy =0. ); 00182 declareProperty("IPz", m_IPz =0. ); 00183 declareProperty("Radius", m_radius =0. ); 00184 declareProperty("ExzCut", m_exzCut ); 00185 declareProperty("CubeProjection", m_cubeProj = true ); 00186 declareProperty("OptimizeSphere", m_sphereOpt = false ); 00187 00188 00189 declareProperty("SwapYZAxis", m_swapYZAxis = false); 00190 declareProperty("ctcut", m_ctcut =0.0 );// according to sin(acos(0.93)) =0.368 00191 // declareProperty("ReadTTR", m_readTTR =false ); 00192 // declareProperty("ReadTR", m_readTTR =false ); 00193 declareProperty("PrintEvent", m_printEvent=10); 00194 declareProperty("PrintMod", m_printMod=100); 00195 declareProperty("RMax", m_rmax = 10000000. ); 00196 declareProperty("ThetaMin", m_thetamin = 0.); 00197 declareProperty("ThetaMax", m_thetamax = 1.); 00198 declareProperty("PhiMin", m_phimin = -1*PI); 00199 declareProperty("PhiMax", m_phimax = PI); 00200 declareProperty("Xposition", m_xpos = 263.5-0.0001); //263.5*cm,263.5*cm,287.5*cm 00201 declareProperty("Yposition", m_ypos = 263.5-0.0001); 00202 declareProperty("Zposition", m_zpos = 287.5-0.0001); 00203 00204 declareProperty("DumpMC", m_dumpMC = 0); 00205 }
CosmicGenerator::~CosmicGenerator | ( | ) |
StatusCode CosmicGenerator::execute | ( | ) |
Definition at line 363 of file CosmicGenerator.cxx.
References alpha, cos(), Bes_Common::DEBUG, calibUtil::ERROR, Bes_Common::FATAL, CosmicGun::GenerateEvent(), generateVertex(), CosmicGun::GetCosmicGun(), CosmicGun::GetMuonCharge(), Bes_Common::INFO, m_accepted, m_center, m_cosmicCharge, m_cosmicE, m_cosmicPhi, m_cosmicTheta, m_cubeProj, m_dumpMC, m_events, m_ffile, m_fourMom, m_fourPos, m_GeV, m_pdgCode, m_polarization, m_radius, m_readfile, m_rejected, m_sphereOpt, m_swapYZAxis, m_tcor, m_tried, m_tuple, m_tuple1, m_xpos, m_ypos, m_zpos, mass, mc_phi, mc_px, mc_py, mc_pz, mc_theta, phi1, sign(), sin(), t(), theta1, v, and x.
00363 { 00364 //--------------------------------------------------------------------------- 00365 MsgStream log(messageService(), name()); 00366 log << MSG::INFO << "CosmicGenerator executing" << endreq; 00367 00368 00369 ++m_events; 00370 // MsgStream log(messageService(), name()); 00371 log << MSG::DEBUG << "Event #" << m_events << endreq; 00372 00373 00374 // clear up the vectors 00375 m_fourPos.clear(); 00376 m_fourMom.clear(); 00377 m_polarization.clear(); 00378 m_pdgCode.clear(); 00379 00380 00381 if(m_readfile) 00382 { 00383 if(!m_ffile.eof()) 00384 { 00385 CosmicEventParser evt; 00386 m_ffile >> evt; 00387 00388 std::cout << evt << std::endl; 00389 00390 double polx = 0; 00391 double poly = 0; 00392 double polz = 0; 00393 Polarization thePolarization; 00394 thePolarization.set_normal3d(HepNormal3D(polx,poly,polz)); 00395 m_polarization.push_back(thePolarization); 00396 00397 // 00398 // units are already converted to MeV's and mm. 00399 // 00400 m_fourPos.push_back(evt.Vertex()); 00401 m_fourMom.push_back(evt.Momentum()); 00402 m_pdgCode.push_back(evt.pdgID()); 00403 00404 } 00405 else 00406 { 00407 log << MSG::FATAL << "End of file reached - stop " << endreq; 00408 exit(1); 00409 return StatusCode::FAILURE; 00410 } 00411 } 00412 00413 else 00414 { 00415 00416 bool accepted=false; 00417 bool projected=false; 00418 HepLorentzVector pp; 00419 CosmicGun* gun = CosmicGun::GetCosmicGun(); 00420 HepLorentzVector vert; 00421 HepLorentzVector vert_proj; 00422 while(!accepted){ 00423 m_tried++; 00424 vert = generateVertex(); 00425 Hep3Vector vert3(vert.x(),vert.y(),vert.z()); 00426 00427 pp = gun->GenerateEvent(); 00428 if(m_dumpMC==1) { 00429 m_cosmicE=pp.e()*m_GeV; 00430 m_cosmicTheta=pp.theta(); 00431 m_cosmicPhi=pp.phi(); 00432 m_cosmicCharge=gun->GetMuonCharge(); 00433 m_tuple->write(); 00434 } 00435 double theta1=pp.theta(); 00436 double phi1=pp.phi(); 00437 double mag1=pp.vect().mag(); 00438 00439 Hep3Vector pp_corr(mag1*sin(theta1)*cos(phi1), 00440 -mag1*cos(theta1), 00441 mag1*sin(theta1)*sin(phi1)); 00442 Hep3Vector direction(pp_corr.x(),pp_corr.y(), pp_corr.z()); 00443 Hep3Vector center_dir=m_center-vert3; 00444 // if optimization activated, check for the direction of the generated muon 00445 00446 00447 if (m_cubeProj) { 00448 double alpha,beta; 00449 if(m_sphereOpt){ 00450 00451 beta=direction.angle(center_dir); 00452 alpha=asin(m_radius/center_dir.mag()); 00453 if(fabs(beta)>alpha){ 00454 log<<MSG::DEBUG<<alpha<<", "<<beta<<" rejected muon due to sphere cut! "<<endreq; 00455 m_rejected++; 00456 continue; 00457 } 00458 } 00459 00460 double l_fight0,l_fight1, l_fight2; 00461 double coor_x0, coor_y0, coor_z0; 00462 double coor_x1, coor_y1, coor_z1; 00463 double coor_x2, coor_y2, coor_z2; 00464 bool isIn0=false; 00465 bool isIn1=false; 00466 bool isIn2=false; 00467 HepPoint3D vert_pro0; 00468 HepPoint3D vert_pro1; 00469 HepPoint3D vert_pro2; 00470 HepPoint3D vert_org(vert3); 00471 00472 coor_y0 = m_ypos; // defined in jobOpt. 00473 coor_x0 = direction.x()*(coor_y0 - vert.y())/direction.y() +vert.x(); 00474 coor_z0 = direction.z()*(coor_y0 - vert.y())/direction.y() +vert.z(); 00475 00476 00477 00478 if( fabs(coor_x0) <= m_xpos && fabs(coor_z0) <= m_zpos){ 00479 vert_pro0=HepPoint3D (coor_x0, coor_y0, coor_z0); 00480 l_fight0=vert_org.distance(vert_pro0); 00481 isIn0=true; 00482 } 00483 else{ 00484 00485 coor_z1 = sign(coor_z0-vert.z())*m_zpos; // defined in jobOpt. 00486 coor_x1 = direction.x()*(coor_z1 - vert.z())/direction.z() +vert.x(); 00487 coor_y1 = direction.y()*(coor_z1 - vert.z())/direction.z() +vert.y(); 00488 00489 vert_pro1=HepPoint3D (coor_x1, coor_y1, coor_z1); 00490 l_fight1=vert_org.distance(vert_pro1); 00491 00492 00493 coor_x2 = sign(coor_x0-vert.x())*m_xpos; // defined in jobOpt. 00494 coor_z2 = direction.z()*(coor_x2 - vert.x())/direction.x() +vert.z(); 00495 coor_y2 = direction.y()*(coor_x2 - vert.x())/direction.x() +vert.y(); 00496 vert_pro2=HepPoint3D (coor_x2, coor_y2, coor_z2); 00497 l_fight2=vert_org.distance(vert_pro2); 00498 if(l_fight1<l_fight2) 00499 isIn1=true; 00500 else 00501 isIn2=true; 00502 } 00503 //reset vert(x,y,z,t), after projection 00504 if(isIn0){ 00505 vert_proj=HepLorentzVector (vert_pro0.x(),vert_pro0.y(),vert_pro0.z() , l_fight0); 00506 00507 projected =true; 00508 accepted = true; 00509 m_accepted++; 00510 } 00511 else if(isIn1){ 00512 vert_proj=HepLorentzVector (vert_pro1.x(),vert_pro1.y(),vert_pro1.z() , l_fight1); 00513 00514 projected =true; 00515 accepted = true; 00516 m_accepted++; 00517 } 00518 else if(isIn2){ 00519 vert_proj=HepLorentzVector (vert_pro2.x(),vert_pro2.y(),vert_pro2.z() , l_fight2); 00520 00521 projected =true; 00522 accepted = true; 00523 m_accepted++; 00524 } else { 00525 log<<MSG::DEBUG<<" rejected muon due to cube projection request! "<<endreq; 00526 m_rejected++; 00527 } 00528 00529 00530 00531 00532 } 00533 else accepted=true; // if no opt required accept the first muon 00534 } 00535 00536 // pp.setX(pp.x()*m_GeV); 00537 // pp.setY(pp.y()*m_GeV); 00538 // pp.setZ(pp.z()*m_GeV); 00539 // pp.setT(pp.t()*m_GeV); 00540 00541 pp.setX(pp.x()); 00542 pp.setY(pp.y()); 00543 pp.setZ(pp.z()); 00544 pp.setT(pp.t()); 00545 // Get the mass of the particle to be generated 00546 int charge = gun->GetMuonCharge(); 00547 // m_pdgCode.push_back(charge*13); 00548 m_pdgCode.push_back(charge*-13); 00549 00550 // const ParticleData* pdtparticle = m_particleTable->particle(ParticleID(abs(m_pdgCode.back()))); 00551 // double mass = pdtparticle->mass().value(); 00552 // double mass =105.5658; 00553 00554 // Compute the kinematic values. First, the vertex 4-vector: 00555 00556 00557 // Set the polarization. Realistically, this is going to be zero 00558 // for most studies, but you never know... 00559 double polx = 0; 00560 double poly = 0; 00561 double polz = 0; 00562 //m_polarization.set_normal3d(HepNormal3D(polx,poly,polz)); 00563 Polarization thePolarization; 00564 00565 // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ? 00566 // if not...do nothing...if so, invert position of y- and z- coordinate 00567 // 00568 // well and don't forget about the direction of the incoming cosmic muon(s) either 00569 // that means: y -> -y 00570 // 00571 if(!m_swapYZAxis){ 00572 // thePolarization.set_normal3d(HepNormal3D(polx,-poly,polz)); 00573 thePolarization.set_normal3d(HepNormal3D(polx,poly,polz)); 00574 } 00575 else 00576 thePolarization.set_normal3d(HepNormal3D(polx,polz,-poly)); 00577 00578 m_polarization.push_back(thePolarization); 00579 00580 00581 // The method of calculating e, theta, and phi depends on the user's 00582 // commands. Let the KinematicManager handle it. 00583 double e = pp.e(); 00584 double theta = pp.theta(); 00585 double phi = pp.phi(); 00586 00587 // At this point, we have e, theta, and phi. Put them together to 00588 // get the four-momentum. 00589 00590 double p2 = e*e - mass*mass; 00591 if ( p2 < 0 ) 00592 { 00593 log << MSG::ERROR 00594 << "Event #" << m_events 00595 << " E=" << e << ", mass=" << mass 00596 << " -- you have generated a tachyon! Increase energy or change particle ID." 00597 << endmsg; 00598 return StatusCode::FAILURE; 00599 } 00600 00601 double p = sqrt(p2); 00602 double px = p*sin(theta)*cos(phi); 00603 double pz = p*sin(theta)*sin(phi); 00604 double py = -p*cos(theta); 00605 00606 00607 00608 // Do we need to swap Y- and Z-axis for the PixelEndCap C Cosmic test ? 00609 // if not...do nothing...if so, invert position of y- and z- coordinate 00610 // 00611 // well and don't forget about the direction of the incoming cosmic muon(s) either 00612 // that means: y -> -y 00613 // 00614 if(!m_swapYZAxis) { 00615 // Line below corrupted py sign and forces muons to be upwards, not downwards. 00616 // m_fourMom.push_back(HepLorentzVector(px,-py,pz,pp.e())); 00617 m_fourMom.push_back(HepLorentzVector(px,py,pz,pp.e())); 00618 } 00619 else 00620 m_fourMom.push_back(HepLorentzVector(px,pz,-py,pp.e())); 00621 double x = vert.x(); 00622 double y = vert.y(); 00623 double z = vert.z(); 00624 double t = vert.t(); 00625 if(projected){ 00626 00627 x = vert_proj.x(); 00628 y = vert_proj.y(); 00629 z = vert_proj.z(); 00630 t = vert.t()-vert_proj.t()/HepLorentzVector(px,py,pz,pp.e()).beta() 00631 +m_tcor; 00632 00633 00634 00635 } 00636 // log << MSG::DEBUG 00637 // << "proj:"<<m_cubeProj<<" ("<<x<<", "<<y<<", "<<z<<", "<<t<<")" 00638 // <<endreq; 00639 // Do we need to swap Y- and Z-axis for the PixelEndCap A Cosmic test ? 00640 // if not...do nothing...if so, invert position of y- and z- coordinate 00641 // 00642 // but not only that...change also the direction of the incoming cosmic muon(s), 00643 // they must go towards the pixel endcap A, i.e. y -> -y 00644 // 00645 00646 00647 if(!m_swapYZAxis) 00648 m_fourPos.push_back(HepLorentzVector(x,y,z,t)); 00649 else 00650 m_fourPos.push_back(HepLorentzVector(x,z,y,t)); 00651 00652 log << MSG::DEBUG 00653 << " (x,y,z,t) = (" 00654 << m_fourPos.back().x() << "," 00655 << m_fourPos.back().y() << "," 00656 << m_fourPos.back().z() << "," 00657 << m_fourPos.back().t() << "), (Px,Py,Pz,E) = (" 00658 << m_fourMom.back().px() << "," 00659 << m_fourMom.back().py() << "," 00660 << m_fourMom.back().pz() << "," 00661 << m_fourMom.back().e() << ")" 00662 << endreq; 00663 log << MSG::DEBUG 00664 << " (theta,phi) = (" << theta << "," << phi << "), " 00665 << "polarization(x,y,z) = (" 00666 << m_polarization.back().normal3d().x() << "," 00667 << m_polarization.back().normal3d().y() << "," 00668 << m_polarization.back().normal3d().z() << ")" 00669 << endreq; 00670 if(m_dumpMC==1) { 00671 // m_cosmicE=pp.e()*m_GeV; 00672 // m_cosmicTheta=pp.theta(); 00673 // m_cosmicPhi=pp.phi(); 00674 // m_cosmicCharge=gun->GetMuonCharge(); 00675 mc_theta=m_fourMom.back().theta(); 00676 mc_phi=m_fourMom.back().phi(); 00677 mc_px=m_fourMom.back().px(); 00678 mc_py=m_fourMom.back().py(); 00679 mc_pz=m_fourMom.back().pz(); 00680 // m_cosmicE=pz; 00681 m_tuple1->write(); 00682 } 00683 } 00684 00685 // fill evt 00686 00687 GenEvent* event = new GenEvent(1,1); 00688 00689 if(m_fourMom.size()==m_fourPos.size()&&m_fourMom.size()==m_polarization.size()){ 00690 00691 for(std::size_t v=0;v<m_fourMom.size();++v){ 00692 00693 // Note: The vertex and particle are owned by the event, so the 00694 // event is responsible for those pointers. 00695 00696 // Create the particle, and specify its polarization. 00697 00698 GenParticle* particle = new GenParticle( m_fourMom[v], m_pdgCode[v], 1); 00699 particle->set_polarization( m_polarization[v] ); 00700 00701 // Create the vertex, and add the particle to the vertex. 00702 GenVertex* vertex = new GenVertex(m_fourPos[v]); 00703 vertex->add_particle_out( particle ); 00704 00705 // Add the vertex to the event. 00706 event->add_vertex( vertex ); 00707 00708 } 00709 00710 event->set_event_number(m_events); // Set the event number 00711 00712 00713 00714 } else { 00715 00716 log<<MSG::ERROR<<"Wrong different number of vertexes/momenta/polaritazions!"<<endreq; 00717 return StatusCode::FAILURE; 00718 } 00719 // Check if the McCollection already exists 00720 SmartDataPtr<McGenEventCol> anMcCol(eventSvc(), "/Event/Gen"); 00721 if (anMcCol!=0) 00722 { 00723 // Add event to existing collection 00724 00725 log << MSG::INFO << "Add McGenEvent to existing collection" << endreq; 00726 McGenEvent* mcEvent = new McGenEvent(event); 00727 anMcCol->push_back(mcEvent); 00728 } 00729 else 00730 { 00731 // Create Collection and add to the transient store 00732 McGenEventCol *mcColl = new McGenEventCol; 00733 McGenEvent* mcEvent = new McGenEvent(event); 00734 mcColl->push_back(mcEvent); 00735 StatusCode sc = eventSvc()->registerObject("/Event/Gen",mcColl); 00736 if (sc != StatusCode::SUCCESS) 00737 { 00738 log << MSG::ERROR << "Could not register McGenEvent" << endreq; 00739 delete mcColl; 00740 delete event; 00741 delete mcEvent; 00742 return StatusCode::FAILURE; 00743 } 00744 00745 } 00746 00747 00748 00749 00750 00751 return StatusCode::SUCCESS; 00752 00753 }
bool CosmicGenerator::exzCut | ( | const Hep3Vector & | pos, | |
const HepLorentzVector & | p | |||
) | [private] |
Definition at line 789 of file CosmicGenerator.cxx.
References cut, m_GeV, and m_rmax.
00790 { 00791 // p is in GeV... 00792 00793 double r =0; 00794 bool cut = false; 00795 if(pos.z()<0){ 00796 r = sqrt((pow(pos.x(),2)+pow(pos.z()+28000,2))) ; 00797 double e = 0.45238*r+5000 ; 00798 cut = p.e()*m_GeV>e; 00799 } 00800 else 00801 { 00802 r = sqrt((pow(pos.x(),2)+pow(pos.z()-20000,2))) ; 00803 if(r<15000) { 00804 cut = true; 00805 } else 00806 { 00807 double e = 0.461538*(r-15000)+10000 ; 00808 cut = p.e()*m_GeV>e; 00809 // std::cout<<"z>0 r , e, p.e = "<<r <<" " <<e <<" " <<p.e()*m_GeV<<std::endl; 00810 } 00811 } 00812 00813 00814 cut = cut && r < m_rmax ; 00815 00816 return cut; 00817 00818 00819 }
StatusCode CosmicGenerator::finalize | ( | ) |
Definition at line 757 of file CosmicGenerator.cxx.
References Bes_Common::INFO, m_accepted, m_cubeProj, m_rejected, and m_tried.
00757 { 00758 //--------------------------------------------------------------------------- 00759 // Get the KinematicManager. 00760 00761 if(m_cubeProj) { 00762 00763 MsgStream log(messageService(), name()); 00764 00765 log << MSG::INFO<<"***************************************************"<<endreq; 00766 log << MSG::INFO <<"** you have run CosmicGenerator with some " <<endreq; 00767 log << MSG::INFO <<"** filters for cosmic muon simulation" <<endreq; 00768 log << MSG::INFO <<"** "<<m_tried<<" muons were tried" <<endreq; 00769 log << MSG::INFO <<"** "<<m_accepted<<" muons were accepted" <<endreq; 00770 log << MSG::INFO <<"** "<<m_rejected<<" muons were rejected" <<endreq; 00771 log << MSG::INFO<<"***************************************************"<<endreq; 00772 std::cout<<"***************************************************"<<std::endl; 00773 std::cout <<"** you have run CosmicGenerator with some " <<std::endl; 00774 std::cout <<"** filters for cosmic muon simulation" <<std::endl; 00775 std::cout <<"** "<<m_tried<<" muons were tried" <<std::endl; 00776 std::cout <<"** "<<m_accepted<<" muons were accepted" <<std::endl; 00777 std::cout <<"** "<<m_rejected<<" muons were rejected" <<std::endl; 00778 std::cout<<"***************************************************"<<std::endl; 00779 00780 } 00781 00782 return StatusCode::SUCCESS; 00783 }
HepLorentzVector CosmicGenerator::generateVertex | ( | void | ) |
Definition at line 334 of file CosmicGenerator.cxx.
References Bes_Common::FATAL, IBesRndmGenSvc::GetEngine(), m_tmax, m_tmin, m_xhig, m_xlow, m_yval, m_zhig, m_zlow, and p_BesRndmGenSvc.
Referenced by execute().
00334 { 00335 00336 MsgStream log(messageService(), name()); 00337 00338 // Get the pointer to the engine of the stream named SINGLE. If the 00339 // stream does not exist is created automaticaly 00340 // HepRandomEngine* engine = CosmicGenerator::p_AtRndmGenSvc->GetEngine("COSMICS"); 00341 HepRandomEngine* engine = CosmicGenerator::p_BesRndmGenSvc->GetEngine("PYTHIA"); 00342 // Generate a random number according to the distribution. 00343 00344 float x_val = RandFlat::shoot(engine, m_xlow, m_xhig); 00345 float z_val = RandFlat::shoot(engine, m_zlow, m_zhig); 00346 00347 // Generate a random number for time offset 00348 float t_val=0; 00349 if(m_tmin < m_tmax){ 00350 t_val = RandFlat::shoot(engine, m_tmin, m_tmax); 00351 } 00352 else if(m_tmin == m_tmax){ 00353 t_val = m_tmin; 00354 } 00355 else log << MSG::FATAL << " You specified m_tmin = " << m_tmin << " and m_tmax " << m_tmax << endreq; 00356 HepLorentzVector p(x_val,m_yval,z_val, t_val*c_light); 00357 00358 return p; 00359 00360 }
StatusCode CosmicGenerator::initialize | ( | ) |
Definition at line 213 of file CosmicGenerator.cxx.
References calibUtil::ERROR, Bes_Common::FATAL, CosmicGun::GetCosmicGun(), Bes_Common::INFO, CosmicGun::InitializeGenerator(), m_accepted, m_center, m_cosmicCharge, m_cosmicE, m_cosmicPhi, m_cosmicTheta, m_ctcut, m_dumpMC, m_emax, m_emin, m_events, m_ffile, m_GeV, m_infile, m_IPx, m_IPy, m_IPz, m_mm, m_printEvent, m_printMod, m_radius, m_readfile, m_rejected, m_tcor, m_tried, m_tuple, m_tuple1, m_xhig, m_xlow, m_xpos, m_ypos, m_yval, m_zhig, m_zlow, m_zpos, mass, mc_phi, mc_px, mc_py, mc_pz, mc_theta, ntupleSvc(), p_BesRndmGenSvc, p_msgSvc, CosmicGun::PrintLevel(), CosmicGun::SetCosCut(), and CosmicGun::SetEnergyRange().
00213 { 00214 //--------------------------------------------------------------------------- 00215 00216 // Initialize event count. 00217 00218 m_events = 0; 00219 m_tried=0; 00220 m_accepted=0; 00221 m_rejected=0; 00222 if(m_emin<0.1) {m_emin=0.1;std::cout<<"WARNING: set emin to 0.1 GeV"<<std::endl;} 00223 m_emin =m_emin *m_GeV; 00224 m_emax =m_emax *m_GeV; 00225 m_xlow =m_xlow *m_mm; 00226 m_xhig =m_xhig *m_mm; 00227 m_zlow =m_zlow *m_mm; 00228 m_zhig =m_zhig *m_mm; 00229 m_yval =m_yval *m_mm; 00230 m_xpos =m_xpos *m_mm; 00231 m_ypos =m_ypos *m_mm; 00232 m_zpos =m_zpos *m_mm; 00233 m_radius=m_radius*m_mm; 00234 m_tcor=m_tcor 00235 +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)) 00236 /(m_emin/sqrt(m_emin*m_emin+mass*mass*m_GeV*m_GeV)); 00237 00238 ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap.h 00239 p_msgSvc = 0; 00240 StatusCode result = svcLocator->service("MessageSvc", p_msgSvc); 00241 00242 if ( !result.isSuccess() || p_msgSvc == 0 ) 00243 { 00244 std::cerr << "VGenCommand(): could not initialize MessageSvc!" << std::endl; 00245 throw GaudiException("Could not initalize MessageSvc","CosmicGenerator",StatusCode::FAILURE); 00246 } 00247 00248 MsgStream log (p_msgSvc,"CosmicGenerator"); 00249 if(m_dumpMC==1) { 00250 StatusCode status; 00251 NTuplePtr nt(ntupleSvc(), "FILE1/comp"); 00252 if(nt) m_tuple = nt; 00253 else { 00254 m_tuple = ntupleSvc()->book ("FILE1/comp", CLID_ColumnWiseTuple, "ks N-Tuple example"); 00255 if ( m_tuple ) { 00256 status = m_tuple->addItem ("cosmic_e", m_cosmicE); 00257 status = m_tuple->addItem ("cosmic_theta", m_cosmicTheta); 00258 status = m_tuple->addItem ("cosmic_phi", m_cosmicPhi); 00259 status = m_tuple->addItem ("cosmic_charge", m_cosmicCharge); 00260 } 00261 else { 00262 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple) << endmsg; 00263 return StatusCode::FAILURE; 00264 } 00265 } 00266 NTuplePtr nt1(ntupleSvc(), "FILE1/compp"); 00267 if(nt1) m_tuple1 = nt1; 00268 else { 00269 m_tuple1 = ntupleSvc()->book ("FILE1/compp", CLID_ColumnWiseTuple, "ks N-Tuple example"); 00270 if ( m_tuple1 ) { 00271 status = m_tuple1->addItem ("mc_theta", mc_theta); 00272 status = m_tuple1->addItem ("mc_phi", mc_phi); 00273 status = m_tuple1->addItem ("mc_px", mc_px); 00274 status = m_tuple1->addItem ("mc_py", mc_py); 00275 status = m_tuple1->addItem ("mc_pz", mc_pz); 00276 } 00277 else { 00278 log << MSG::ERROR << " Cannot book N-tuple:" << long(m_tuple1) << endmsg; 00279 return StatusCode::FAILURE; 00280 } 00281 } 00282 } 00283 if(m_infile=="NONE") 00284 00285 { 00286 // Initialize the BES random number services. 00287 00288 static const bool CREATEIFNOTTHERE(true); 00289 StatusCode RndmStatus = svcLocator->service("BesRndmGenSvc", p_BesRndmGenSvc, CREATEIFNOTTHERE); 00290 00291 if (!RndmStatus.isSuccess() || 0 == p_BesRndmGenSvc) 00292 { 00293 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq; 00294 return RndmStatus; 00295 } 00296 00297 CosmicGun* gun = CosmicGun::GetCosmicGun(); 00298 00299 gun->SetEnergyRange(m_emin/m_GeV,m_emax/m_GeV); 00300 gun->SetCosCut(m_ctcut); 00301 gun->PrintLevel(m_printEvent, m_printMod); 00302 float flux_withCT = gun->InitializeGenerator(); 00303 00304 log << MSG::INFO << "Initialisation cosmic gun done." << endreq; 00305 log << MSG::INFO << "Accepted diff flux after E and cos(theta) cuts = " << 00306 flux_withCT << " /cm^2/s" << endreq; 00307 log << MSG::INFO << "Accepted total flux after E and cos(theta) cuts = " << 00308 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << endreq; 00309 00310 std::cout<< "Accepted diff flux after E and cos(theta) cuts = " << 00311 flux_withCT << " /cm^2/s" << std::endl; 00312 std::cout << "Accepted total flux after E and cos(theta) cuts = " << 00313 flux_withCT*(m_xhig-m_xlow)/m_mm*(m_zhig-m_zlow)/m_mm << " /s" << std::endl; 00314 00315 } 00316 else 00317 { 00318 log << MSG::INFO << "Cosmics are read from file " << m_infile << endreq; 00319 m_ffile.open(m_infile.c_str()); 00320 if(!m_ffile) 00321 { 00322 log << MSG::FATAL << "Could not open input file - stop! " << endreq; 00323 return StatusCode::FAILURE; 00324 } 00325 m_readfile = true; 00326 } 00327 00328 m_center=Hep3Vector(m_IPx, m_IPy, m_IPz); 00329 log << MSG::INFO << "Initialisation done" << endreq; 00330 return StatusCode::SUCCESS; 00331 00332 }
long int CosmicGenerator::m_accepted [private] |
Definition at line 81 of file CosmicGenerator.h.
Referenced by execute(), finalize(), and initialize().
ActiveStoreSvc* CosmicGenerator::m_activeStore [private] |
Definition at line 67 of file CosmicGenerator.h.
Hep3Vector CosmicGenerator::m_center [private] |
NTuple::Item<double> CosmicGenerator::m_cosmicCharge [private] |
NTuple::Item<double> CosmicGenerator::m_cosmicE [private] |
NTuple::Item<double> CosmicGenerator::m_cosmicPhi [private] |
NTuple::Item<double> CosmicGenerator::m_cosmicTheta [private] |
float CosmicGenerator::m_ctcut [private] |
bool CosmicGenerator::m_cubeProj [private] |
Definition at line 92 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and finalize().
int CosmicGenerator::m_dumpMC [private] |
Definition at line 95 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
float CosmicGenerator::m_emax [private] |
float CosmicGenerator::m_emin [private] |
long int CosmicGenerator::m_events [private] |
bool CosmicGenerator::m_exzCut [private] |
std::ifstream CosmicGenerator::m_ffile [private] |
std::vector<HepLorentzVector> CosmicGenerator::m_fourMom [private] |
std::vector<HepLorentzVector> CosmicGenerator::m_fourPos [private] |
float CosmicGenerator::m_GeV [private] |
Definition at line 98 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), exzCut(), and initialize().
std::string CosmicGenerator::m_infile [private] |
Definition at line 103 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), and initialize().
float CosmicGenerator::m_IPx [private] |
float CosmicGenerator::m_IPy [private] |
float CosmicGenerator::m_IPz [private] |
float CosmicGenerator::m_mm [private] |
std::vector<int> CosmicGenerator::m_pdgCode [private] |
float CosmicGenerator::m_phimax [private] |
float CosmicGenerator::m_phimin [private] |
std::vector<HepMC::Polarization> CosmicGenerator::m_polarization [private] |
long int CosmicGenerator::m_printEvent [private] |
long int CosmicGenerator::m_printMod [private] |
float CosmicGenerator::m_radius [private] |
Definition at line 86 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
bool CosmicGenerator::m_readfile [private] |
Definition at line 102 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
long int CosmicGenerator::m_rejected [private] |
Definition at line 81 of file CosmicGenerator.h.
Referenced by execute(), finalize(), and initialize().
float CosmicGenerator::m_rmax [private] |
long int CosmicGenerator::m_selection [private] |
Definition at line 94 of file CosmicGenerator.h.
StoreGateSvc* CosmicGenerator::m_sgSvc [private] |
Definition at line 66 of file CosmicGenerator.h.
bool CosmicGenerator::m_sphereOpt [private] |
bool CosmicGenerator::m_swapYZAxis [private] |
float CosmicGenerator::m_tcor [private] |
Definition at line 88 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
float CosmicGenerator::m_thetamax [private] |
float CosmicGenerator::m_thetamin [private] |
float CosmicGenerator::m_tmax [private] |
Definition at line 88 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), and generateVertex().
float CosmicGenerator::m_tmin [private] |
Definition at line 88 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), and generateVertex().
long int CosmicGenerator::m_tried [private] |
Definition at line 81 of file CosmicGenerator.h.
Referenced by execute(), finalize(), and initialize().
NTuple::Tuple* CosmicGenerator::m_tuple [private] |
Definition at line 68 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
NTuple::Tuple* CosmicGenerator::m_tuple1 [private] |
Definition at line 73 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
float CosmicGenerator::m_xhig [private] |
Definition at line 85 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), generateVertex(), and initialize().
float CosmicGenerator::m_xlow [private] |
Definition at line 85 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), generateVertex(), and initialize().
float CosmicGenerator::m_xpos [private] |
Definition at line 86 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
float CosmicGenerator::m_ypos [private] |
Definition at line 86 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
float CosmicGenerator::m_yval [private] |
Definition at line 85 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), generateVertex(), and initialize().
float CosmicGenerator::m_zhig [private] |
Definition at line 85 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), generateVertex(), and initialize().
float CosmicGenerator::m_zlow [private] |
Definition at line 85 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), generateVertex(), and initialize().
float CosmicGenerator::m_zpos [private] |
Definition at line 86 of file CosmicGenerator.h.
Referenced by CosmicGenerator(), execute(), and initialize().
NTuple::Item<double> CosmicGenerator::mc_phi [private] |
NTuple::Item<double> CosmicGenerator::mc_px [private] |
NTuple::Item<double> CosmicGenerator::mc_py [private] |
NTuple::Item<double> CosmicGenerator::mc_pz [private] |
NTuple::Item<double> CosmicGenerator::mc_theta [private] |
IBesRndmGenSvc * CosmicGenerator::p_BesRndmGenSvc = 0 [static] |
Definition at line 62 of file CosmicGenerator.h.
Referenced by cosmicrndm_(), generateVertex(), and initialize().
IMessageSvc* CosmicGenerator::p_msgSvc [private] |