#include <G4HepMCInterface.h>
Inheritance diagram for G4HepMCInterface:
Public Member Functions | |
void | Boost (HepMC::GenEvent *hepmcevt) |
void | Boost (HepMC::GenEvent *hepmcevt) |
G4int | CheckType (const HepMC::GenEvent *hepmcevt) |
G4int | CheckType (const HepMC::GenEvent *hepmcevt) |
G4HepMCInterface () | |
G4HepMCInterface () | |
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) |
virtual | ~G4HepMCInterface () |
virtual | ~G4HepMCInterface () |
Public Attributes | |
std::vector< G4HepMCParticle * > | HPlist |
std::vector< G4HepMCParticle * > | HPlist |
Protected Member Functions | |
virtual HepMC::GenEvent * | GenerateHepMCEvent () |
virtual HepMC::GenEvent * | GenerateHepMCEvent () |
Protected Attributes | |
HepMC::GenEvent * | hepmcEvent |
HepMC::GenEvent * | hepmcEvent |
G4int | m_logLevel |
Private Attributes | |
RealizationSvc * | m_RealizationSvc |
RealizationSvc * | m_RealizationSvc |
|
00056 : hepmcEvent(0) 00058 { 00059 00060 IRealizationSvc *tmpReal; 00061 ISvcLocator* svcLocator = Gaudi::svcLocator(); 00062 StatusCode status = svcLocator->service("RealizationSvc",tmpReal); 00063 if (!status.isSuccess()) 00064 { 00065 std::cout << " Could not initialize Realization Service in G4HepMCInterface" << std::endl; 00066 } else { 00067 m_RealizationSvc=dynamic_cast<RealizationSvc*>(tmpReal); 00068 } 00069 00070 }
|
|
00075 {
00076 delete hepmcEvent;
00077 }
|
|
|
|
|
|
|
|
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 }
|
|
Reimplemented in BesHepMCInterface, and BesHepMCInterface. |
|
Reimplemented in BesHepMCInterface, and BesHepMCInterface. 00430 { 00431 HepMC::GenEvent* aevent= new HepMC::GenEvent(); 00432 return aevent; 00433 }
|
|
|
|
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;}
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|