BOSS6.5.5: BesHepMCInterface Class Reference

00013 {
00014   p_evtSvc = 0;
00015 }

00017 {
00018         //std::cout<< "\b the BesHepMCInterface is being destroyed "<<std::endl;
00019 }

00148 {
00149   //suppose relavtive vectoy is v(vx,vy,vz),position are (x,y,z,t)in CMS and(X,Y,Z,T)in Lab,the formulae for boost vertex position from Cms to Lab in natural units are following:
00150        // v2 = vx*vx + vy*vy + vz*vz;  
00151        // gamma = 1.0 / sqrt(1.0 - v2);
00152        // bp = vx*x + vy*y + vz*z;
00153        // X=x + gamma2*bp*vx + gamma*vx*t;
00154        // Y=y + gamma2*bp*vy + gamma*vy*t;
00155        // Z=z + gamma2*bp*vz + gamma*vz*t;
00156        // T=gamma*(t + bp);
00157   // the function of these formulae is the same as xvtx.boost(ThreeVector v)      
00158   //For the E_cms,need to loop and splice the momentum of all the final state particles
00159   G4LorentzVector ptot=0;
00160   double E_cms,p_Mag,e_Mass2,M_cms,theta;
00161   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00162       vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
00163     
00164     for(HepMC::GenVertex::particle_iterator
00165           vpitr = (*vitr)->particles_begin(HepMC::children);
00166         vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
00167       
00168       if((*vpitr)->status() !=1)continue;// Not final state particle
00169       
00170       G4LorentzVector p= (*vpitr)-> momentum();
00171       ptot = p + ptot; 
00172     }
00173   }
00174   E_cms=ptot.e()*GeV;
00175   
00176   //get G4Svc, allow user to access the beam momentum shift in JobOptions
00177     ISvcLocator* svcLocator = Gaudi::svcLocator();
00178   IG4Svc* tmpSvc;
00179   G4Svc* m_G4Svc;
00180   StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
00181   m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00182   
00183  
00184   G4ThreeVector v=0;   //velocity of cms as seen from lab
00185   //for cms velocity 
00186   theta=11*mrad;     
00187   //theta=(m_G4Svc->GetBeamAngle())*mrad;
00188   //theta is half of the angle between the two beams,namely,is the angle between the beam and Z axis.
00189   M_cms=E_cms;  //for p_cms=0 in the CMS
00190   e_Mass2=electron_mass_c2*electron_mass_c2; //square of electron mass 
00191   p_Mag=sqrt((E_cms*E_cms-4*e_Mass2)/(2*(1-cos(pi-2*theta))));
00192   G4ThreeVector p_Positron(p_Mag*sin(theta),0,p_Mag*cos(theta));
00193   G4ThreeVector p_Electron(p_Mag*sin(pi-theta),0,p_Mag*cos(pi-theta));
00194   G4ThreeVector vee=p_Positron+p_Electron;
00195   G4ThreeVector beamshift(m_G4Svc->GetBeamShiftPx(),m_G4Svc->GetBeamShiftPy(),m_G4Svc->GetBeamShiftPz());
00196   if(m_G4Svc-> GetSetBeamShift()==true && fabs(M_cms-3686)<50.0) {vee = beamshift;} 
00197 
00198     double pee=vee.r(); 
00199     M_cms = sqrt(M_cms*M_cms + pee*pee);
00200     v=vee/M_cms;
00201 
00202   
00203   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00204       vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
00205     
00206     // Boost vertex position from cms to lab 
00207     G4LorentzVector xvtx= (*vitr)-> position();
00208     xvtx.boost(v);
00209     (*vitr)->set_position(xvtx);
00210     
00211     // Boost the particles four momentum from cms to lab
00212     // the transform formulae are similary to the position boost       
00213     for (HepMC::GenVertex::particle_iterator 
00214            vpitr = (*vitr)->particles_begin(HepMC::children);
00215          vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
00216       G4LorentzVector p= (*vpitr)-> momentum();
00217       p=p.boost(v);
00218       (*vpitr)->set_momentum(p);
00219     }
00220   }
00221 }

