BDSIM: /scratch0/jsnuveri/BDSIM/BDSIMgit/bdsim/src/BDSPhysicsList.cc Source File

00001 //   BDSIM, (C) 2001-2007 
00002 //    
00003 //   version 0.4 
00004 //   last modified : 10 Sept 2007 by malton@pp.rhul.ac.uk
00005 //  
00006 
00007 //
00008 //    Physics lists
00009 //
00010 
00011 #include <iomanip>   
00012 
00013 #include "BDSDebug.hh"
00014 #include "BDSExecOptions.hh"
00015 #include "BDSGlobalConstants.hh" 
00016 #include "BDSPhysicsList.hh"
00017 
00018 #include "globals.hh"
00019 #include "G4ParticleDefinition.hh"
00020 #include "G4ProcessManager.hh"
00021 #include "G4ProcessVector.hh"
00022 #include "G4ParticleTypes.hh"
00023 #include "G4ParticleTable.hh"
00024 #include "G4FastSimulationManagerProcess.hh"
00025 
00026 #include "G4Version.hh"
00027 
00028 #include "G4Material.hh"
00029 #include "G4MaterialTable.hh"
00030 #include "G4ios.hh"
00031 #include "QGSP_BERT.hh"
00032 #include "QGSP_BERT_HP.hh"
00033 #if G4VERSION_NUMBER < 1000
00034 #include "HadronPhysicsQGSP_BERT.hh"
00035 #include "HadronPhysicsQGSP_BERT_HP.hh"
00036 #include "HadronPhysicsFTFP_BERT.hh"
00037 #include "HadronPhysicsLHEP.hh"
00038 #else
00039 #include "G4HadronPhysicsQGSP_BERT.hh"
00040 #include "G4HadronPhysicsQGSP_BERT_HP.hh"
00041 #include "G4HadronPhysicsQGSP_BIC_HP.hh"
00042 #include "G4HadronPhysicsFTFP_BERT.hh"
00043 #endif
00044 #include "G4Decay.hh"
00045 #include "G4eeToHadrons.hh"
00046 
00047 #include "G4EmStandardPhysics.hh"
00048 #include "G4EmLivermorePhysics.hh"
00049 #include "G4EmPenelopePhysics.hh"
00050 
00051 // physics processes
00052 
00053 //
00054 // EM
00055 //
00056 
00057 // gamma
00058 #include "G4ComptonScattering.hh"
00059 #include "G4GammaConversion.hh"
00060 #include "G4GammaConversionToMuons.hh"
00061 #include "G4PhotoElectricEffect.hh"
00062 #include "BDSXSBias.hh" //added by L.D. 16/11/09
00063 #include "BDSGammaConversion_LPB.hh" //added by M.D. Salt, R.B. Appleby, 15/10/09
00064 
00065 // charged particles
00066 #include "G4eMultipleScattering.hh"
00067 #include "G4MuMultipleScattering.hh"
00068 #include "G4hMultipleScattering.hh"
00069 
00070 //Optical processes
00071 #include "G4Cerenkov.hh"
00072 #include "G4Scintillation.hh"
00073 #include "G4OpAbsorption.hh"
00074 #include "G4OpRayleigh.hh"
00075 #include "G4OpMieHG.hh"
00076 #include "G4OpBoundaryProcess.hh"
00077 #include "G4CoulombScattering.hh"
00078 
00079 // e
00080 #include "G4eIonisation.hh"
00081 #include "G4eBremsstrahlung.hh"
00082 #include "G4eplusAnnihilation.hh"
00083 #include "BDSeBremsstrahlung_LPB.hh" //added by M.D. Salt, R.B.Appleby,15/10/09
00084 
00085 // mu
00086 #include "G4MuIonisation.hh"
00087 #include "G4MuBremsstrahlung.hh"
00088 #include "G4MuPairProduction.hh"
00089 #if G4VERSION_NUMBER < 953
00090 #include "G4MuNuclearInteraction.hh"
00091 #else 
00092 #include "G4MuonVDNuclearModel.hh"
00093 #endif
00094 
00095 //ions
00096 #include "G4hIonisation.hh"
00097 #include "G4ionIonisation.hh"
00098 
00099 
00100 //Decay of unstable particles
00101 #include "G4Decay.hh"
00102 
00103 //
00104 // Low EM
00105 //
00106 
00107 //gamma
00108 #include "G4RayleighScattering.hh"
00109 #include "G4LivermoreRayleighModel.hh"
00110 #include "G4PhotoElectricEffect.hh"
00111 #include "G4LivermorePhotoElectricModel.hh"
00112 #include "G4ComptonScattering.hh"
00113 #include "G4LivermoreComptonModel.hh"
00114 #include "G4GammaConversion.hh"
00115 #include "G4LivermoreGammaConversionModel.hh"
00116 
00117 //e
00118 #include "G4eIonisation.hh"
00119 #include "G4LivermoreIonisationModel.hh"
00120 #include "G4UniversalFluctuation.hh"
00121 
00122 #include "G4eBremsstrahlung.hh"
00123 #include "G4LivermoreBremsstrahlungModel.hh"
00124 
00125 #include "G4AnnihiToMuPair.hh"
00126 
00127 //ions
00128 #include "BDSLaserCompton.hh"
00129 #include "BDSSynchrotronRadiation.hh"
00130 #include "BDSContinuousSR.hh"
00131 #include "G4StepLimiter.hh"
00132 #include "G4UserSpecialCuts.hh"
00133 
00134 //
00135 // Hadronic
00136 //
00137 #include "G4TheoFSGenerator.hh"
00138 #include "G4GeneratorPrecompoundInterface.hh"
00139 #include "G4QGSModel.hh"
00140 #include "G4GammaParticipants.hh"
00141 #include "G4QGSMFragmentation.hh"
00142 #include "G4ExcitedStringDecay.hh"
00143 
00144 #if G4VERSION_NUMBER < 1000
00145 #include "G4GammaNuclearReaction.hh"
00146 #include "G4ElectroNuclearReaction.hh"
00147 #endif
00148 #include "G4PhotoNuclearProcess.hh"
00149 #include "G4ElectronNuclearProcess.hh"
00150 #include "G4PositronNuclearProcess.hh"
00151 
00152 //Planck scattering
00153 #include "BDSPlanckScatterBuilder.hh"
00154 
00155 //
00156 // particle definition
00157 //
00158 
00159 // Bosons
00160 #include "G4ChargedGeantino.hh"
00161 #include "G4Geantino.hh"
00162 #include "G4Gamma.hh"
00163 #include "G4OpticalPhoton.hh"
00164 
00165 // leptons
00166 #include "G4MuonPlus.hh"
00167 #include "G4MuonMinus.hh"
00168 #include "G4NeutrinoMu.hh"
00169 #include "G4AntiNeutrinoMu.hh"
00170 
00171 #include "G4TauPlus.hh"
00172 #include "G4TauMinus.hh"
00173 #include "G4NeutrinoTau.hh"
00174 #include "G4AntiNeutrinoTau.hh"
00175 
00176 #include "G4Electron.hh"
00177 #include "G4Positron.hh"
00178 #include "G4NeutrinoE.hh"
00179 #include "G4AntiNeutrinoE.hh"
00180 
00181 // Hadrons
00182 #include "G4MesonConstructor.hh"
00183 #include "G4BaryonConstructor.hh"
00184 #include "G4IonConstructor.hh"
00185 
00186 //ShortLived
00187 #include "G4ShortLivedConstructor.hh"
00188 
00189 
00190 
00191 
00192 BDSPhysicsList::BDSPhysicsList():  G4VUserPhysicsList()
00193 {
00194   theCerenkovProcess           = NULL;
00195   theScintillationProcess      = NULL;
00196   theAbsorptionProcess         = NULL;
00197   theRayleighScatteringProcess = NULL;
00198   theMieHGScatteringProcess    = NULL;
00199   theBoundaryProcess           = NULL;
00200 
00201   theReferenceHadronicPhysList = NULL;
00202   theReferenceEmPhysList       = NULL;
00203   theBDSIMPhysList             = NULL;
00204   theHadPhysList1              = NULL;
00205   theHadPhysList2              = NULL;
00206 
00207   verbose = BDSExecOptions::Instance()->GetVerbose();
00208 
00209   // construct particles
00210 
00211   //defaultCutValue = 0.7*CLHEP::mm;  
00212   defaultCutValue = BDSGlobalConstants::Instance()->GetDefaultRangeCut()*CLHEP::m;  
00213   SetDefaultCutValue(BDSGlobalConstants::Instance()->GetDefaultRangeCut()*CLHEP::m);
00214 
00215   G4cout  << __METHOD_NAME__ << "Charged Thresholdcut=" 
00216           << BDSGlobalConstants::Instance()->GetThresholdCutCharged()/CLHEP::GeV<<" GeV"<<G4endl;
00217   G4cout  << __METHOD_NAME__ << "Photon Thresholdcut=" 
00218           << BDSGlobalConstants::Instance()->GetThresholdCutPhotons()/CLHEP::GeV<<" GeV"<<G4endl;
00219   G4cout  << __METHOD_NAME__ << "Default range cut=" 
00220           << BDSGlobalConstants::Instance()->GetDefaultRangeCut()/CLHEP::m<<" m"<<G4endl;
00221 
00222   //This is the GEANT4 physics list verbose level.
00223   SetVerboseLevel(1);   
00224 }
00225 
00226 BDSPhysicsList::~BDSPhysicsList()
00227 {
00228   delete theReferenceHadronicPhysList;
00229   delete theReferenceEmPhysList;
00230   delete theBDSIMPhysList;
00231   delete theHadPhysList1;
00232   delete theHadPhysList2;
00233 }
00234 
00235 void BDSPhysicsList::ConstructProcess()
00236 { 
00237 #if BDSDEBUG
00238   G4cout << __METHOD_NAME__ << G4endl;
00239 #endif 
00240 
00241   bool plistFound=false;
00242   //standard physics lists
00243   if(BDSGlobalConstants::Instance()->GetPhysListName() == "QGSP_BERT"){
00244 #if G4VERSION_NUMBER < 1000
00245     theReferenceHadronicPhysList = new HadronPhysicsQGSP_BERT();
00246 #else
00247     theReferenceHadronicPhysList = new G4HadronPhysicsQGSP_BERT();
00248 #endif
00249     theReferenceEmPhysList = new G4EmStandardPhysics();
00250     theReferenceHadronicPhysList->ConstructProcess();
00251     theReferenceEmPhysList->ConstructProcess();
00252     plistFound=true;
00253   } else if(BDSGlobalConstants::Instance()->GetPhysListName() == "QGSP_BERT_HP"){
00254 #if G4VERSION_NUMBER < 1000
00255     theReferenceHadronicPhysList = new HadronPhysicsQGSP_BERT_HP();
00256 #else
00257     theReferenceHadronicPhysList = new G4HadronPhysicsQGSP_BERT_HP();
00258 #endif
00259     theReferenceHadronicPhysList->ConstructProcess();
00260     theReferenceEmPhysList = new G4EmStandardPhysics();
00261     theReferenceEmPhysList->ConstructProcess();
00262     plistFound=true;
00263   } else if (BDSGlobalConstants::Instance()->GetPhysListName() == "QGSP_BERT_HP_muon"){ //Modified standard physics lists
00264 #if G4VERSION_NUMBER < 1000
00265     theReferenceHadronicPhysList = new HadronPhysicsQGSP_BERT_HP();
00266 #else
00267     theReferenceHadronicPhysList = new G4HadronPhysicsQGSP_BERT_HP();
00268 #endif
00269     theReferenceHadronicPhysList->ConstructProcess();
00270     ConstructMuon();
00271     plistFound=true;
00272   } else if (BDSGlobalConstants::Instance()->GetPhysListName() == "QGSP_BERT_muon"){
00273 #if G4VERSION_NUMBER < 1000
00274     theReferenceHadronicPhysList = new HadronPhysicsQGSP_BERT();
00275 #else
00276     theReferenceHadronicPhysList = new G4HadronPhysicsQGSP_BERT();
00277 #endif
00278     theReferenceHadronicPhysList->ConstructProcess();
00279     ConstructMuon();
00280     plistFound=true;
00281   } else if(BDSGlobalConstants::Instance()->GetPhysListName() == "QGSP_BERT_HP_muon_em_low"){
00282 #if G4VERSION_NUMBER < 1000
00283     theReferenceHadronicPhysList = new HadronPhysicsQGSP_BERT_HP();
00284 #else
00285     theReferenceHadronicPhysList = new G4HadronPhysicsQGSP_BERT_HP();
00286 #endif
00287     theReferenceHadronicPhysList->ConstructProcess();
00288     ConstructMuon();
00289     ConstructEM_Low_Energy();
00290     plistFound=true;
00291   } else if(BDSGlobalConstants::Instance()->GetPhysListName() == "livermore"){
00292     G4EmLivermorePhysics* physList = new G4EmLivermorePhysics;
00293     physList->ConstructProcess();
00294     plistFound=true;
00295   }else if(BDSGlobalConstants::Instance()->GetPhysListName() == "penelope"){
00296     G4EmPenelopePhysics* physList = new G4EmPenelopePhysics;
00297     physList->ConstructProcess();
00298     plistFound=true;
00299   }else if(BDSGlobalConstants::Instance()->GetPhysListName() == "G4EmStandard"){
00300     G4EmStandardPhysics* physList = new G4EmStandardPhysics;
00301     physList->ConstructProcess();
00302     plistFound=true;
00303   }
00304 
00305 
00306 
00307   
00308   //===========================================
00309   //Some options
00310   //-------------------------------------------
00311   //Build planck scattering if option is set
00312   if(BDSGlobalConstants::Instance()->GetDoPlanckScattering()){
00313     BDSPlanckScatterBuilder* psbuild = new BDSPlanckScatterBuilder();
00314     psbuild->Build();
00315     // psbuild is leaked. This can't be deleted as its BDSPlanckScatterProcess is registered
00316   }
00317   //A flag to switch on hadronic lead particle biasing
00318   if (BDSGlobalConstants::Instance()->GetUseHadLPB() ){
00319     setenv("SwitchLeadBiasOn","1",1); 
00320   }
00321   //Synchrotron radiation
00322   if(BDSGlobalConstants::Instance()->GetSynchRadOn()) {
00323 #ifdef BDSDEBUG
00324     G4cout << __METHOD_NAME__ << "synch. rad. is turned on" << G4endl;
00325 #endif
00326     ConstructSR();
00327   } else {
00328 #ifdef BDSDEBUG
00329     G4cout << __METHOD_NAME__ << "synch. rad. is turned OFF!" << G4endl;
00330 #endif
00331   }
00332   //Particle decay
00333   if(BDSGlobalConstants::Instance()->GetDecayOn()) ConstructDecay();
00334   //Optical processes
00335   ConstructOptical();
00336   //============================================
00337   //Need to add transportation and step limiter in all cases.
00338   AddTransportation();
00339   theParticleIterator->reset();
00340   while( (*theParticleIterator)() ){
00341     G4ParticleDefinition* particle = theParticleIterator->value();
00342     G4ProcessManager *pmanager = particle->GetProcessManager();
00343     if((particle->GetParticleName()=="gamma")||
00344        (particle->GetParticleName()=="e-")||
00345        (particle->GetParticleName()=="e+")||
00346        (particle->GetParticleName()=="proton")){
00347       particle->SetApplyCutsFlag(true);
00348     }
00349     pmanager->AddProcess(new G4StepLimiter,-1,-1,1);
00350 #ifndef NOUSERSPECIALCUTS
00351     pmanager->AddDiscreteProcess(new G4UserSpecialCuts);
00352 #endif
00353   }
00354   //Always add parameterisation
00355   AddParameterisation();
00356   
00357   if(plistFound) return;
00358   //Search BDSIM physics lists
00359   if (BDSGlobalConstants::Instance()->GetPhysListName() == "standard") return;
00360   // register physics processes here
00361   // standard e+/e-/gamma electromagnetic interactions
00362   if(BDSGlobalConstants::Instance()->GetPhysListName() == "em_standard") 
00363     {
00364       ConstructEM();
00365     }
00366   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "em_single_scatter") 
00367     {
00368       ConstructEMSingleScatter();
00369     }
00370   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "merlin") 
00371     {
00372       ConstructMerlin();
00373     }
00374   
00375   // low energy em processes
00376   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "em_low") 
00377     {
00378       ConstructEM_Low_Energy();
00379     }
00380       
00381   // standard electromagnetic + muon
00382   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "em_muon") 
00383     {
00384       ConstructEM();
00385       ConstructMuon();
00386     }
00387   // standard hadronic - photo-nuclear etc.
00388   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_standard") 
00389     {
00390       ConstructEM();
00391       ConstructHadronic();
00392     }
00393       
00394   // standard electromagnetic + muon + hadronic
00395   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_muon") 
00396     {
00397       ConstructEM();
00398       ConstructMuon();
00399       ConstructHadronic();
00400     }
00401   
00402   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_QGSP_BERT") 
00403     {
00404       ConstructEM();
00405 #if G4VERSION_NUMBER < 1000
00406       theBDSIMPhysList = new HadronPhysicsQGSP_BERT("hadron");
00407 #else
00408       theBDSIMPhysList = new G4HadronPhysicsQGSP_BERT("hadron");
00409 #endif
00410       theBDSIMPhysList->ConstructProcess();
00411     }
00412   
00413   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_QGSP_BERT_muon") 
00414     {
00415       ConstructEM();
00416       ConstructMuon();
00417 #if G4VERSION_NUMBER < 1000
00418       theBDSIMPhysList = new HadronPhysicsQGSP_BERT("hadron");
00419 #else
00420       theBDSIMPhysList = new G4HadronPhysicsQGSP_BERT("hadron");
00421 #endif
00422       theBDSIMPhysList->ConstructProcess();
00423     }
00424   
00425   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_FTFP_BERT"){
00426     ConstructEM();
00427 #if G4VERSION_NUMBER < 1000
00428     theBDSIMPhysList = new HadronPhysicsFTFP_BERT;
00429 #else
00430     theBDSIMPhysList = new G4HadronPhysicsFTFP_BERT;
00431 #endif
00432     theBDSIMPhysList->ConstructProcess();
00433   }
00434   
00435   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_QGSP_BERT_HP_muon"){
00436     ConstructEM();
00437     ConstructMuon();
00438     ConstructHadronic();
00439 #if G4VERSION_NUMBER < 1000
00440     theBDSIMPhysList = new HadronPhysicsQGSP_BERT_HP;
00441 #else
00442     theBDSIMPhysList = new G4HadronPhysicsQGSP_BERT_HP;
00443 #endif
00444     theBDSIMPhysList->ConstructProcess();
00445   }
00446       
00447   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "hadronic_FTFP_BERT_muon"){
00448     G4cout << __METHOD_NAME__ << "Using hadronic_FTFP_BERT_muon" << G4endl;
00449     ConstructEM();
00450     ConstructMuon();
00451 #if G4VERSION_NUMBER < 1000
00452     theBDSIMPhysList = new HadronPhysicsFTFP_BERT;
00453 #else
00454     theBDSIMPhysList = new G4HadronPhysicsFTFP_BERT;
00455 #endif
00456     theBDSIMPhysList->ConstructProcess();
00457   }
00458   // physics list for laser wire - standard em stuff +
00459   // weighted compton scattering from laser wire
00460   else if(BDSGlobalConstants::Instance()->GetPhysListName() == "lw") {
00461     ConstructEM();
00462     ConstructLaserWire();
00463   }
00464   else {
00465     G4cerr<<"WARNING : Unknown physics list "<<BDSGlobalConstants::Instance()->GetPhysListName()<<G4endl;
00466     exit(1);
00467   }
00468 }
00469 
00470 void BDSPhysicsList::ConstructParticle()
00471 {
00472 #ifdef BDSDEBUG
00473   G4cout << __METHOD_NAME__ << G4endl;
00474 #endif
00475   //standard physics lists
00476   if (theReferenceHadronicPhysList || theReferenceEmPhysList) {
00477     if (theReferenceHadronicPhysList) {
00478       theReferenceHadronicPhysList->ConstructParticle();
00479     }
00480     if (theReferenceEmPhysList) {
00481       theReferenceEmPhysList->ConstructParticle();
00482     }
00483   } else {
00484     // pseudo-particles
00485     G4Geantino::GeantinoDefinition();
00486     G4ChargedGeantino::ChargedGeantinoDefinition();
00487     
00488     // gamma
00489     G4Gamma::GammaDefinition();
00490     
00491     // optical photon
00492     G4OpticalPhoton::OpticalPhotonDefinition();
00493     
00494     // leptons
00495     G4Electron::ElectronDefinition();
00496     G4Positron::PositronDefinition();
00497     G4MuonPlus::MuonPlusDefinition();
00498     G4MuonMinus::MuonMinusDefinition();
00499     G4TauPlus::TauPlusDefinition();
00500     G4TauMinus::TauMinusDefinition();
00501 
00502     G4NeutrinoE::NeutrinoEDefinition();
00503     G4AntiNeutrinoE::AntiNeutrinoEDefinition();
00504     G4NeutrinoMu::NeutrinoMuDefinition();
00505     G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
00506     G4NeutrinoTau::NeutrinoTauDefinition();
00507     G4AntiNeutrinoTau::AntiNeutrinoTauDefinition();
00508 
00509     // mesons
00510     G4MesonConstructor mConstructor;
00511     mConstructor.ConstructParticle();
00512     
00513     // baryons
00514     G4BaryonConstructor bConstructor;
00515     bConstructor.ConstructParticle();
00516     
00517     // ions
00518     G4IonConstructor iConstructor;
00519     iConstructor.ConstructParticle();
00520     
00521     //  Construct resonances and quarks
00522     G4ShortLivedConstructor pShortLivedConstructor;
00523     pShortLivedConstructor.ConstructParticle();
00524   }
00525   
00526   // set primary particle definition and kinetic beam parameters other than total energy
00527   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
00528   BDSGlobalConstants::Instance()->SetParticleDefinition(particleTable->
00529                                     FindParticle(BDSGlobalConstants::Instance()->GetParticleName()) );  
00530   
00531   if(!BDSGlobalConstants::Instance()->GetParticleDefinition()) 
00532     {
00533       G4Exception("Particle not found, quitting!", "-1", FatalException, "");
00534       exit(1);
00535     }
00536   
00537   // set kinetic beam parameters other than total energy
00538   BDSGlobalConstants::Instance()->SetBeamMomentum( sqrt(pow(BDSGlobalConstants::Instance()->GetBeamTotalEnergy(),2)-
00539                                     pow(BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGMass(),2)) );
00540   
00541   BDSGlobalConstants::Instance()->SetBeamKineticEnergy(BDSGlobalConstants::Instance()->GetBeamTotalEnergy() - 
00542                                    BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGMass() );
00543 
00544 
00545   BDSGlobalConstants::Instance()->SetParticleMomentum( sqrt(pow(BDSGlobalConstants::Instance()->GetParticleTotalEnergy(),2)-
00546                                     pow(BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGMass(),2)) );
00547   
00548   BDSGlobalConstants::Instance()->SetParticleKineticEnergy(BDSGlobalConstants::Instance()->GetParticleTotalEnergy() - 
00549                                    BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGMass() );
00550   
00551   G4cout << __METHOD_NAME__ << "Beam properties:"<<G4endl;
00552   G4cout << __METHOD_NAME__ << "Particle : " 
00553          << BDSGlobalConstants::Instance()->GetParticleDefinition()->GetParticleName()<<G4endl;
00554   G4cout << __METHOD_NAME__ << "Mass : " 
00555          << BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGMass()/CLHEP::GeV<< " GeV"<<G4endl;
00556   G4cout << __METHOD_NAME__ << "Charge : " 
00557          << BDSGlobalConstants::Instance()->GetParticleDefinition()->GetPDGCharge()<< " e"<<G4endl;
00558   G4cout << __METHOD_NAME__ << "Total Energy : "
00559          << BDSGlobalConstants::Instance()->GetBeamTotalEnergy()/CLHEP::GeV<<" GeV"<<G4endl;
00560   G4cout << __METHOD_NAME__ << "Kinetic Energy : "
00561          << BDSGlobalConstants::Instance()->GetBeamKineticEnergy()/CLHEP::GeV<<" GeV"<<G4endl;
00562   G4cout << __METHOD_NAME__ << "Momentum : "
00563          << BDSGlobalConstants::Instance()->GetBeamMomentum()/CLHEP::GeV<<" GeV"<<G4endl;
00564 }
00565 
00566 #include "G4Region.hh"
00567 #include "G4ProductionCuts.hh"
00568 void BDSPhysicsList::SetCuts()
00569 {
00570   if (verbose){
00571     G4cout << __METHOD_NAME__ << " setting cuts\n";
00572     
00573   }
00574   
00575   SetCutsWithDefault();   
00576 
00577 
00578   
00579     if(BDSGlobalConstants::Instance()->GetProdCutPhotons()>0)
00580       SetCutValue(BDSGlobalConstants::Instance()->GetProdCutPhotons(),G4ProductionCuts::GetIndex("gamma"));
00581   
00582    if(BDSGlobalConstants::Instance()->GetProdCutElectrons()>0)
00583      SetCutValue(BDSGlobalConstants::Instance()->GetProdCutElectrons(),G4ProductionCuts::GetIndex("e-"));
00584   
00585   if(BDSGlobalConstants::Instance()->GetProdCutPositrons()>0)
00586     SetCutValue(BDSGlobalConstants::Instance()->GetProdCutPositrons(),G4ProductionCuts::GetIndex("e+"));
00587   
00588 
00589     
00590   if(1)
00591     DumpCutValuesTable(); 
00592 
00593 }
00594 
00595 
00596 // particular physics process constructors
00597 
00598 void BDSPhysicsList::ConstructEM(){
00599 #ifdef BDSDEBUG
00600   G4cout << __METHOD_NAME__ << G4endl;
00601 #endif
00602   ConstructEMMisc();
00603   ConstructMultipleScattering();
00604 }
00605 
00606 void BDSPhysicsList::ConstructEMSingleScatter(){
00607 #ifdef BDSDEBUG
00608   G4cout << __METHOD_NAME__ << G4endl;
00609 #endif
00610   ConstructEMMisc();
00611   ConstructCoulombScattering();
00612 }
00613 
00614 void BDSPhysicsList::ConstructEMMisc()
00615 {
00616 #ifdef BDSDEBUG
00617   G4cout << __METHOD_NAME__ << G4endl;
00618 #endif
00619   theParticleIterator->reset();
00620   while( (*theParticleIterator)() ){
00621     G4ParticleDefinition* particle = theParticleIterator->value();
00622     G4ProcessManager* pmanager = particle->GetProcessManager();
00623     G4String particleName = particle->GetParticleName();
00624     if (particleName == "gamma") {
00625       // gamma         
00626       pmanager->AddDiscreteProcess(new G4PhotoElectricEffect);
00627       pmanager->AddDiscreteProcess(new G4ComptonScattering);
00628       // if(0){
00629       //        G4GammaConversion* gammaconversion = new G4GammaConversion();
00630       //        gammaconversion->SetLambdaFactor(1/1.0e-20);
00631       //        BDSXSBias* gammaconversion_xsbias = new BDSXSBias();
00632       //        gammaconversion_xsbias->RegisterProcess(gammaconversion);
00633       //        gammaconversion_xsbias->eFactor(1e-20);
00634       //        pmanager->AddDiscreteProcess(gammaconversion_xsbias);
00635         
00636       // } else 
00637       if (BDSGlobalConstants::Instance()->GetUseEMLPB()){ //added by M.D. Salt, R.B. Appleby, 15/10/09
00638           G4GammaConversion* gammaconversion = new G4GammaConversion();
00639           GammaConversion_LPB* gammaconversion_lpb = new GammaConversion_LPB();
00640           gammaconversion_lpb->RegisterProcess(gammaconversion);
00641           pmanager->AddDiscreteProcess(gammaconversion_lpb);
00642       } else {
00643         pmanager->AddDiscreteProcess(new G4GammaConversion);
00644       }
00645     
00646       
00647     } else if (particleName == "e-") {
00648       pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
00649       // if(0){
00650       //        G4eBremsstrahlung* ebremsstrahlung = new G4eBremsstrahlung();
00651       //        ebremsstrahlung->SetLambdaFactor(1/1.0e-20);
00652       //        BDSXSBias* ebremsstrahlung_xsbias = new BDSXSBias();
00653       //        ebremsstrahlung_xsbias->RegisterProcess(ebremsstrahlung);
00654       //        ebremsstrahlung_xsbias->eFactor(1e-20);
00655       //        pmanager->AddDiscreteProcess(ebremsstrahlung_xsbias);     
00656       // }      else 
00657       if(BDSGlobalConstants::Instance()->GetUseEMLPB()){ //added by M.D. Salt, R.B. Appleby, 15/10/09
00658           
00659         G4eBremsstrahlung* ebremsstrahlung = new G4eBremsstrahlung();
00660         eBremsstrahlung_LPB* ebremsstrahlung_lpb = new eBremsstrahlung_LPB();
00661         ebremsstrahlung_lpb->RegisterProcess(ebremsstrahlung);
00662         pmanager->AddProcess(ebremsstrahlung_lpb,     -1,-1,3);
00663       } else {
00664         G4eBremsstrahlung* ebremsstrahlung = new G4eBremsstrahlung();
00665         pmanager->AddProcess(ebremsstrahlung,   -1, 3,3);     
00666       }
00667             
00668       if(BDSGlobalConstants::Instance()->GetTurnOnCerenkov()){
00669         G4Cerenkov* theCerenkovProcess = new G4Cerenkov;
00670         pmanager->AddProcess(theCerenkovProcess);
00671         pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
00672       }
00673       
00674     } else if (particleName == "e+") {
00675       //positron
00676       pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
00677       // if(0){
00678       //        G4eBremsstrahlung* ebremsstrahlung = new G4eBremsstrahlung();
00679       //        ebremsstrahlung->SetLambdaFactor(1/1.0e-20);
00680       //        BDSXSBias* ebremsstrahlung_xsbias = new BDSXSBias();
00681       //        ebremsstrahlung_xsbias->RegisterProcess(ebremsstrahlung);
00682       //        ebremsstrahlung_xsbias->eFactor(1e-20);
00683       //        pmanager->AddDiscreteProcess(ebremsstrahlung_xsbias);      
00684       // } else 
00685       if (BDSGlobalConstants::Instance()->GetUseEMLPB()){
00686         G4eBremsstrahlung* ebremsstrahlung = new G4eBremsstrahlung();
00687         eBremsstrahlung_LPB* ebremsstrahlung_lpb = new eBremsstrahlung_LPB();
00688         ebremsstrahlung_lpb->RegisterProcess(ebremsstrahlung);
00689         pmanager->AddProcess(ebremsstrahlung_lpb,     -1,-1,3);
00690       } else {
00691         pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);
00692       }
00693       pmanager->AddProcess(new G4eplusAnnihilation,  0,-1,4);
00694       if(BDSGlobalConstants::Instance()->GetTurnOnCerenkov()){      
00695         G4Cerenkov* theCerenkovProcess = new G4Cerenkov;
00696         pmanager->AddProcess(theCerenkovProcess);
00697         pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
00698       }
00699     } else if ((!particle->IsShortLived()) &&
00700                (particle->GetPDGCharge() != 0.0) && 
00701                (particle->GetParticleName() != "chargedgeantino")) {
00702       //all others charged particles except geantino
00703       pmanager->AddProcess(new G4hIonisation,       -1, 2,2);
00704            if(BDSGlobalConstants::Instance()->GetTurnOnCerenkov()){
00705         G4Cerenkov* theCerenkovProcess = new G4Cerenkov;
00706         pmanager->AddProcess(theCerenkovProcess);
00707         pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
00708       }
00709     }
00710   }
00711 }
00712 
00713 void BDSPhysicsList::ConstructMultipleScattering(){
00714 #ifdef BDSDEBUG
00715   G4cout << __METHOD_NAME__ << G4endl;
00716 #endif
00717   theParticleIterator->reset();
00718   while( (*theParticleIterator)() ){
00719     G4ParticleDefinition* particle     = theParticleIterator->value();
00720     G4ProcessManager*     pmanager     = particle->GetProcessManager();
00721     G4String              particleName = particle->GetParticleName();
00722     if (particleName == "e-") {
00723       //electron
00724       pmanager->AddProcess(new G4eMultipleScattering,-1, 1,1);
00725 
00726       
00727     } else if (particleName == "e+") {
00728       
00729       pmanager->AddProcess(new G4eMultipleScattering,-1, 1,1);
00730     } else if ((!particle->IsShortLived()) &&
00731                (particle->GetPDGCharge() != 0.0) && 
00732                (particle->GetParticleName() != "chargedgeantino")) {
00733       
00734       //all others charged particles except geantino
00735       pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
00736     }
00737   }
00738 }
00739 
00740 void BDSPhysicsList::ConstructCoulombScattering(){
00741 #ifdef BDSDEBUG
00742   G4cout << __METHOD_NAME__ << G4endl;
00743 #endif
00744   theParticleIterator->reset();
00745   while( (*theParticleIterator)() ){
00746     G4ParticleDefinition* particle = theParticleIterator->value();
00747     G4ProcessManager* pmanager = particle->GetProcessManager();
00748     G4String particleName = particle->GetParticleName(); 
00749     if (particleName == "e-") {
00750       pmanager->AddDiscreteProcess(new G4CoulombScattering);
00751     } else if (particleName == "e+") {
00752       pmanager->AddDiscreteProcess(new G4CoulombScattering);
00753     }else if ((!particle->IsShortLived()) &&
00754               (particle->GetPDGCharge() != 0.0) && 
00755               (particle->GetParticleName() != "chargedgeantino")) {
00756       //all others charged particles except geantino
00757       pmanager->AddDiscreteProcess(new G4CoulombScattering);
00758     } 
00759   }
00760 }
00761 
00762 // particular physics process constructors
00763 void BDSPhysicsList::ConstructMuon()
00764 {
00765 #ifdef BDSDEBUG
00766   G4cout << __METHOD_NAME__ << G4endl;
00767 #endif
00768   theParticleIterator->reset();
00769   while( (*theParticleIterator)() ){
00770     G4ParticleDefinition* particle = theParticleIterator->value();
00771     G4ProcessManager* pmanager = particle->GetProcessManager();
00772     G4String particleName = particle->GetParticleName();
00773 
00774     if (particleName == "gamma") {
00775       // gamma         
00776       G4GammaConversionToMuons* gammaconversiontomuons = new G4GammaConversionToMuons();
00777       BDSXSBias* gammaconversiontomuon_xsbias = new BDSXSBias();
00778       gammaconversiontomuons->SetCrossSecFactor(BDSGlobalConstants::Instance()->GetGammaToMuFe());
00779       gammaconversiontomuon_xsbias->RegisterProcess(gammaconversiontomuons);
00780       gammaconversiontomuon_xsbias->eFactor(BDSGlobalConstants::Instance()->GetGammaToMuFe());
00781       pmanager->AddDiscreteProcess(gammaconversiontomuon_xsbias);
00782 #ifdef BDSDEBUG
00783       G4cout << __METHOD_NAME__ << "GammaToMuFe = " << BDSGlobalConstants::Instance()->GetGammaToMuFe() << G4endl;
00784 #endif
00785     } else if (particleName == "e+") {
00786       //positron
00787       //===========ee to hadrons in development================
00788       G4eeToHadrons* eetohadrons = new G4eeToHadrons();
00789       // BDSXSBias* eetohadrons_xsbias = new BDSXSBias();
00790       // G4cout << "eeToHadronsXSBias = " << BDSGlobalConstants::Instance()->GetEeToHadronsFe() << G4endl;
00791       eetohadrons->SetCrossSecFactor(BDSGlobalConstants::Instance()->GetEeToHadronsFe());
00792       //eetohadrons_xsbias->RegisterProcess(eetohadrons);
00793       //eetohadrons_xsbias->eFactor(BDSGlobalConstants::Instance()->GetEeToHadronsFe());
00794       //pmanager->AddDiscreteProcess(eetohadrons_xsbias);
00795       pmanager->AddDiscreteProcess(eetohadrons);
00796       //-------------------------------------------------------
00797       G4AnnihiToMuPair* annihitomupair = new G4AnnihiToMuPair();
00798       BDSXSBias* annihitomupair_xsbias = new BDSXSBias();
00799       annihitomupair->SetCrossSecFactor(BDSGlobalConstants::Instance()->GetAnnihiToMuFe());
00800       annihitomupair_xsbias->RegisterProcess(annihitomupair);
00801       annihitomupair_xsbias->eFactor(BDSGlobalConstants::Instance()->GetAnnihiToMuFe());
00802       pmanager->AddDiscreteProcess(annihitomupair_xsbias); 
00803 #ifdef BDSDEBUG
00804       G4cout << __METHOD_NAME__ << "AnnihiToMuFe = " << BDSGlobalConstants::Instance()->GetAnnihiToMuFe() << G4endl;
00805 #endif    
00806     } else if( particleName == "mu+" || 
00807                particleName == "mu-"    ) {
00808       //muon  
00809       pmanager->AddProcess(new G4MuMultipleScattering,-1, 1,1);
00810       pmanager->AddProcess(new G4MuIonisation,      -1, 2,2);
00811       pmanager->AddProcess(new G4MuBremsstrahlung,  -1, 3,3);
00812       pmanager->AddProcess(new G4MuPairProduction,  -1, 4,4);
00813       if(BDSGlobalConstants::Instance()->GetTurnOnCerenkov()){
00814         G4Cerenkov* theCerenkovProcess = new G4Cerenkov;
00815         pmanager->AddProcess(theCerenkovProcess);
00816         pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
00817 #if G4VERSION_NUMBER < 953
00818         pmanager->AddDiscreteProcess(new G4MuNuclearInteraction);     
00819 #else 
00820         /*      pmanager->AddDiscreteProcess(new G4MuonVDNuclearModel); */
00821 #endif
00822       }
00823     }
00824   }
00825 }
00826    
00827 
00828 void BDSPhysicsList::ConstructDecay()
00829 {
00830 #ifdef BDSDEBUG
00831   G4cout << __METHOD_NAME__ << G4endl;
00832 #endif
00833   theParticleIterator->reset();
00834   G4Decay* theDecayProcess = new G4Decay();
00835   while( (*theParticleIterator)() ){
00836     G4ParticleDefinition* particle = theParticleIterator->value();
00837     G4ProcessManager* pmanager = particle->GetProcessManager();
00838     G4String particleName = particle->GetParticleName();
00839     
00840     if (theDecayProcess->IsApplicable(*particle)) { 
00841       pmanager -> AddProcess(theDecayProcess);
00842       pmanager -> SetProcessOrdering(theDecayProcess, idxPostStep);
00843       pmanager -> SetProcessOrdering(theDecayProcess, idxAtRest);
00844     }
00845   }
00846 }
00847 
00848 
00849 void BDSPhysicsList::ConstructOptical()
00850 {
00851 #ifdef BDSDEBUG
00852   G4cout << __METHOD_NAME__ << G4endl;
00853 #endif
00854   bool bCerOn=BDSGlobalConstants::Instance()->GetTurnOnCerenkov();
00855   bool bBirksOn=BDSGlobalConstants::Instance()->GetTurnOnBirksSaturation();
00856 
00857 //  theCerenkovProcess->DumpPhysicsTable();
00858 //  theScintillationProcess->DumpPhysicsTable();
00859 //  theRayleighScatteringProcess->DumpPhysicsTable();
00860 
00861   if(bCerOn){
00862     if(!theCerenkovProcess){
00863       theCerenkovProcess = new G4Cerenkov("Cerenkov");
00864     }
00865   }
00866   
00867   theScintillationProcess        = new G4Scintillation("Scintillation");
00868   if(BDSGlobalConstants::Instance()->GetTurnOnOpticalAbsorption()){
00869     theAbsorptionProcess         = new G4OpAbsorption();
00870   }
00871   if(BDSGlobalConstants::Instance()->GetTurnOnRayleighScattering()){
00872     theRayleighScatteringProcess = new G4OpRayleigh();
00873   }
00874   if(BDSGlobalConstants::Instance()->GetTurnOnMieScattering()){
00875     theMieHGScatteringProcess    = new G4OpMieHG();
00876   }
00877   if(BDSGlobalConstants::Instance()->GetTurnOnOpticalSurface()){
00878     theBoundaryProcess           = new G4OpBoundaryProcess();
00879 #if G4VERSION_NUMBER < 960
00880     G4OpticalSurfaceModel themodel = unified;
00881     theBoundaryProcess->SetModel(themodel);
00882 #endif
00883   }
00884 
00885   SetVerboseLevel(1);
00886   theScintillationProcess->SetScintillationYieldFactor(BDSGlobalConstants::Instance()->GetScintYieldFactor());
00887   theScintillationProcess->SetTrackSecondariesFirst(true);
00888 
00889   // Use Birks Correction in the Scintillation process
00890 
00891   if(bBirksOn){
00892     G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
00893     theScintillationProcess->AddSaturation(emSaturation);
00894   }
00895 
00896   theParticleIterator->reset();
00897   while( (*theParticleIterator)() ){
00898     G4ParticleDefinition* particle = theParticleIterator->value();
00899     G4ProcessManager* pmanager = particle->GetProcessManager();
00900     G4String particleName = particle->GetParticleName();
00901     if(bCerOn){
00902       if (theCerenkovProcess->IsApplicable(*particle)) {
00903         pmanager->AddProcess(theCerenkovProcess);
00904         pmanager->SetProcessOrdering(theCerenkovProcess,idxPostStep);
00905         theCerenkovProcess->SetMaxNumPhotonsPerStep(20);
00906         theCerenkovProcess->SetMaxBetaChangePerStep(10.0);
00907         theCerenkovProcess->SetTrackSecondariesFirst(true);
00908       }
00909     }
00910     if (theScintillationProcess->IsApplicable(*particle)) {
00911       pmanager->AddProcess(theScintillationProcess);
00912       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxAtRest);
00913       pmanager->SetProcessOrderingToLast(theScintillationProcess, idxPostStep);
00914     }
00915     if (particleName == "opticalphoton") {
00916 #ifdef BDSDEBUG
00917       G4cout << "AddDiscreteProcess to OpticalPhoton " << G4endl;
00918 #endif
00919       if(BDSGlobalConstants::Instance()->GetTurnOnOpticalAbsorption()){
00920         pmanager->AddDiscreteProcess(theAbsorptionProcess);
00921       }
00922       if(BDSGlobalConstants::Instance()->GetTurnOnRayleighScattering()){
00923         pmanager->AddDiscreteProcess(theRayleighScatteringProcess);
00924       }
00925       if(BDSGlobalConstants::Instance()->GetTurnOnMieScattering()){
00926         pmanager->AddDiscreteProcess(theMieHGScatteringProcess);
00927       }
00928       if(BDSGlobalConstants::Instance()->GetTurnOnOpticalSurface()){
00929         pmanager->AddDiscreteProcess(theBoundaryProcess);
00930       }
00931     }
00932   }
00933 }
00934 
00935 
00936 void BDSPhysicsList::ConstructMerlin()
00937 {
00938 #ifdef BDSDEBUG
00939   G4cout << __METHOD_NAME__ << G4endl;
00940 #endif
00941   theParticleIterator->reset();
00942   while( (*theParticleIterator)() ){
00943     G4ParticleDefinition* particle = theParticleIterator->value();
00944     G4ProcessManager* pmanager = particle->GetProcessManager();
00945     G4String particleName = particle->GetParticleName();
00946     
00947     if (particleName == "e-") {
00948       //electron
00949       pmanager->AddProcess(new G4eMultipleScattering,-1, 1,1);
00950       pmanager->AddProcess(new G4eIonisation,       -1, 2,2);
00951       pmanager->AddProcess(new G4eBremsstrahlung,   -1, 3,3);      
00952     } 
00953   }
00954 }
00955 
00956 void BDSPhysicsList::ConstructEM_Low_Energy()
00957 {
00958 #ifdef BDSDEBUG
00959   G4cout << __METHOD_NAME__ << G4endl;
00960 #endif
00961   //Applicability range for Livermore models                                                                                                                                
00962   //for higher energies, the Standard models are used                                                                                                                       
00963   G4double highEnergyLimit = 1*CLHEP::GeV;
00964 
00965   theParticleIterator->reset();
00966   while( (*theParticleIterator)() ){
00967     G4ParticleDefinition* particle = theParticleIterator->value();
00968     G4ProcessManager* pmanager = particle->GetProcessManager();
00969     G4String particleName = particle->GetParticleName();
00970      
00971     if (particleName == "gamma") {
00972       G4RayleighScattering* rayl = new G4RayleighScattering();
00973       G4LivermoreRayleighModel*
00974         raylModel = new G4LivermoreRayleighModel();
00975       raylModel->SetHighEnergyLimit(highEnergyLimit);
00976       rayl->AddEmModel(0, raylModel);
00977       pmanager->AddDiscreteProcess(rayl);
00978 
00979       G4PhotoElectricEffect* phot = new G4PhotoElectricEffect();
00980       G4LivermorePhotoElectricModel*
00981         photModel = new G4LivermorePhotoElectricModel();
00982       photModel->SetHighEnergyLimit(highEnergyLimit);
00983       phot->AddEmModel(0, photModel);
00984       pmanager->AddDiscreteProcess(phot);
00985 
00986       G4ComptonScattering* compt = new G4ComptonScattering();
00987       G4LivermoreComptonModel*
00988         comptModel = new G4LivermoreComptonModel();
00989       comptModel->SetHighEnergyLimit(highEnergyLimit);
00990       compt->AddEmModel(0, comptModel);
00991       pmanager->AddDiscreteProcess(compt);
00992       
00993       G4GammaConversion* conv = new G4GammaConversion();
00994       G4LivermoreGammaConversionModel*
00995         convModel = new G4LivermoreGammaConversionModel();
00996       convModel->SetHighEnergyLimit(highEnergyLimit);
00997       conv->AddEmModel(0, convModel);
00998       pmanager->AddDiscreteProcess(conv);
00999       
01000     } else if (particleName == "e-") {
01001         pmanager->AddProcess(new G4eMultipleScattering,-1, 1,1);
01002         G4eIonisation* eIoni = new G4eIonisation();
01003         G4LivermoreIonisationModel*
01004           eIoniModel = new G4LivermoreIonisationModel();
01005         eIoniModel->SetHighEnergyLimit(highEnergyLimit);
01006         eIoni->AddEmModel(0, eIoniModel, new G4UniversalFluctuation() );
01007         pmanager->AddProcess(eIoni,                   -1,-1, 1);
01008         
01009         G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
01010         G4LivermoreBremsstrahlungModel*
01011           eBremModel = new G4LivermoreBremsstrahlungModel();
01012         eBremModel->SetHighEnergyLimit(highEnergyLimit);
01013         eBrem->AddEmModel(0, eBremModel);
01014         pmanager->AddProcess(eBrem,                   -1,-1, 2);
01015             
01016     } else if (particleName == "e+") {
01017         pmanager->AddProcess(new G4eMultipleScattering,-1, 1,1);
01018         pmanager->AddProcess(new G4eIonisation,        -1, 2,2);
01019         pmanager->AddProcess(new G4eBremsstrahlung,    -1, 3,3);
01020         pmanager->AddProcess(new G4eplusAnnihilation,   0,-1,4);
01021     } else if( particleName == "mu+" || 
01022                particleName == "mu-"    ) {
01023         pmanager->AddProcess(new G4MuMultipleScattering,-1, 1,1);
01024         pmanager->AddProcess(new G4MuIonisation,      -1, 2,2);
01025         pmanager->AddProcess(new G4MuBremsstrahlung,  -1, 3,3);
01026         pmanager->AddProcess(new G4MuPairProduction,  -1, 4,4);       
01027 
01028     } else if (particleName == "GenericIon") {
01029         pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
01030         pmanager->AddProcess(new G4ionIonisation,     -1,-1, 1);
01031 
01032     } else if ((!particle->IsShortLived()) &&
01033                (particle->GetPDGCharge() != 0.0) && 
01034                (particle->GetParticleName() != "chargedgeantino")) {
01035 
01036       pmanager->AddProcess(new G4hMultipleScattering,-1, 1,1);
01037       pmanager->AddProcess(new G4hIonisation,       -1,-1, 1);
01038     }
01039   }
01040 }
01041 
01042 void BDSPhysicsList::ConstructLaserWire()
01043 {
01044 #ifdef BDSDEBUG
01045   G4cout << __METHOD_NAME__ << G4endl;
01046 #endif
01047   G4cout << "Constructing laser-wire" << G4endl;
01048 
01049   theParticleIterator->reset();
01050 
01051   BDSLaserCompton* lwProcess = new BDSLaserCompton;
01052 
01053   while( (*theParticleIterator)() ){
01054     G4ParticleDefinition* particle = theParticleIterator->value();
01055     G4ProcessManager* pmanager = particle->GetProcessManager();
01056     G4String particleName = particle->GetParticleName();
01057     
01058     if (particleName == "e-") {
01059       pmanager->AddProcess(lwProcess);
01060       pmanager->SetProcessOrderingToLast(lwProcess,idxPostStep);
01061     }
01062 
01063     if (particleName == "e+") {
01064       pmanager->AddProcess(lwProcess);
01065       pmanager->SetProcessOrderingToLast(lwProcess,idxPostStep);
01066     }
01067 
01068   }
01069 
01070 }
01071 
01072 
01073 // Hadron Processes
01074 
01075 #include "G4HadronElasticProcess.hh"
01076 #include "G4HadronFissionProcess.hh"
01077 #include "G4HadronCaptureProcess.hh"
01078 
01079 #include "G4PionPlusInelasticProcess.hh"
01080 #include "G4PionMinusInelasticProcess.hh"
01081 #include "G4KaonPlusInelasticProcess.hh"
01082 #include "G4KaonZeroSInelasticProcess.hh"
01083 #include "G4KaonZeroLInelasticProcess.hh"
01084 #include "G4KaonMinusInelasticProcess.hh"
01085 #include "G4ProtonInelasticProcess.hh"
01086 #include "G4AntiProtonInelasticProcess.hh"
01087 #include "G4NeutronInelasticProcess.hh"
01088 #include "G4AntiNeutronInelasticProcess.hh"
01089 #include "G4LambdaInelasticProcess.hh"
01090 #include "G4AntiLambdaInelasticProcess.hh"
01091 #include "G4SigmaPlusInelasticProcess.hh"
01092 #include "G4SigmaMinusInelasticProcess.hh"
01093 #include "G4AntiSigmaPlusInelasticProcess.hh"
01094 #include "G4AntiSigmaMinusInelasticProcess.hh"
01095 #include "G4XiZeroInelasticProcess.hh"
01096 #include "G4XiMinusInelasticProcess.hh"
01097 #include "G4AntiXiZeroInelasticProcess.hh"
01098 #include "G4AntiXiMinusInelasticProcess.hh"
01099 #include "G4OmegaMinusInelasticProcess.hh"
01100 #include "G4AntiOmegaMinusInelasticProcess.hh"
01101 
01102 // Low-energy Models
01103 #if G4VERSION_NUMBER < 1000
01104 #include "G4LCapture.hh"
01105 #else
01106 #include "G4HadronElastic.hh"
01107 #include "G4NeutronRadCapture.hh"
01108 #endif
01109 #include "G4LFission.hh"
01110 
01111 #if G4VERSION_NUMBER < 1000
01112 #include "G4LEPionPlusInelastic.hh"
01113 #include "G4LEPionMinusInelastic.hh"
01114 #include "G4LEKaonPlusInelastic.hh"
01115 #include "G4LEKaonZeroSInelastic.hh"
01116 #include "G4LEKaonZeroLInelastic.hh"
01117 #include "G4LEKaonMinusInelastic.hh"
01118 #include "G4LEProtonInelastic.hh"
01119 #include "G4LEAntiProtonInelastic.hh"
01120 #include "G4LENeutronInelastic.hh"
01121 #include "G4LEAntiNeutronInelastic.hh"
01122 #include "G4LELambdaInelastic.hh"
01123 #include "G4LEAntiLambdaInelastic.hh"
01124 #include "G4LESigmaPlusInelastic.hh"
01125 #include "G4LESigmaMinusInelastic.hh"
01126 #include "G4LEAntiSigmaPlusInelastic.hh"
01127 #include "G4LEAntiSigmaMinusInelastic.hh"
01128 #include "G4LEXiZeroInelastic.hh"
01129 #include "G4LEXiMinusInelastic.hh"
01130 #include "G4LEAntiXiZeroInelastic.hh"
01131 #include "G4LEAntiXiMinusInelastic.hh"
01132 #include "G4LEOmegaMinusInelastic.hh"
01133 #include "G4LEAntiOmegaMinusInelastic.hh"
01134 #else
01135 #include "G4CascadeInterface.hh"
01136 #include "G4BinaryLightIonReaction.hh"
01137 #endif
01138 
01139 // -- generator models
01140 #include "G4TheoFSGenerator.hh"
01141 #include "G4ExcitationHandler.hh"
01142 #include "G4GeneratorPrecompoundInterface.hh"
01143 #include "G4StringModel.hh"
01144 #include "G4PreCompoundModel.hh"
01145 #include "G4QGSMFragmentation.hh"
01146 #include "G4ExcitedStringDecay.hh"
01147 
01148 void BDSPhysicsList::ConstructHadronic()
01149 {
01150 #ifdef BDSDEBUG
01151   G4cout << __METHOD_NAME__ << G4endl;
01152 #endif
01153   
01154 #if G4VERSION_NUMBER < 1000
01155   G4NeutronBuilder* theNeutrons=new G4NeutronBuilder;
01156   theNeutrons->RegisterMe(new G4LHEPNeutronBuilder);
01157 
01158   G4ProtonBuilder * thePro=new G4ProtonBuilder;
01159   thePro->RegisterMe(new G4LHEPProtonBuilder);
01160 
01161   G4PiKBuilder * thePiK=new G4PiKBuilder;
01162   thePiK->RegisterMe(new G4LHEPPiKBuilder);
01163 
01164   theNeutrons->Build();
01165   thePro->Build();
01166   thePiK->Build();
01167 
01168   // Photonuclear processes
01169 
01170   G4PhotoNuclearProcess * thePhotoNuclearProcess = new G4PhotoNuclearProcess;
01171   G4GammaNuclearReaction * theGammaReaction = new G4GammaNuclearReaction;
01172   G4ElectronNuclearProcess * theElectronNuclearProcess = new G4ElectronNuclearProcess;
01173   G4PositronNuclearProcess * thePositronNuclearProcess = new G4PositronNuclearProcess;
01174   G4ElectroNuclearReaction * theElectroReaction = new G4ElectroNuclearReaction;
01175   G4TheoFSGenerator * theModel = new G4TheoFSGenerator;
01176   
01177   G4QGSModel< G4GammaParticipants > * theStringModel = new G4QGSModel< G4GammaParticipants >;
01178   G4ExcitedStringDecay * theStringDecay = new G4ExcitedStringDecay(/*theFragmentation=*/new G4QGSMFragmentation);
01179   theStringModel->SetFragmentationModel(theStringDecay);
01180   
01181   G4GeneratorPrecompoundInterface * theCascade = new G4GeneratorPrecompoundInterface;
01182   
01183   theModel->SetTransport(theCascade);
01184   theModel->SetHighEnergyGenerator(theStringModel);
01185 
01186   G4ProcessManager * aProcMan = 0;
01187   
01188   aProcMan = G4Gamma::Gamma()->GetProcessManager();
01189   theGammaReaction->SetMaxEnergy(3.5*CLHEP::GeV);
01190   thePhotoNuclearProcess->RegisterMe(theGammaReaction);
01191   theModel->SetMinEnergy(3.*CLHEP::GeV);
01192   theModel->SetMaxEnergy(100*CLHEP::TeV);
01193   thePhotoNuclearProcess->RegisterMe(theModel);
01194   aProcMan->AddDiscreteProcess(thePhotoNuclearProcess);
01195 
01196   aProcMan = G4Electron::Electron()->GetProcessManager();
01197   theElectronNuclearProcess->RegisterMe(theElectroReaction);
01198   aProcMan->AddDiscreteProcess(theElectronNuclearProcess);
01199   
01200   aProcMan = G4Positron::Positron()->GetProcessManager();
01201   thePositronNuclearProcess->RegisterMe(theElectroReaction);
01202   aProcMan->AddDiscreteProcess(thePositronNuclearProcess);
01203 #else
01204     /*
01205     From (14-7-14)
01206     http://geant4.cern.ch/support/proc_mod_catalog/physics_lists/useCases.shtml  
01207 
01208     LHC neutron fluxes
01209 
01210     The Simulation of neutron fluxes needs, in addition to the simulation of hadronic showers, a transport for low energy neutron down to thermal energies. Recommended: QGSP_BERT_HP, QGSP_BIC_HP
01211     
01212     Linear collider neutron fluxes
01213 
01214     Recommended: QGSP_BERT_HP, QGSP_BIC_HP
01215 
01216     Shielding applications (all energies)
01217 
01218     Recommended: QGSP_BERT_HP, QGSP_BIC_HP, QGSP_INCLXX, Shielding
01219   */
01220 
01221   /* based on these recommendations QGSP_BERT_HP, QGSP_BIC_HP are chosen 
01222      with QGSP_INCLXX and Shielding perhaps optional for certain studies
01223    */
01224 
01225   theHadPhysList1 = new G4HadronPhysicsQGSP_BERT_HP();
01226   theHadPhysList1->ConstructProcess();
01227 
01228   theHadPhysList2 = new G4HadronPhysicsQGSP_BIC_HP();
01229   theHadPhysList2->ConstructProcess();
01230 
01231 #endif
01232 }
01233 
01234 void BDSPhysicsList::ConstructSR()
01235 {
01236   // BDSIM's version of Synchrotron Radiation
01237   BDSSynchrotronRadiation* srProcess = new BDSSynchrotronRadiation;
01238   BDSContinuousSR *contSR = new BDSContinuousSR(); // contin. energy loss process
01239 
01240   // G4's version of Synchrotron Radiation - not used because does not have
01241   // Multiplicity or MeanFreeFactor capability
01242   // G4SynchrotronRadiation* srProcess = new G4SynchrotronRadiation;
01243 
01244   theParticleIterator->reset();
01245 
01246   while( (*theParticleIterator)() ){
01247     G4ParticleDefinition* particle = theParticleIterator->value();
01248     G4ProcessManager* pmanager = particle->GetProcessManager();
01249     G4String particleName = particle->GetParticleName();
01250     
01251     if (particleName == "e-") {
01252       pmanager->AddProcess(srProcess);
01253       pmanager->SetProcessOrderingToLast(srProcess,idxPostStep);
01254 
01255       G4int idx = pmanager->AddProcess(contSR);
01256       pmanager->SetProcessOrderingToLast(contSR,idxPostStep);
01257       pmanager->SetProcessActivation(idx, false);
01258     }
01259     
01260     if (particleName == "e+") {
01261       pmanager->AddProcess(srProcess);
01262       pmanager->SetProcessOrderingToLast(srProcess,idxPostStep);
01263 
01264       // JS : why is contSR switched off for positrons?
01265       
01266       //G4int idx = pmanager->AddProcess(contSR);
01267       //      pmanager->SetProcessOrderingToLast(contSR,idxPostStep);
01268       //      pmanager->SetProcessActivation(idx, false);
01269     }
01270     
01271   }
01272   return; 
01273 }
01274 
01275 void BDSPhysicsList::AddParameterisation()
01276 {
01277   G4FastSimulationManagerProcess*
01278     theFastSimulationManagerProcess =
01279     new G4FastSimulationManagerProcess();
01280   G4cout << "FastSimulationManagerProcess" <<G4endl;
01281   theParticleIterator->reset();
01282   //G4cout<<"---"<<G4endl;                                                                                                                                              
01283   while( (*theParticleIterator)() ){
01284     //G4cout<<"+++"<<G4endl;                                                                                                                                            
01285 
01286     G4ParticleDefinition* particle = theParticleIterator->value();
01287     // G4cout<<"--- particle "<<particle->GetParticleName()<<G4endl;                                                                                                    
01288     G4ProcessManager* pmanager = particle->GetProcessManager();
01289     // The fast simulation process becomes a discrete process only since 9.0:                                                                                                 
01290     pmanager->AddDiscreteProcess(theFastSimulationManagerProcess);
01291   }
01292 }