CosmicGenerator Class Reference

#include <CosmicGenerator.h>

List of all members.

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

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


Detailed Description

Definition at line 51 of file CosmicGenerator.h.


Constructor & Destructor Documentation

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 (  ) 

Definition at line 208 of file CosmicGenerator.cxx.

00210 {}


Member Function Documentation

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 }


Member Data Documentation

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]

Definition at line 109 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::m_cosmicCharge [private]

Definition at line 72 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::m_cosmicE [private]

Definition at line 69 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::m_cosmicPhi [private]

Definition at line 71 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::m_cosmicTheta [private]

Definition at line 70 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

float CosmicGenerator::m_ctcut [private]

Definition at line 84 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

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]

Definition at line 83 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

float CosmicGenerator::m_emin [private]

Definition at line 83 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

long int CosmicGenerator::m_events [private]

Definition at line 81 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

bool CosmicGenerator::m_exzCut [private]

Definition at line 119 of file CosmicGenerator.h.

Referenced by CosmicGenerator().

std::ifstream CosmicGenerator::m_ffile [private]

Definition at line 104 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

std::vector<HepLorentzVector> CosmicGenerator::m_fourMom [private]

Definition at line 108 of file CosmicGenerator.h.

Referenced by execute().

std::vector<HepLorentzVector> CosmicGenerator::m_fourPos [private]

Definition at line 107 of file CosmicGenerator.h.

Referenced by execute().

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]

Definition at line 85 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

float CosmicGenerator::m_IPy [private]

Definition at line 85 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

float CosmicGenerator::m_IPz [private]

Definition at line 85 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

float CosmicGenerator::m_mm [private]

Definition at line 99 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

std::vector<int> CosmicGenerator::m_pdgCode [private]

Definition at line 82 of file CosmicGenerator.h.

Referenced by execute().

float CosmicGenerator::m_phimax [private]

Definition at line 96 of file CosmicGenerator.h.

Referenced by CosmicGenerator().

float CosmicGenerator::m_phimin [private]

Definition at line 96 of file CosmicGenerator.h.

Referenced by CosmicGenerator().

std::vector<HepMC::Polarization> CosmicGenerator::m_polarization [private]

Definition at line 110 of file CosmicGenerator.h.

Referenced by execute().

long int CosmicGenerator::m_printEvent [private]

Definition at line 93 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

long int CosmicGenerator::m_printMod [private]

Definition at line 93 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and initialize().

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]

Definition at line 121 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and exzCut().

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]

Definition at line 89 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and execute().

bool CosmicGenerator::m_swapYZAxis [private]

Definition at line 91 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), and execute().

float CosmicGenerator::m_tcor [private]

Definition at line 88 of file CosmicGenerator.h.

Referenced by CosmicGenerator(), execute(), and initialize().

float CosmicGenerator::m_thetamax [private]

Definition at line 96 of file CosmicGenerator.h.

Referenced by CosmicGenerator().

float CosmicGenerator::m_thetamin [private]

Definition at line 96 of file CosmicGenerator.h.

Referenced by CosmicGenerator().

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]

Definition at line 74 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::mc_px [private]

Definition at line 76 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::mc_py [private]

Definition at line 77 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::mc_pz [private]

Definition at line 78 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

NTuple::Item<double> CosmicGenerator::mc_theta [private]

Definition at line 75 of file CosmicGenerator.h.

Referenced by execute(), and initialize().

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]

Definition at line 113 of file CosmicGenerator.h.

Referenced by initialize().


Generated on Tue Nov 29 23:18:09 2016 for BOSS_7.0.2 by  doxygen 1.4.7