00122 {
00123   G4int flag =0;
00124   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00125           vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
00126     for (HepMC::GenVertex::particle_iterator
00127              pitr= (*vitr)->particles_begin(HepMC::children);
00128                 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00129       if((*pitr)->end_vertex())
00130         {flag = 1; break;}
00131     }
00132     if(flag) break;
00133   }
00134   if(m_logLevel <= 4)
00135   {
00136     if(flag == 0)
00137       G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
00138     else
00139       G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
00140   }
00141   return flag;
00142 }

00022 {
00023   
00024   if (p_evtSvc == 0) {     
00025     //std::cout<<" standard interface to EvtDataSvc for retrieving HepMC events"<<std::endl;
00026     ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap
00027     StatusCode sc=svcLocator->service("EventDataSvc", p_evtSvc);
00028     if (sc.isFailure()) 
00029     {
00030             //std::cout<<"BesHepMCInterface could not access EventDataSvc!!"<<std::endl;
00031     }
00032   }
00033         int n = 0;
00034   //std::cout <<" BesHepMCInterface::GenerateAnEvent "<<std::endl;
00035 
00036         SmartDataPtr<McGenEventCol> mcCollptr( p_evtSvc, "/Event/Gen");
00037 
00038         if ( mcCollptr != 0 ) {
00039           //std::cout <<" could retrieve the collection "<<std::endl;
00040 
00041           n = mcCollptr->size();
00042           
00043           //std::cout <<" nr of events "<<n<<std::endl;
00044 
00045           McGenEventCol::const_iterator it = mcCollptr->begin();
00046 
00047           McGenEvent* mcEvent = (McGenEvent* ) (*it);   
00048           // getting only the first event here.         
00049           // std::cout <<"iterator"<<std::endl;
00050                 
00051           HepMC::GenEvent *p_evt;
00052                 
00053           //std::cout <<"hepmc event"<<std::endl;
00054           p_evt = mcEvent->getGenEvt();
00055           //std::cout << "eventNumber = " << p_evt-> event_number()  << std::endl;
00056 
00057           //std::cout << " BesHepMCInterface:: --> " <<std::endl;
00058           //std::cout << " particles_size = " << p_evt->particles_size() 
00059                 //    << " vertices_size = " << p_evt->vertices_size()
00060                 //    << std::endl;
00061 
00062           //std::cout <<"got it"<<std::endl;
00063           
00064           return p_evt;
00065         }
00066 
00067     else {
00068           std::cout << "no McGenEventCollection found."  << std::endl;
00069         }     
00070           return 0;
00071 }

00438 {
00439   // delete previous event object
00440   delete hepmcEvent;
00441 
00442   // generate next event
00443   hepmcEvent= GenerateHepMCEvent();
00444   if(! hepmcEvent) {
00445     G4Exception("G4HepMCInterface","GeneratePrimaryVertex",RunMustBeAborted,
00446               "HepMCInterface: no generated particles. run terminated..." );
00447     return;
00448   }
00449   HepMC2G4(hepmcEvent, anEvent);
00450 }

00095 { 
00096   return hepmcEvent; 
00097 }

00055 {return m_logLevel;}

00055 {return m_logLevel;}

