00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #ifndef WIN32 // Temporarly disabled on Windows, until CLHEP
00034
00035
00036 #include "G4Svc/G4HepMCInterface.h"
00037
00038 #include "globals.hh"
00039 #include "G4LorentzVector.hh"
00040 #include "G4Event.hh"
00041 #include "G4PrimaryParticle.hh"
00042 #include "G4PrimaryVertex.hh"
00043 #include "Randomize.hh"
00044
00045 #include "G4Svc/IG4Svc.h"
00046 #include "G4Svc/G4Svc.h"
00047 #include "GaudiKernel/IDataProviderSvc.h"
00048 #include "GaudiKernel/AlgFactory.h"
00049 #include "GaudiKernel/MsgStream.h"
00050 #include "GaudiKernel/SvcFactory.h"
00051 #include "GaudiKernel/ISvcLocator.h"
00052 #include "GaudiKernel/SmartDataPtr.h"
00053 #include "GaudiKernel/Bootstrap.h"
00055 G4HepMCInterface::G4HepMCInterface()
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 }
00071
00073 G4HepMCInterface::~G4HepMCInterface()
00075 {
00076 delete hepmcEvent;
00077 }
00078 void G4HepMCInterface::Print(const HepMC::GenEvent* hepmcevt)
00079 {
00080 G4int i=0;
00081 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00082 vitr != hepmcevt->vertices_end(); ++vitr ) {
00083 G4cout<<G4endl<<"vertex:"<<i<<" barcode:"<<(*vitr)->barcode()<<" ";
00084 i++;
00085
00086
00087
00088
00089 G4LorentzVector xvtx;
00090 xvtx.setX( (*vitr)-> position().x());
00091 xvtx.setY( (*vitr)-> position().y());
00092 xvtx.setZ( (*vitr)-> position().z());
00093 xvtx.setT( (*vitr)-> position().t());
00094 G4cout<<"x:"<<xvtx.x()<<" y:"<<xvtx.y()<<" z:"<<xvtx.z()
00095 <<" t:"<<xvtx.t()*mm/c_light<<G4endl;
00096
00097 G4int j=0;
00098 for (HepMC::GenVertex::particle_iterator
00099 pitr= (*vitr)->particles_begin(HepMC::children);
00100 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00101 G4int pdgcode= (*pitr)-> pdg_id();
00102 G4LorentzVector p;
00103 p.setX( (*pitr)-> momentum().x());
00104 p.setY( (*pitr)-> momentum().y());
00105 p.setZ( (*pitr)-> momentum().z());
00106 p.setT( (*pitr)-> momentum().t());
00107 G4cout<<"particle:"<<j<<" pdgcode:"<<pdgcode<<" ";
00108 G4cout<<"p:"<<p.x()<<" "<<p.y()<<" "<<p.z()<<" ";
00109 if((*pitr)->end_vertex())
00110 G4cout<<"endVtx:"<<(*pitr)->end_vertex()->barcode()<<" ";
00111 G4cout<<"status:"<<(*pitr)->status()<<G4endl;
00112 j++;
00113 }
00114 }
00115 }
00116
00117 void G4HepMCInterface::PrintHPlist()
00118 {
00119 G4cout<<G4endl;
00120 G4cout<<"particles sent to Geant4: "<<G4endl;
00121 for( size_t ii=0; ii<HPlist.size(); ii++ )
00122 {
00123 G4int pdgcode = HPlist[ii]->GetTheParticle()->GetPDGcode();
00124 G4cout<<pdgcode<<" ";
00125 }
00126 G4cout<<G4endl;
00127 }
00128
00129 G4int G4HepMCInterface::CheckType(const HepMC::GenEvent* hepmcevt)
00130 {
00131 G4int flag =0;
00132 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00133 vitr != hepmcevt->vertices_end(); ++vitr ) {
00134 for (HepMC::GenVertex::particle_iterator
00135 pitr= (*vitr)->particles_begin(HepMC::children);
00136 pitr != (*vitr)->particles_end(HepMC::children); ++pitr) {
00137 if((*pitr)->end_vertex())
00138 {flag = 1; break;}
00139 }
00140 if(flag) break;
00141 }
00142 if(m_logLevel <= 4)
00143 {
00144 if(flag == 0)
00145 G4cout<<G4endl<<"generator is GENBES!"<<G4endl;
00146 else
00147 G4cout<<G4endl<<"generator is not GENBES!"<<G4endl;
00148 }
00149 return flag;
00150 }
00151
00152
00154 void G4HepMCInterface::Boost(HepMC::GenEvent* hepmcevt)
00156 {
00157
00158
00159
00160
00161
00162
00163
00164
00165
00166
00167 G4LorentzVector ptot=0;
00168 double E_cms,p_Mag,e_Mass2,M_cms,theta;
00169 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00170 vitr != hepmcevt->vertices_end(); ++vitr ) {
00171
00172 for(HepMC::GenVertex::particle_iterator
00173 vpitr = (*vitr)->particles_begin(HepMC::children);
00174 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
00175
00176 if((*vpitr)->status() !=1)continue;
00177
00178 G4LorentzVector p;
00179 p.setX( (*vpitr)-> momentum().x());
00180 p.setY( (*vpitr)-> momentum().y());
00181 p.setZ( (*vpitr)-> momentum().z());
00182 p.setT( (*vpitr)-> momentum().t());
00183 ptot = p + ptot;
00184 }
00185 }
00186 E_cms=ptot.e()*GeV;
00187
00188
00189 ISvcLocator* svcLocator = Gaudi::svcLocator();
00190 IG4Svc* tmpSvc;
00191 G4Svc* m_G4Svc;
00192 StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
00193 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00194
00195
00196 G4ThreeVector v=0;
00197
00198 theta=11*mrad;
00199
00200
00201 M_cms=E_cms;
00202 e_Mass2=electron_mass_c2*electron_mass_c2;
00203 p_Mag=sqrt((E_cms*E_cms-4*e_Mass2)/(2*(1-cos(pi-2*theta))));
00204 G4ThreeVector p_Positron(p_Mag*sin(theta),0,p_Mag*cos(theta));
00205 G4ThreeVector p_Electron(p_Mag*sin(pi-theta),0,p_Mag*cos(pi-theta));
00206 G4ThreeVector vee=p_Positron+p_Electron;
00207 G4ThreeVector beamshift(m_G4Svc->GetBeamShiftPx(),m_G4Svc->GetBeamShiftPy(),m_G4Svc->GetBeamShiftPz());
00208 if(m_G4Svc-> GetSetBeamShift()==true && fabs(M_cms-3686)<50.0) {vee = beamshift;}
00209
00210 double pee=vee.r();
00211 M_cms = sqrt(M_cms*M_cms + pee*pee);
00212 v=vee/M_cms;
00213
00214
00215 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00216 vitr != hepmcevt->vertices_end(); ++vitr ) {
00217
00218
00219 G4LorentzVector xvtx;
00220 xvtx.setX( (*vitr)-> position().x());
00221 xvtx.setY( (*vitr)-> position().y());
00222 xvtx.setZ( (*vitr)-> position().z());
00223 xvtx.setT( (*vitr)-> position().t());
00224 xvtx.boost(v);
00225 (*vitr)->set_position(xvtx);
00226
00227
00228
00229 for (HepMC::GenVertex::particle_iterator
00230 vpitr = (*vitr)->particles_begin(HepMC::children);
00231 vpitr != (*vitr)->particles_end(HepMC::children); ++vpitr) {
00232 G4LorentzVector p;
00233 p.setX( (*vpitr)-> momentum().x());
00234 p.setY( (*vpitr)-> momentum().y());
00235 p.setZ( (*vpitr)-> momentum().z());
00236 p.setT( (*vpitr)-> momentum().t());
00237 p=p.boost(v);
00238 (*vpitr)->set_momentum(p);
00239 }
00240 }
00241 }
00242
00244 void G4HepMCInterface::HepMC2G4(HepMC::GenEvent* hepmcevt,
00245 G4Event* g4event)
00247 {
00248
00249 if(m_logLevel <=4 )
00250 Print(hepmcevt);
00251
00252
00253 ISvcLocator* svcLocator = Gaudi::svcLocator();
00254 IG4Svc* tmpSvc;
00255 G4Svc* m_G4Svc;
00256 StatusCode sc=svcLocator->service("G4Svc", tmpSvc);
00257 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00258
00259
00260
00261 if(m_G4Svc->GetBoostLab()== true)
00262 Boost(hepmcevt);
00263
00264
00265 if(m_logLevel <=2 && m_G4Svc->GetBoostLab()== true)
00266 Print(hepmcevt);
00267
00268
00269 G4int flag = CheckType(hepmcevt);
00270 if(flag)
00271 {
00272 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00273 vitr != hepmcevt->vertices_end(); ++vitr )
00274 {
00275 G4int vtxFlag=0;
00276 G4int pOut = (*vitr)->particles_out_size();
00277 HepMC::GenVertex::particle_iterator pitr = (*vitr)->particles_begin(HepMC::children);
00278 G4int pdgcode= (*pitr)-> pdg_id();
00279 if(pOut == 1 && abs(pdgcode) == 11)
00280 vtxFlag=1;
00281 G4int barcodeVtx = (*vitr)->barcode();
00282
00283 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
00284 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
00285 {
00286 G4int pdgcode = (*pitr)->pdg_id();
00287 if(vtxFlag==0)
00288 {
00289 if(pdgcode==-22) pdgcode=22;
00290 G4LorentzVector p;
00291 p.setX( (*pitr)-> momentum().x());
00292 p.setY( (*pitr)-> momentum().y());
00293 p.setZ( (*pitr)-> momentum().z());
00294 p.setT( (*pitr)-> momentum().t());
00295 G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
00296
00297
00298 G4int ISTHEP = (*pitr)->status();
00299
00300 G4int barcodeEndVtx = 0;
00301 if((*pitr)->end_vertex())
00302 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00303
00304 G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
00305 HPlist.push_back(hepmcParticle);
00306
00307 for( size_t ii=0; ii<HPlist.size(); ii++ )
00308 {
00309 if(barcodeVtx == HPlist[ii]->GetBarcodeEndVtx())
00310 {
00311 HPlist[ii]->GetTheParticle()->SetDaughter(particle);
00312 hepmcParticle->Done();
00313 break;
00314 }
00315 }
00316 }
00317 }
00318 }
00319 }
00320 else
00321 {
00322 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00323 vitr != hepmcevt->vertices_end(); ++vitr )
00324 {
00325 for (HepMC::GenVertex::particle_iterator pitr= (*vitr)->particles_begin(HepMC::children);
00326 pitr != (*vitr)->particles_end(HepMC::children); ++pitr)
00327 {
00328 G4int ISTHEP = (*pitr)->status();
00329 G4LorentzVector p;
00330 p.setX( (*pitr)-> momentum().x());
00331 p.setY( (*pitr)-> momentum().y());
00332 p.setZ( (*pitr)-> momentum().z());
00333 p.setT( (*pitr)-> momentum().t());
00334
00335 G4int pdgcode = (*pitr)->pdg_id();
00336 G4int barcodeEndVtx = 0;
00337 if((*pitr)->end_vertex())
00338 barcodeEndVtx = (*pitr)->end_vertex()->barcode();
00339 G4PrimaryParticle* particle = new G4PrimaryParticle(pdgcode, p.x()*GeV, p.y()*GeV, p.z()*GeV);
00340
00341
00342
00343 G4HepMCParticle* hepmcParticle = new G4HepMCParticle(particle,ISTHEP,barcodeEndVtx);
00344 HPlist.push_back(hepmcParticle);
00345 if(ISTHEP>1)
00346 hepmcParticle->Done();
00347 }
00348 }
00349 }
00350
00351
00352 if(m_logLevel<=4)
00353 PrintHPlist();
00354
00355
00356 G4double pmPosX,pmPosY,pmPosZ,pmTime;
00357 G4double tmpPosX,tmpPosY,tmpPosZ,tmpT;
00358 G4double beamPosX,beamPosY,beamPosZ,beamSizeX,beamSizeY,beamSizeZ;
00359
00360 if(m_RealizationSvc->UseDBFlag() == false) {
00361 beamPosX = m_G4Svc->GetBeamPosX();
00362 beamPosY = m_G4Svc->GetBeamPosY();
00363 beamPosZ = m_G4Svc->GetBeamPosZ();
00364
00365 beamSizeX = m_G4Svc->GetBeamSizeX();
00366 beamSizeY = m_G4Svc->GetBeamSizeY();
00367 beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
00368 }
00369 if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == true) {
00370 beamPosX = m_RealizationSvc->getBunchPosX();
00371 beamPosY = m_RealizationSvc->getBunchPosY();
00372 beamPosZ = m_RealizationSvc->getBunchPosZ();
00373
00374 beamSizeX = m_RealizationSvc->getBunchSizeX();
00375 beamSizeY = m_RealizationSvc->getBunchSizeY();
00376 beamSizeZ = m_RealizationSvc->getBunchSizeZ();
00377 }
00378 if(m_RealizationSvc->UseDBFlag() == true && m_RealizationSvc->ifReadBunch() == false) {
00379 beamPosX = m_G4Svc->GetBeamPosX();
00380 beamPosY = m_G4Svc->GetBeamPosY();
00381 beamPosZ = m_G4Svc->GetBeamPosZ();
00382
00383 beamSizeX = m_G4Svc->GetBeamSizeX();
00384 beamSizeY = m_G4Svc->GetBeamSizeY();
00385 beamSizeZ = m_G4Svc->GetBeamSizeZ()/sqrt(2);
00386 }
00387
00388 G4double gaussX = G4RandGauss::shoot();
00389 G4double gaussY = G4RandGauss::shoot();
00390 G4double gaussZ = G4RandGauss::shoot();
00391 G4double gaussT = G4RandGauss::shoot();
00392
00393 G4double beamStartTime = m_G4Svc->GetBeamStartTime() * ns;
00394 G4double beamDeltaTime = m_G4Svc->GetBeamDeltaTime() * ns;
00395 G4double nBunch = m_G4Svc->GetNBunch();
00396 G4double bunchTimeSigma = m_G4Svc->GetBunchTimeSigma() * ns;
00397
00398 G4double ran=G4UniformRand();
00399 G4double beamTime = bunchTimeSigma*G4RandGauss::shoot() + beamStartTime + beamDeltaTime*int(ran*nBunch);
00400
00401 tmpPosX = (beamPosX + beamSizeX*gaussX ) *mm;
00402 tmpPosY = (beamPosY + beamSizeY*gaussY ) *mm;
00403 tmpPosZ = (beamPosZ + beamSizeZ*gaussZ ) *mm;
00404 tmpT = (beamSizeZ * gaussT ) * mm/c_light +beamTime;
00405
00406 G4LorentzVector tmpv(tmpPosX,tmpPosY,tmpPosZ,tmpT*c_light/mm);
00407
00408 HepMC::GenEvent::vertex_const_iterator vitr0= hepmcevt->vertices_begin();
00409 G4LorentzVector xvtx0 ;
00410 xvtx0.setX( (*vitr0)-> position().x());
00411 xvtx0.setY( (*vitr0)-> position().y());
00412 xvtx0.setZ( (*vitr0)-> position().z());
00413 xvtx0.setT( (*vitr0)-> position().t());
00414 pmPosX = xvtx0.x()*mm + tmpPosX;
00415 pmPosY = xvtx0.y()*mm + tmpPosY;
00416 pmPosZ = xvtx0.z()*mm + tmpPosZ;
00417 pmTime = xvtx0.t()*mm/c_light + tmpT;
00418
00419 if(m_logLevel<=4)
00420 {
00421 G4cout<<G4endl;
00422 G4cout<<xvtx0.x()<<" "<<xvtx0.y()<<" "<<xvtx0.z()<<" "
00423 <<beamSizeX*gaussX<<" "
00424 <<beamSizeY*gaussY<<" "
00425 <<beamSizeZ*gaussZ<<" "<<G4endl;
00426 G4cout<<xvtx0.t()* mm/c_light<<" "
00427 <<beamSizeZ * gaussT/sqrt(2)*mm/c_light<<" "
00428 <<beamTime<<G4endl;
00429 }
00430 for(HepMC::GenEvent::vertex_const_iterator vitr= hepmcevt->vertices_begin();
00431 vitr != hepmcevt->vertices_end(); ++vitr )
00432 {
00433 G4LorentzVector xvtx;
00434 xvtx.setX((*vitr)-> position().x());
00435 xvtx.setY((*vitr)-> position().y());
00436 xvtx.setZ((*vitr)-> position().z());
00437 xvtx.setT((*vitr)-> position().t());
00438 (*vitr)->set_position(xvtx+tmpv);
00439 }
00440 m_G4Svc->SetBeamTime(pmTime);
00441
00442 G4PrimaryVertex* g4vtx= new G4PrimaryVertex(pmPosX,pmPosY,pmPosZ,pmTime);
00443
00444
00445 for( size_t ii=0; ii<HPlist.size(); ii++ )
00446 {
00447 if( HPlist[ii]->GetISTHEP() > 0 )
00448
00449 {
00450 G4PrimaryParticle* initialParticle = HPlist[ii]->GetTheParticle();
00451 g4vtx->SetPrimary( initialParticle );
00452 }
00453 }
00454
00455
00456 for(size_t iii=0;iii<HPlist.size();iii++)
00457 { delete HPlist[iii]; }
00458 HPlist.clear();
00459
00460 g4event->AddPrimaryVertex(g4vtx);
00461 }
00462
00463
00465 HepMC::GenEvent* G4HepMCInterface::GenerateHepMCEvent()
00467 {
00468 HepMC::GenEvent* aevent= new HepMC::GenEvent();
00469 return aevent;
00470 }
00471
00473 void G4HepMCInterface::GeneratePrimaryVertex(G4Event* anEvent)
00475 {
00476
00477 delete hepmcEvent;
00478
00479
00480 hepmcEvent= GenerateHepMCEvent();
00481 if(! hepmcEvent) {
00482 G4Exception("G4HepMCInterface","GeneratePrimaryVertex",RunMustBeAborted,
00483 "HepMCInterface: no generated particles. run terminated..." );
00484 return;
00485 }
00486 HepMC2G4(hepmcEvent, anEvent);
00487 }
00488
00489 #endif