00227 {
00228   //Print HepMC Event before boost
00229   if(m_logLevel <=4 )
00230     Print(hepmcevt);
00231 
00232   //get G4Svc
00233   ISvcLocator* svcLocator = Gaudi::svcLocator();
00234   IG4Svc* tmpSvc;
00235   G4Svc* m_G4Svc;
00236   StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
00237   m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc); 
00238 
00239   //considering 22rmad beam clossing angle ,need to boost the Cms to Lab 
00240   // need to boost ?
00241   if(m_G4Svc->GetBoostLab()== true) 
00242     Boost(hepmcevt); 
00243   
00244   //Print HepMC Event after boost
00245   if(m_logLevel <=2  && m_G4Svc->GetBoostLab()== true)
00246     Print(hepmcevt);
00247 
00248   //*********************send particles to HPlist*************************//
00249   G4int flag = CheckType(hepmcevt);
00250   if(flag)
00251   {
00252     for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00253         vitr != hepmcevt->vertices_end(); ++vitr ) 
00254     { // loop for vertex ...
00255       G4int vtxFlag=0;
00256       G4int pOut = (*vitr)->particles_out_size();
00257       HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
00258       G4int pdgcode= (*pitr)-> pdg_id();
00259       if(pOut == 1 && pdgcode == -11)
00260         vtxFlag=1;
00261       G4int barcodeVtx = (*vitr)->barcode();
00262 
00263       for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
00264            pitr != (*vitr)->particles_end(HepMC::children); ++pitr) 
00265       {// loop for particle from this vertex 
00266         G4int pdgcode = (*pitr)->pdg_id();
00267         if(vtxFlag==0)
00268         {
00269           if(pdgcode==-22) pdgcode=22;
00270           G4LorentzVector p= (*pitr)-> momentum();    
00271           G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
00272           //G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
00273           //particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
00274           G4int ISTHEP = (*pitr)->status();
00275   
00276           G4int barcodeEndVtx = 0;
00277           if((*pitr)->end_vertex())
00278             barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00279   
00280           G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
00281           HPlist.push_back(hepmcParticle);
00282           
00283           for( size_t ii=0; ii<HPlist.size(); ii++ )
00284           {  
00285             if(barcodeVtx == HPlist[ii]->GetBarcodeEndVtx())
00286             {
00287               HPlist[ii]->GetTheParticle()->SetDaughter(particle); 
00288               hepmcParticle->Done();
00289               break; 
00290             }
00291           }
00292         }
00293       }
00294     }
00295   }
00296   else 
00297   {
00298     for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00299                 vitr != hepmcevt->vertices_end(); ++vitr )
00300     {
00301       for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
00302                      pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
00303       {
00304         G4int ISTHEP = (*pitr)->status();
00305         G4LorentzVector p= (*pitr)-> momentum();
00306         G4int pdgcode = (*pitr)->pdg_id();
00307         G4int barcodeEndVtx = 0;
00308         if((*pitr)->end_vertex())
00309           barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00310         G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
00311         //G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode);
00312         //particle->Set4Momentum(p.x()*GeV,p.y()*GeV,p.z()*GeV,p.t()*GeV);
00313 
00314         G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
00315         HPlist.push_back(hepmcParticle);
00316         if(ISTHEP>1)
00317           hepmcParticle->Done();
00318       }
00319     }
00320   }
00321 
00322   //print particles in HPlist
00323   if(m_logLevel<=4)
00324     PrintHPlist();
00325 
00326   //get time and position information from G4Svc
00327   G4double pmPosX,pmPosY,pmPosZ,pmTime;
00328   G4double tmpPosX,tmpPosY,tmpPosZ,tmpT;
00329   G4double beamPosX,beamPosY,beamPosZ,beamSizeX,beamSizeY,beamSizeZ;
00330 
00331   if(m_RealizationSvc->UseDBFlag() == false) {
00332     beamPosX = m_G4Svc->GetBeamPosX();
00333     beamPosY = m_G4Svc->GetBeamPosY();
00334     beamPosZ = m_G4Svc->GetBeamPosZ();
00335 
00336     beamSizeX = m_G4Svc->GetBeamSizeX();
00337     beamSizeY = m_G4Svc->GetBeamSizeY();
00338     beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
00339   } 
00340   if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == true) {
00341     beamPosX = m_RealizationSvc->getBunchPosX();
00342     beamPosY = m_RealizationSvc->getBunchPosY();
00343     beamPosZ = m_RealizationSvc->getBunchPosZ();
00344 
00345     beamSizeX = m_RealizationSvc->getBunchSizeX();
00346     beamSizeY = m_RealizationSvc->getBunchSizeY();
00347     beamSizeZ = m_RealizationSvc->getBunchSizeZ();
00348   }
00349   if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == false) {
00350     beamPosX = m_G4Svc->GetBeamPosX();
00351     beamPosY = m_G4Svc->GetBeamPosY();
00352     beamPosZ = m_G4Svc->GetBeamPosZ();
00353 
00354     beamSizeX = m_G4Svc->GetBeamSizeX();
00355     beamSizeY = m_G4Svc->GetBeamSizeY();
00356     beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
00357   }
00358 
00359   G4double gaussX = G4RandGauss::shoot();
00360   G4double gaussY = G4RandGauss::shoot();
00361   G4double gaussZ = G4RandGauss::shoot();
00362   G4double gaussT = G4RandGauss::shoot();
00363 
00364   G4double beamStartTime = m_G4Svc->GetBeamStartTime() * ns;
00365   G4double beamDeltaTime = m_G4Svc->GetBeamDeltaTime() * ns;
00366   G4double nBunch = m_G4Svc->GetNBunch();
00367   G4double bunchTimeSigma = m_G4Svc->GetBunchTimeSigma() * ns;
00368   
00369   G4double ran=G4UniformRand();
00370   G4double beamTime = bunchTimeSigma*G4RandGauss::shoot() + beamStartTime + beamDeltaTime*int(ran*nBunch); 
00371 
00372   tmpPosX = (beamPosX + beamSizeX*gaussX ) *mm;
00373   tmpPosY = (beamPosY + beamSizeY*gaussY ) *mm;
00374   tmpPosZ = (beamPosZ + beamSizeZ*gaussZ ) *mm;
00375   tmpT = (beamSizeZ * gaussT ) * mm/c_light +beamTime;
00376 
00377   G4LorentzVector tmpv(tmpPosX,tmpPosY,tmpPosZ,tmpT*c_light/mm);
00378 
00379   HepMC::GenEvent::vertex_const_iterator vitr0= hepmcevt->vertices_begin();
00380   G4LorentzVector xvtx0 = (*vitr0)-> position();
00381   pmPosX = xvtx0.x()*mm + tmpPosX;
00382   pmPosY = xvtx0.y()*mm + tmpPosY;
00383   pmPosZ = xvtx0.z()*mm + tmpPosZ;
00384   pmTime = xvtx0.t()*mm/c_light + tmpT;
00385 
00386   if(m_logLevel<=4)
00387   {
00388     G4cout<<G4endl;
00389     G4cout<<xvtx0.x()<<" "<<xvtx0.y()<<" "<<xvtx0.z()<<" "
00390           <<beamSizeX*gaussX<<" "
00391           <<beamSizeY*gaussY<<" "
00392           <<beamSizeZ*gaussZ<<" "<<G4endl;
00393     G4cout<<xvtx0.t()* mm/c_light<<" "
00394           <<beamSizeZ * gaussT/sqrt(2)*mm/c_light<<" "
00395           <<beamTime<<G4endl;
00396   }
00397   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00398               vitr != hepmcevt->vertices_end(); ++vitr )
00399   {
00400     G4LorentzVector xvtx= (*vitr)-> position();
00401     (*vitr)->set_position(xvtx+tmpv);
00402   }
00403   m_G4Svc->SetBeamTime(pmTime);
00404 
00405   G4PrimaryVertex* g4vtx= new G4PrimaryVertex(pmPosX,pmPosY,pmPosZ,pmTime);
00406 
00407   // put initial particles to the vertex
00408   for( size_t ii=0; ii<HPlist.size(); ii++ )
00409   {
00410     if( HPlist[ii]->GetISTHEP() > 0 ) // ISTHEP of daughters had been 
00411                                       // set to negative
00412     {
00413       G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle();
00414       g4vtx->SetPrimary( initialParticle );
00415     }
00416   }
00417   
00418   //clear G4HepMCInterface
00419   for(size_t iii=0;iii<HPlist.size();iii++)
00420   { delete HPlist[iii]; }
00421   HPlist.clear();
00422 
00423   g4event->AddPrimaryVertex(g4vtx);
00424 }

00079 {
00080   G4int i=0;
00081   for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00082       vitr != hepmcevt->vertices_end(); ++vitr ) { // loop for vertex ...
00083     G4cout<<G4endl<<"vertex:"<<i<<"  barcode:"<<(*vitr)->barcode()<<" "; 
00084     i++;
00085     //if((*vitr)->mother())
00086     //  G4cout<<"mother particle: "<<(*vitr)->mother()-> pdg_id()<<"  ";
00087     //if((*vitr)->secondMother())
00088     //  G4cout<<" second mother: "<<(*vitr)->secondMother()-> pdg_id()<<G4endl;
00089     G4LorentzVector xvtx= (*vitr)-> position();
00090     G4cout<<"x:"<<xvtx.x()<<" y:"<<xvtx.y()<<" z:"<<xvtx.z()
00091           <<" t:"<<xvtx.t()*mm/c_light<<G4endl;
00092     
00093     G4int j=0;
00094     for (HepMC::GenVertex::particle_iterator 
00095            pitr= (*vitr)->particles_begin(HepMC::children);
00096          pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00097       G4int pdgcode= (*pitr)-> pdg_id();
00098       G4LorentzVector p= (*pitr)-> momentum();
00099       G4cout<<"particle:"<<j<<" pdgcode:"<<pdgcode<<" "; 
00100       G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" ";
00101       if((*pitr)->end_vertex()) 
00102         G4cout<<"endVtx:"<<(*pitr)->end_vertex()->barcode()<<" ";
00103       G4cout<<"status:"<<(*pitr)->status()<<G4endl;
00104       j++;
00105     }
00106   }
00107 }

00110 {
00111   G4cout<<G4endl;
00112   G4cout<<"particles sent to Geant4: "<<G4endl;
00113   for( size_t ii=0; ii<HPlist.size(); ii++ )
00114   {  
00115     G4int pdgcode = HPlist[ii]->GetTheParticle()->GetPDGcode(); 
00116     G4cout<<pdgcode<<" ";
00117   }
00118   G4cout<<G4endl;
00119 }

00056 {m_logLevel = level;}

00056 {m_logLevel = level;}


Public Member Functions
	BesHepMCInterface ()
	BesHepMCInterface ()
void	Boost (HepMC::GenEvent *hepmcevt)
void	Boost (HepMC::GenEvent *hepmcevt)
G4int	CheckType (const HepMC::GenEvent *hepmcevt)
G4int	CheckType (const HepMC::GenEvent *hepmcevt)
HepMC::GenEvent *	GenerateHepMCEvent ()
HepMC::GenEvent *	GenerateHepMCEvent ()
virtual void	GeneratePrimaryVertex (G4Event *anEvent)
virtual void	GeneratePrimaryVertex (G4Event *anEvent)
HepMC::GenEvent *	GetHepMCGenEvent () const
HepMC::GenEvent *	GetHepMCGenEvent () const
G4int	GetLogLevel ()
G4int	GetLogLevel ()
void	HepMC2G4 (HepMC::GenEvent hepmcevt, G4Event g4event)
void	HepMC2G4 (HepMC::GenEvent hepmcevt, G4Event g4event)
void	Print (const HepMC::GenEvent *hepmcevt)
void	Print (const HepMC::GenEvent *hepmcevt)
void	PrintHPlist ()
void	PrintHPlist ()
void	SetLogLevel (G4int level)
void	SetLogLevel (G4int level)
	~BesHepMCInterface ()
	~BesHepMCInterface ()
Public Attributes
std::vector< G4HepMCParticle * >	HPlist
std::vector< G4HepMCParticle * >	HPlist
Protected Attributes
HepMC::GenEvent *	hepmcEvent
HepMC::GenEvent *	hepmcEvent
G4int	m_logLevel
Private Attributes
IDataProviderSvc *	p_evtSvc
IDataProviderSvc *	p_evtSvc

BesHepMCInterface Class Reference

Public Member Functions

Public Attributes

Protected Attributes

Private Attributes

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation