diff --git a/geant4/LEMuSR/G4Modified/update.sh b/geant4/LEMuSR/G4Modified/update.sh index 50a80c1..a4d4f0a 100644 --- a/geant4/LEMuSR/G4Modified/update.sh +++ b/geant4/LEMuSR/G4Modified/update.sh @@ -17,6 +17,9 @@ cp G4El_UsualEqRhs.cc $G4INSTALL/source/geometry/magneticfield/src/ cp G4El_UsualEqRhs.hh $G4INSTALL/source/geometry/magneticfield/include/ # Enable the muon decay channel for "Mu" particle +# Set parent polarization +cp G4VDecayChannel.cc $G4INSTALL/source/particles/management/src/ +cp G4VDecayChannel.hh $G4INSTALL/source/particles/management/include/ cp G4MuonDecayChannel.cc $G4INSTALL/source/particles/management/src/ cp G4MuonDecayChannel.hh $G4INSTALL/source/particles/management/include/ diff --git a/geant4/LEMuSR/LEMuSR.cc b/geant4/LEMuSR/LEMuSR.cc index 9a4ef1c..102375c 100644 --- a/geant4/LEMuSR/LEMuSR.cc +++ b/geant4/LEMuSR/LEMuSR.cc @@ -100,7 +100,9 @@ int main(int argc,char** argv)//argc:: defines the user interface runManager ->SetUserAction( new LEMuSRRunAction()); runManager ->SetUserAction( new LEMuSREventAction());// scintillators, sensitive detectors runManager ->SetUserAction( new LEMuSRTrackingAction()); + runManager ->SetUserAction( new LEMuSRStackingAction()); + //! Optionnal stepping action: enable one at once only diff --git a/geant4/LEMuSR/LastModifications.txt b/geant4/LEMuSR/LastModifications.txt new file mode 100644 index 0000000..197aa49 --- /dev/null +++ b/geant4/LEMuSR/LastModifications.txt @@ -0,0 +1,28 @@ +LAST MODIFICATIONS + +Stacking Manager: + The stacking manager is now responsible for the tracks selection. If the tracks in the staks are not desired partiles they are killed directly. This was done before by the stepping manager and costed more CPU time. + The gamma particles are killed. It seems that some processes for gamma particles have a bug in the 8.0.2.p01 version of geant4. The consequence is a global time that goes back. One should contact them to check it. + +DetectorConstruction: + The ring anode was entirely built in the half or quarter detector view modes (cf command Detector/View ). This caused few problems with the geometry controler. + + +AtRestSpinPrecession: + The spin precession at rest is not calculated if the magnetic field is (0,0,0). This "mute" problem caused crashed in the simulation under rare special conditions. The typical error when a (0,0,0) vector is used or when a vector is divided by zero leads to a "ZMxpvInifiteVector thrown" message. + + +PhysicsList: + For simplicity the modular physics list was removed and everything is now in LEMuSRPhyiscsList. + + + +SteppingAction: + Some data about the killed particles (killed for looping in fields) are now printed out. + + +LEMuSRMultipleScattering: + This class is now G4 independant and has its own ParticleChange. + + + diff --git a/geant4/LEMuSR/include/LEMuSRCryoSD.hh b/geant4/LEMuSR/include/LEMuSRCryoSD.hh index 2ac6647..bcd39e7 100644 --- a/geant4/LEMuSR/include/LEMuSRCryoSD.hh +++ b/geant4/LEMuSR/include/LEMuSRCryoSD.hh @@ -146,8 +146,27 @@ typedef struct Int_t muon,positron,gamma, runid; } cryoHit ; + +typedef struct +{ + + Int_t muon,positron,gamma, runid; + +} particleHit ; + + +typedef struct +{ + Int_t grid, guards, mcp, cryo; + Float_t RAL, RAR, L3, TDet1, TDet2, TDet3; + Float_t Bmag; +} detectorData ; + //! Root hit object. cryoHit theHit; + particleHit theParticle; + detectorData detector; + //! Filling the Root hit object. /*! diff --git a/geant4/LEMuSR/include/LEMuSRIonPhysics.hh b/geant4/LEMuSR/include/LEMuSRIonPhysics.hh deleted file mode 100644 index 9d4a3d4..0000000 --- a/geant4/LEMuSR/include/LEMuSRIonPhysics.hh +++ /dev/null @@ -1,80 +0,0 @@ -// ------------------------------------------- -// History -// first version 12 Nov. 2000 by H.Kurashige -// ------------------------------------------------------------ -#ifndef LEMuSRIonPhysics_h -#define LEMuSRIonPhysics_h 1 - -#include "globals.hh" -#include "G4ios.hh" - -#include "G4VPhysicsConstructor.hh" - -#include "G4HadronElasticProcess.hh" -#include "G4LElastic.hh" - -#include "G4DeuteronInelasticProcess.hh" -#include "G4LEDeuteronInelastic.hh" - -#include "G4TritonInelasticProcess.hh" -#include "G4LETritonInelastic.hh" - -#include "G4AlphaInelasticProcess.hh" -#include "G4LEAlphaInelastic.hh" - -#include "G4hIonisation.hh" -#include "G4ionIonisation.hh" -#include "G4MultipleScattering.hh" - -class LEMuSRIonPhysics : public G4VPhysicsConstructor -{ - public: - LEMuSRIonPhysics(const G4String& name="ion"); - virtual ~LEMuSRIonPhysics(); - - public: - // This method will be invoked in the Construct() method. - // each particle type will be instantiated - virtual void ConstructParticle(); - - // This method will be invoked in the Construct() method. - // each physics process will be instantiated and - // registered to the process manager of each particle type - virtual void ConstructProcess(); - - protected: - // Elastic Process - G4HadronElasticProcess theElasticProcess; - G4LElastic* theElasticModel; - - // Generic Ion physics - G4MultipleScattering fIonMultipleScattering; - G4ionIonisation fIonIonisation; - - // Deuteron physics - G4MultipleScattering fDeuteronMultipleScattering; - G4hIonisation fDeuteronIonisation; - G4DeuteronInelasticProcess fDeuteronProcess; - G4LEDeuteronInelastic* fDeuteronModel; - - // Triton physics - G4MultipleScattering fTritonMultipleScattering; - G4hIonisation fTritonIonisation; - G4TritonInelasticProcess fTritonProcess; - G4LETritonInelastic* fTritonModel; - - // Alpha physics - G4MultipleScattering fAlphaMultipleScattering; - G4hIonisation fAlphaIonisation; - G4AlphaInelasticProcess fAlphaProcess; - G4LEAlphaInelastic* fAlphaModel; - - // He3 physics - G4MultipleScattering fHe3MultipleScattering; - G4hIonisation fHe3Ionisation; - -}; - - -#endif - diff --git a/geant4/LEMuSR/include/LEMuSRMSC.hh b/geant4/LEMuSR/include/LEMuSRMSC.hh index 89002a1..40efe57 100644 --- a/geant4/LEMuSR/include/LEMuSRMSC.hh +++ b/geant4/LEMuSR/include/LEMuSRMSC.hh @@ -64,7 +64,7 @@ #include "G4GPILSelection.hh" #include "G4PhysicsLogVector.hh" #include "G4VPhysicalVolume.hh" -#include "G4ParticleChangeForMSC.hh" +#include "LEMuSRParticleChangeForMSC.hh" #include "G4UnitsTable.hh" #include "G4MaterialCutsCouple.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -212,7 +212,7 @@ class LEMuSRMSC : public G4VContinuousDiscreteProcess G4double NuclCorrPar; G4double FactPar; - G4ParticleChangeForMSC fParticleChange; + LEMuSRParticleChangeForMSC fParticleChange; G4double alfa1,alfa2,alfa3,b,xsi,c0,facxsi ; // angle distr. parameters // facxsi : some tuning diff --git a/geant4/LEMuSR/include/LEMuSRMagneticField.hh b/geant4/LEMuSR/include/LEMuSRMagneticField.hh index d11c138..ef90701 100644 --- a/geant4/LEMuSR/include/LEMuSRMagneticField.hh +++ b/geant4/LEMuSR/include/LEMuSRMagneticField.hh @@ -20,11 +20,11 @@ // Magnetic Field //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$// -#include "G4UniformMagField.hh" +#include "G4MagneticField.hh" #include"G4ThreeVector.hh" #include"G4ios.hh" -class LEMuSRMagneticField : public G4UniformMagField +class LEMuSRMagneticField : public G4MagneticField { public: diff --git a/geant4/LEMuSR/include/LEMuSRStackingAction.hh b/geant4/LEMuSR/include/LEMuSRStackingAction.hh index 0b583be..1643a61 100644 --- a/geant4/LEMuSR/include/LEMuSRStackingAction.hh +++ b/geant4/LEMuSR/include/LEMuSRStackingAction.hh @@ -64,7 +64,12 @@ private: public: - /* inline void SetNRequestMuon(G4int val) { reqMuon = val; } + G4bool kill_e, kill_gamma, kill_nu_e, kill_nu_mu; + + void KillUnwanted(); + + + /* inline void SetNRequestMuon(G4int val) { reqMuon = val; } inline G4int GetNRequestMuon() const { return reqMuon; } inline void SetNRequestIsoMuon(G4int val) { reqIsoMuon = val; } inline G4int GetNRequestIsoMuon() const { return reqIsoMuon; } diff --git a/geant4/LEMuSR/include/LEMuSRStackingActionMessenger.hh b/geant4/LEMuSR/include/LEMuSRStackingActionMessenger.hh index 0d4cf61..7e841fa 100644 --- a/geant4/LEMuSR/include/LEMuSRStackingActionMessenger.hh +++ b/geant4/LEMuSR/include/LEMuSRStackingActionMessenger.hh @@ -25,8 +25,30 @@ #define LEMuSRStackingActionMessenger_h 1 class LEMuSRStackingAction; -class G4UIcmdWithAnInteger; + +#include "G4UIdirectory.hh" +#include "G4UIcmdWith3VectorAndUnit.hh" +#include "G4UIcmdWith3Vector.hh" +#include "G4UIcmdWithADoubleAndUnit.hh" +#include "G4UIcmdWithAString.hh" +#include "G4UIcmdWithADouble.hh" + +#include "G4UIcmdWithAnInteger.hh" +#include "G4UIcmdWithoutParameter.hh" +#include "G4UIcommand.hh" +#include "G4UImanager.hh" +#include "G4UIterminal.hh" +#include "G4UItcsh.hh" + +class G4UIcommand; +class G4UIdirectory; +class G4UIcmdWithADouble; class G4UIcmdWithADoubleAndUnit; +class G4UIcmdWith3VectorAndUnit; +class G4UIcmdWith3Vector; +class G4UIcmdWithAnInteger; +class G4UIcmdWithAString; +class G4UIcmdWithoutParameter; #include "G4UImessenger.hh" #include "globals.hh" @@ -49,12 +71,13 @@ class LEMuSRStackingActionMessenger: public G4UImessenger LEMuSRStackingAction * myAction; private: //commands - /* G4UIcmdWithAnInteger * muonCmd; - G4UIcmdWithAnInteger * isomuonCmd; - G4UIcmdWithAnInteger * isoCmd; - G4UIcmdWithADoubleAndUnit * roiCmd; + /* + G4UIcmdWithAnInteger * muonCmd; + G4UIcmdWithAnInteger * isomuonCmd; + G4UIcmdWithAnInteger * isoCmd; + G4UIcmdWithADoubleAndUnit * roiCmd; */ - + }; #endif diff --git a/geant4/LEMuSR/include/LEMuSRTrack.hh b/geant4/LEMuSR/include/LEMuSRTrack.hh deleted file mode 100644 index b2de162..0000000 --- a/geant4/LEMuSR/include/LEMuSRTrack.hh +++ /dev/null @@ -1,11 +0,0 @@ -#include "G4Track.hh" -#include "G4DynamicParticle.hh" - -class LEMuSRTrack : public G4Track -{ - -public: - void SetDynamicParticle(G4DynamicParticle ); - - -}; diff --git a/geant4/LEMuSR/src/LEMuSRAtRestSpinRotation.cc b/geant4/LEMuSR/src/LEMuSRAtRestSpinRotation.cc index 932575c..e344427 100644 --- a/geant4/LEMuSR/src/LEMuSRAtRestSpinRotation.cc +++ b/geant4/LEMuSR/src/LEMuSRAtRestSpinRotation.cc @@ -82,7 +82,7 @@ G4VParticleChange* LEMuSRAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack, { theParticleChange.Initialize(theTrack); - theParticleChange.ProposeTrackStatus(fAlive); + theParticleChange.ProposeTrackStatus(fStopButAlive); theParticleChange.ProposePosition(theTrack.GetPosition()); theParticleChange.ProposeMomentumDirection(theTrack.GetMomentumDirection()); theParticleChange.ProposeEnergy(theTrack.GetKineticEnergy()); @@ -128,24 +128,27 @@ G4VParticleChange* LEMuSRAtRestSpinRotation::AtRestDoIt(const G4Track& theTrack, #ifdef G4SRVERBOSE G4cout <<"AT REST::: MAGNETIC FIELD B="<< B[0] <<" " << B[1] <<" " << B[2] <GetCurrentRun()->GetRunID(); + + + theParticle.muon = mu; + theParticle.positron = e; + theParticle.gamma = g; + theParticle.runid = G4RunManager::GetRunManager()->GetCurrentRun()->GetRunID(); + } // ROOT METHODS @@ -211,6 +218,8 @@ void LEMuSRCryoSD::BookRoot() tree->Branch("cryoHit",&theHit.kenergy,"kenergy/F:tenergy/F:edeposit/F:localtime/F:globaltime:propertime/F:positionx/F:positiony:positionz:momdirx:momdiry:momdirz/F:foil/F:muon/I:muonium/I:gamma/I:runID/I"); + tree->Branch("particle",&theParticle.muon,"muon/I:muonium/I:gamma/I:runID/I"); + } void LEMuSRCryoSD::FillRoot() diff --git a/geant4/LEMuSR/src/LEMuSRDetectorConstruction.cc b/geant4/LEMuSR/src/LEMuSRDetectorConstruction.cc index 83b50af..3614fae 100644 --- a/geant4/LEMuSR/src/LEMuSRDetectorConstruction.cc +++ b/geant4/LEMuSR/src/LEMuSRDetectorConstruction.cc @@ -42,6 +42,7 @@ #include "G4PropagatorInField.hh" #include "G4ChordFinder.hh" #include "G4El_UsualEqRhs.hh" +#include "G4Mag_UsualEqRhs.hh" #include "G4ClassicalRK4.hh" #include "G4SimpleHeum.hh" #include "G4SimpleRunge.hh" @@ -169,13 +170,13 @@ G4VPhysicalVolume* LEMuSRDetectorConstruction::lemuDetector() // !mind the V in //! Definition of the top mothe volume: the laboratory. Refered to as Wolrd. // +++++++++++++++++++++++++++++++DEFINE THE MOTHER VOLUME: THE LABORATOY++++++++++++++++++++++ // solid - G4double LABO_x = 1*m; - G4double LABO_y = 1*m; - G4double LABO_z = 2*m; + G4double LABO_x = 2*m; + G4double LABO_y = 2*m; + G4double LABO_z = 4*m; LABO_box = new G4Box("World_box",LABO_x,LABO_y,LABO_z); // logical volume - LABO_material = G4Material::GetMaterial("Air"); + LABO_material = G4Material::GetMaterial("vacuum"); lv_LABO = new G4LogicalVolume(LABO_box,LABO_material,"lv_World",0,0,0); // physical volume @@ -354,7 +355,7 @@ void LEMuSRDetectorConstruction::lemuMCP2() // gate valve chamber - GATV_tube= new G4Tubs( "sGATV",7.65*cm, 10.325*cm, 9.25*cm, dSPhi, dEPhi); + GATV_tube= new G4Tubs( "sGATV",7.655*cm, 10.32*cm, 9.25*cm, dSPhi, dEPhi); GATS_tube= new G4Tubs( "sGATS", 10.325*cm, 12.65*cm, 9.25*cm, dSPhi, dEPhi); G4Transform3D tr_gatv; @@ -391,13 +392,12 @@ void LEMuSRDetectorConstruction::lemuMCP2() + if(AsymCheck!=1){ pv_MCPS = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, 0*cm),lv_MCPS ,"pv_MCPS",lv_LABO, false, 0 ); pv_F160 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, -15.10*cm),lv_F160 ,"pv_F160",lv_LABO, false, 0 ); - pv_F100 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, +16.2*cm),lv_F100 ,"pv_F100",lv_LABO, false, 0 ); + pv_F100 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, +16.2*cm),lv_F100 ,"pv_F100",lv_MCPV, false, 0 ); - G4VPhysicalVolume* pv_F101; - pv_F101 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, -124.*cm),lv_F100 ,"pv_F101",lv_LABO, false, 0 ); GATVz= -25.45; pv_F200 = new G4PVPlacement( 0, G4ThreeVector( 0.*cm, 0.*cm, 8.05*cm),lv_F200 ,"pv_F200",lv_GATV, false, 0 ); @@ -416,7 +416,7 @@ void LEMuSRDetectorConstruction::lemuMCP2() // Apply attributes and user limit - lv_MCPV->SetVisAttributes(G4VisAttributes::Invisible); + lv_MCPV->SetVisAttributes(VTBB_style);//G4VisAttributes::Invisible); lv_MCPS->SetVisAttributes(Blue_style); lv_F160->SetVisAttributes(fBlue_style); lv_F100->SetVisAttributes(fBlue_style); @@ -428,7 +428,6 @@ void LEMuSRDetectorConstruction::lemuMCP2() // user limits G4UserLimits* RAV_lim = new G4UserLimits(); RAV_lim -> SetMaxAllowedStep(FieldStepLim ); - //RAV_lim->SetUserMaxTrackLength(1.2*m); RAV_lim-> SetUserMinEkine(1.*eV); lv_MCPV->SetUserLimits(RAV_lim); // @@ -453,13 +452,13 @@ void LEMuSRDetectorConstruction::lemuANODE() RA_Ebox = new G4Box("RA_Ebox", 1.*cm,12.5*cm,2.25*cm ); RA_Mbox = new G4Box("RA_Mbox", 1.*cm,12.5*cm,2.25*cm ); - RA_G_tube= new G4Tubs( "sRA_G", 5.8*cm,6.25*cm,5.55*cm, dSPhi, dEPhi); + RA_G_tube = new G4Tubs( "sRA_G", 5.8*cm,6.25*cm,5.55*cm, dSPhi, dEPhi); // logical volumes - lv_RA_E = new G4LogicalVolume( RA_E_cone , SSteel ,"lv_RA_E",0,0,0); - lv_RA_M = new G4LogicalVolume( RA_M_cone , SSteel ,"lv_RA_M",0,0,0); - lv_RA_G = new G4LogicalVolume( RA_G_tube , SSteel ,"lv_RA_G",0,0,0); + lv_RA_E = new G4LogicalVolume( RA_E_cone ,G4Material::GetMaterial("stainless_steel") ,"lv_RA_E",0,0,0); + lv_RA_M = new G4LogicalVolume( RA_M_cone ,G4Material::GetMaterial("stainless_steel") ,"lv_RA_M",0,0,0); + lv_RA_G = new G4LogicalVolume( RA_G_tube , G4Material::GetMaterial("copper") ,"lv_RA_G",0,0,0); // physical volumes @@ -469,8 +468,12 @@ void LEMuSRDetectorConstruction::lemuANODE() // Mother Volume is lv_MCPV!!!! - pv_RA_E = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Ez*cm-mcpv_z),lv_RA_E ,"pv_RA_E",lv_MCPV, false, 0 ); - pv_RA_M = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Mz*cm-mcpv_z),lv_RA_M ,"pv_RA_M",lv_MCPV, false, 0 ); + if(halfview==0) + { + pv_RA_E = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Ez*cm-mcpv_z),lv_RA_E ,"pv_RA_E",lv_MCPV, false, 0 ); + pv_RA_M = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Mz*cm-mcpv_z),lv_RA_M ,"pv_RA_M",lv_MCPV, false, 0 ); + } + pv_RA_G = new G4PVPlacement( 0, G4ThreeVector( 0*cm, 0*cm, RA_Gz*cm-mcpv_z),lv_RA_G ,"pv_RA_G",lv_MCPV, false, 0 ); rotation1=G4RotateZ3D(180*deg) ; @@ -493,7 +496,6 @@ void LEMuSRDetectorConstruction::lemuANODE() lv_RA_E->SetVisAttributes(lBlue_style); lv_RA_M->SetVisAttributes(dBlue_style); lv_RA_G->SetVisAttributes(lBlue_style); - } @@ -534,7 +536,7 @@ void LEMuSRDetectorConstruction:: lemuCRYO() Guard_Rings = new G4Tubs ("sGuard",3.0*cm,3.5*cm,0.25*cm, dSPhi, dEPhi); - // logicals volumes + // logicals volumes // sample MAGNETIC field @@ -1164,7 +1166,7 @@ void LEMuSRDetectorConstruction::buildAnodeField() // 2 case of electro magnetic field - G4UniformMagField* mcField; + G4MagneticField* mcField; if(magfield==1) { @@ -1235,7 +1237,7 @@ void LEMuSRDetectorConstruction::buildAnodeField() if(magfield==1) { - G4UniformMagField* mcField; + G4MagneticField* mcField; #ifdef LEMU_TRANSVERSE_FIELD diff --git a/geant4/LEMuSR/src/LEMuSRDetectorMessenger.cc b/geant4/LEMuSR/src/LEMuSRDetectorMessenger.cc index 8568059..2164df4 100644 --- a/geant4/LEMuSR/src/LEMuSRDetectorMessenger.cc +++ b/geant4/LEMuSR/src/LEMuSRDetectorMessenger.cc @@ -403,8 +403,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu { if(newvalue=="total") { - theDetector->dSPhi=0*deg; - theDetector->dEPhi=360*deg; + theDetector->dSPhi=0.*deg; + theDetector->dEPhi=360.*deg; theDetector->halfview=0; newDetector = theDetector->Construct(); @@ -413,8 +413,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu } else if(newvalue=="half") { - theDetector->dSPhi=+90*deg; - theDetector->dEPhi=180*deg; + theDetector->dSPhi=90.*deg; + theDetector->dEPhi=180.*deg; theDetector->halfview=1; newDetector = theDetector->Construct(); @@ -423,8 +423,8 @@ void LEMuSRDetectorMessenger::SetNewValue(G4UIcommand* command, G4String newvalu } else if(newvalue=="quarter") { - theDetector->dSPhi=0*deg; - theDetector->dEPhi=270*deg; + theDetector->dSPhi=0.*deg; + theDetector->dEPhi=270.*deg; theDetector->halfview=1; newDetector = theDetector->Construct(); diff --git a/geant4/LEMuSR/src/LEMuSRIonPhysics.cc b/geant4/LEMuSR/src/LEMuSRIonPhysics.cc deleted file mode 100644 index e1607d3..0000000 --- a/geant4/LEMuSR/src/LEMuSRIonPhysics.cc +++ /dev/null @@ -1,103 +0,0 @@ -// $Id$ -// GEANT4 tag $Name: $ -// -// - -#include "LEMuSRIonPhysics.hh" - -#include "globals.hh" -#include "G4ios.hh" -#include - - -LEMuSRIonPhysics::LEMuSRIonPhysics(const G4String& name) - : G4VPhysicsConstructor(name) -{ -} - -LEMuSRIonPhysics::~LEMuSRIonPhysics() -{ -} - -#include "G4ParticleDefinition.hh" -#include "G4ParticleTable.hh" - -// Nuclei -#include "G4IonConstructor.hh" - -void LEMuSRIonPhysics::ConstructParticle() -{ - // Construct light ions - G4IonConstructor pConstructor; - pConstructor.ConstructParticle(); -} - - -#include "G4ProcessManager.hh" - - -void LEMuSRIonPhysics::ConstructProcess() -{ - G4ProcessManager * pManager = 0; - - // Elastic Process - theElasticModel = new G4LElastic(); - theElasticProcess.RegisterMe(theElasticModel); - - // Generic Ion - pManager = G4GenericIon::GenericIon()->GetProcessManager(); - // add process - pManager->AddDiscreteProcess(&theElasticProcess); - - pManager->AddProcess(&fIonMultipleScattering, -1, 1, 1); - pManager->AddProcess(&fIonIonisation, -1, 2, 2); - - // Deuteron - pManager = G4Deuteron::Deuteron()->GetProcessManager(); - // add process - pManager->AddDiscreteProcess(&theElasticProcess); - - fDeuteronModel = new G4LEDeuteronInelastic(); - fDeuteronProcess.RegisterMe(fDeuteronModel); - pManager->AddDiscreteProcess(&fDeuteronProcess); - - pManager->AddProcess(&fDeuteronMultipleScattering, -1, 1, 1); - pManager->AddProcess(&fDeuteronIonisation, -1, 2, 2); - - // Triton - pManager = G4Triton::Triton()->GetProcessManager(); - // add process - pManager->AddDiscreteProcess(&theElasticProcess); - - fTritonModel = new G4LETritonInelastic(); - fTritonProcess.RegisterMe(fTritonModel); - pManager->AddDiscreteProcess(&fTritonProcess); - - pManager->AddProcess(&fTritonMultipleScattering, -1, 1, 1); - pManager->AddProcess(&fTritonIonisation, -1, 2, 2); - - // Alpha - pManager = G4Alpha::Alpha()->GetProcessManager(); - // add process - pManager->AddDiscreteProcess(&theElasticProcess); - - fAlphaModel = new G4LEAlphaInelastic(); - fAlphaProcess.RegisterMe(fAlphaModel); - pManager->AddDiscreteProcess(&fAlphaProcess); - - pManager->AddProcess(&fAlphaMultipleScattering, -1, 1, 1); - pManager->AddProcess(&fAlphaIonisation, -1, 2, 2); - - // He3 - pManager = G4He3::He3()->GetProcessManager(); - // add process - pManager->AddDiscreteProcess(&theElasticProcess); - - pManager->AddProcess(&fHe3MultipleScattering, -1, 1, 1); - pManager->AddProcess(&fHe3Ionisation, -1, 2, 2); - - -} - - - diff --git a/geant4/LEMuSR/src/LEMuSRMUONIUM.cc b/geant4/LEMuSR/src/LEMuSRMUONIUM.cc index c360b4d..485e387 100644 --- a/geant4/LEMuSR/src/LEMuSRMUONIUM.cc +++ b/geant4/LEMuSR/src/LEMuSRMUONIUM.cc @@ -67,7 +67,7 @@ G4VParticleChange* LEMuSRMUONIUM::PostStepDoIt( { fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ; } - + // G4cout<<"MU Formation"<-20*cm) - { //G4cout<<"true!"; - factor=exp((-Z*Z)/(10*cm*10*cm)); - field[0]= BField.x()*factor ;//TAO - field[1]= BField.y()*factor ; - field[2]= BField.z()*factor ; - - //G4cout<<"true!"<-20*cm) + { //G4cout<<"true!";*/ + factor=exp((-Z*Z)/(10*cm*10*cm)); + field[0]= BField.x()*factor*0 ;//TAO + field[1]= BField.y()*factor*0 ; + field[2]= BField.z()*factor*0 ; + + // G4cout<<"true: "<SetDefaultValue(Hep3Vector(0., 0., -114)); momCmd = new G4UIcmdWith3Vector("/lemuGun/MomentumDirection",this); - momCmd->SetGuidance("Set momentum direction >> norm must be one."); + momCmd->SetGuidance("Set momentum direction >> sum must equal one."); momCmd->SetParameterName("momX ","momY ","momZ ",true,true); momCmd->SetDefaultValue(Hep3Vector(0., 0., 1.)); diff --git a/geant4/LEMuSR/src/LEMuSRPhysicsList.cc b/geant4/LEMuSR/src/LEMuSRPhysicsList.cc index 75763b7..bb6a61e 100644 --- a/geant4/LEMuSR/src/LEMuSRPhysicsList.cc +++ b/geant4/LEMuSR/src/LEMuSRPhysicsList.cc @@ -51,8 +51,8 @@ LEMuSRPhysicsList::LEMuSRPhysicsList(): G4VUserPhysicsList() { - defaultCutValue = 1.0*mm; - SetVerboseLevel(1); + defaultCutValue = 5.0*mm; + SetVerboseLevel(2); } LEMuSRPhysicsList::~LEMuSRPhysicsList() @@ -78,11 +78,13 @@ void LEMuSRPhysicsList::ConstructParticle() // for all particles which you want to use. // This ensures that objects of these particle types will be // created in the program. - + SetCuts(); ConstructBosons(); ConstructLeptons(); ConstructBaryons(); // ConstructMesons(); + + } @@ -169,12 +171,14 @@ void LEMuSRPhysicsList::ConstructBaryons() #include "G4PhotoElectricEffect.hh" #include "G4MultipleScattering.hh" +#include "G4MultipleScattering52.hh" #include "G4eIonisation.hh" #include "G4eBremsstrahlung.hh" #include "G4eplusAnnihilation.hh" #include "G4MuIonisation.hh" +#include "G4MuIonisation52.hh" #include "G4MuBremsstrahlung.hh" #include "G4MuPairProduction.hh" @@ -210,7 +214,9 @@ void LEMuSRPhysicsList::ConstructProcess() */ ConstructEM(); - /*! Construction of the Decay process when applicable. */ + /*! Construction of the Decay process when applicable. + * + * Construction of User Cuts and Step Limitation. */ ConstructGeneral(); } @@ -291,7 +297,7 @@ void LEMuSRPhysicsList::ConstructEM() } - + /*! \anchor muplusphysics*/ //! Muon Plus Physics is implemented here. else if( particleName == "mu+" ) { @@ -299,7 +305,7 @@ void LEMuSRPhysicsList::ConstructEM() // PROCESSES FOR MUON PLUS //! Multiple Scattering including Meyer'Algorithm: LEMuSRMSC - G4VProcess* aMultipleScattering = new LEMuSRMSC(); + G4VProcess* aMultipleScattering = new LEMuSRMSC();//new G4MultipleScattering52(); //! Muonium formation: LEMuSRMUONIUM G4VProcess* aMuoniumFormation = new LEMuSRMUONIUM(); @@ -313,7 +319,7 @@ void LEMuSRPhysicsList::ConstructEM() G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung(); G4VProcess* aPairProduction = new G4MuPairProduction(); G4VProcess* anIonisation = new G4MuIonisation(); - + // ADD WANTED PROCESSES TO PROCESS MANAGER @@ -332,20 +338,20 @@ void LEMuSRPhysicsList::ConstructEM() // SET PROCESSES ORDERING (ALONG STEP. POST STEP. AT REST). // set ordering for AlongStepDoIt pmanager->SetProcessOrdering(aMultipleScattering, idxAlongStep,1); - pmanager->SetProcessOrdering(anIonisation, idxAlongStep,2); - pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,3); - pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,4); + pmanager->SetProcessOrdering(anIonisation, idxAlongStep,3); + pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,4); + pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,5); // set ordering for PostStepDoIt pmanager->SetProcessOrdering(aMultipleScattering, idxPostStep,1); pmanager->SetProcessOrdering(aMuoniumFormation, idxPostStep,2); - pmanager->SetProcessOrdering(anIonisation, idxPostStep,3); + pmanager->SetProcessOrdering(anIonisation, idxPostStep,3); pmanager->SetProcessOrdering(aBremsstrahlung, idxPostStep,4); pmanager->SetProcessOrdering(aPairProduction, idxPostStep,5); // set ordering for AtRestDoIt - pmanager->SetProcessOrdering(aRestSpinRotation, idxAtRest,1); + pmanager->SetProcessOrderingToFirst(aRestSpinRotation, idxAtRest); } @@ -360,42 +366,19 @@ void LEMuSRPhysicsList::ConstructEM() //! Spin Precession Process At Rest: LEMuSRAtRestSpinRotation G4VProcess* aRestSpinRotation = new LEMuSRAtRestSpinRotation(); G4VProcess* aMuScatt = new LEMuSRMUONIUMScatt(); - - /* - //! Unchanged \gf processes - G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung(); - G4VProcess* aPairProduction = new G4MuPairProduction(); - G4VProcess* anIonisation = new G4MuIonisation(); - - */ // ADD WANTED PROCESSES TO PROCESS MANAGER // add lemu processes pmanager->AddProcess(aMuScatt); pmanager->AddProcess(aRestSpinRotation); - - // add \gf processes - // pmanager->AddProcess(aBremsstrahlung); - // pmanager->AddProcess(anIonisation); - // pmanager->AddProcess(aPairProduction); - - - // SET PROCESSES ORDERING (ALONG STEP. POST STEP. AT REST). - // set ordering for AlongStepDoIt - // pmanager->SetProcessOrdering(anIonisation, idxAlongStep,1); - // pmanager->SetProcessOrdering(aBremsstrahlung, idxAlongStep,2); - // pmanager->SetProcessOrdering(aPairProduction, idxAlongStep,3); // set ordering for PostStepDoIt pmanager->SetProcessOrdering(aMuScatt, idxPostStep,1); - // pmanager->SetProcessOrdering(anIonisation, idxPostStep,2); - // pmanager->SetProcessOrdering(aBremsstrahlung, idxPostStep,3); - // pmanager->SetProcessOrdering(aPairProduction, idxPostStep,4); // set ordering for AtRestDoIt - pmanager->SetProcessOrdering(aRestSpinRotation, idxAtRest,1); + pmanager->SetProcessOrderingToFirst(aRestSpinRotation, idxAtRest); } @@ -405,7 +388,7 @@ void LEMuSRPhysicsList::ConstructEM() //! Muon Minus Physics is implemented here. else if( particleName == "mu-" ) { //muon minus - G4VProcess* aMultipleScattering = new G4MultipleScattering(); + G4VProcess* aMultipleScattering = new G4MultipleScattering52(); G4VProcess* aBremsstrahlung = new G4MuBremsstrahlung(); G4VProcess* aPairProduction = new G4MuPairProduction(); G4VProcess* anIonisation = new G4MuIonisation(); @@ -468,7 +451,8 @@ void LEMuSRPhysicsList::ConstructEM() #include "G4Decay.hh" #include "G4UserSpecialCuts.hh" #include "G4DecayTable.hh" - +// Step Limiter Process +#include "G4StepLimiter.hh" /*! \anchor constructdecayprocess*/ //$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ @@ -481,9 +465,12 @@ void LEMuSRPhysicsList::ConstructEM() void LEMuSRPhysicsList::ConstructGeneral() { /*! Declare Decay Processes (from \lemu and \gf).*/ - G4Decay* theDecayProcess = new G4Decay(); + // G4Decay* theDecayProcess = new G4Decay(); LEMuSRDecay* theLEMuSRDecayProcess = new LEMuSRDecay(); - + + /*! Declare User Cut.*/ + G4StepLimiter* aStepLimiter = new G4StepLimiter(); + G4UserSpecialCuts* aUserCut = new G4UserSpecialCuts(); /*! Initialization of the particle iterator.*/ theParticleIterator->reset(); @@ -512,10 +499,14 @@ void LEMuSRPhysicsList::ConstructGeneral() if (theLEMuSRDecayProcess->IsApplicable(*particle)) { pmanager ->AddProcess(theLEMuSRDecayProcess); - // set ordering for PostStepDoIt and AtRestDoIt - pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxPostStep); - pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxAtRest); + pmanager ->AddProcess(aStepLimiter); + + // set ordering + pmanager ->SetProcessOrderingToLast(aStepLimiter, idxPostStep); + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep); + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest); } + pmanager->DumpInfo(); } else if (particle->GetParticleName()=="Mu") { @@ -534,11 +525,15 @@ void LEMuSRPhysicsList::ConstructGeneral() if (theLEMuSRDecayProcess->IsApplicable(*particle)) { pmanager ->AddProcess(theLEMuSRDecayProcess); - // set ordering for PostStepDoIt and AtRestDoIt - pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxPostStep); - pmanager ->SetProcessOrdering(theLEMuSRDecayProcess, idxAtRest); + pmanager ->AddProcess(aStepLimiter); + + // set ordering + // pmanager ->SetProcessOrdering(aStepLimiter, idxPostStep, 1); // not needed for the muonium + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep); + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest); } + pmanager->DumpInfo(); } @@ -553,6 +548,8 @@ void LEMuSRPhysicsList::ConstructGeneral() else { G4ProcessManager* pmanager = particle->GetProcessManager(); + pmanager ->AddProcess(aUserCut); + pmanager ->SetProcessOrdering(aUserCut, idxPostStep); if (theLEMuSRDecayProcess->IsApplicable(*particle)) { pmanager ->AddProcess(theLEMuSRDecayProcess); @@ -580,6 +577,9 @@ void LEMuSRPhysicsList::SetCuts() // These values are used as the default production thresholds // for the world volume. SetCutsWithDefault(); + SetCutValue(1.*mm, "e+"); + SetCutValue(1.*mm, "e-"); + // SetCutValue(1*mm, "proton"); } diff --git a/geant4/LEMuSR/src/LEMuSRPrimaryGeneratorAction.cc b/geant4/LEMuSR/src/LEMuSRPrimaryGeneratorAction.cc index 10f2195..15ff4cd 100644 --- a/geant4/LEMuSR/src/LEMuSRPrimaryGeneratorAction.cc +++ b/geant4/LEMuSR/src/LEMuSRPrimaryGeneratorAction.cc @@ -78,7 +78,7 @@ LEMuSRPrimaryGeneratorAction::LEMuSRPrimaryGeneratorAction() lemuParticleGun->SetNumberOfParticles(1); X=0; Y=0; - Z=0.; + Z=-114.; momX=0.; momY=0.; momZ=1; scan=0; diff --git a/geant4/LEMuSR/src/LEMuSRStackingAction.cc b/geant4/LEMuSR/src/LEMuSRStackingAction.cc index f56e314..b29a073 100644 --- a/geant4/LEMuSR/src/LEMuSRStackingAction.cc +++ b/geant4/LEMuSR/src/LEMuSRStackingAction.cc @@ -38,6 +38,8 @@ LEMuSRStackingAction::LEMuSRStackingAction() :ScintHits(0), stage(0) { theMessenger = new LEMuSRStackingActionMessenger(this); + + kill_e = true; kill_gamma = true; kill_nu_e = kill_nu_mu = true; } LEMuSRStackingAction::~LEMuSRStackingAction() @@ -48,10 +50,50 @@ LEMuSRStackingAction::~LEMuSRStackingAction() G4ClassificationOfNewTrack -LEMuSRStackingAction::ClassifyNewTrack(const G4Track *) +LEMuSRStackingAction::ClassifyNewTrack(const G4Track *track) { G4ClassificationOfNewTrack classification; classification = fUrgent; + + G4String p_name; + p_name = track->GetDefinition()->GetParticleName(); // particle name + + + + if (p_name=="e-") + { + if (kill_e) classification=fKill; + } + + else if (p_name=="gamma") + { + if (kill_gamma) classification=fKill; + } + + else if (p_name=="nu_e") + { + if (kill_nu_e) classification=fKill; + } + + else if (p_name=="anti_nu_e") + { + if (kill_nu_e) classification=fKill; + } + + else if (p_name=="nu_mu") + { + if (kill_nu_mu) classification=fKill; + } + + else if (p_name=="anti_nu_mu") + { + if (kill_nu_mu) classification=fKill; + } + + + + + return classification; } diff --git a/geant4/LEMuSR/src/LEMuSRSteppingAction.cc b/geant4/LEMuSR/src/LEMuSRSteppingAction.cc index 19aa557..9594098 100644 --- a/geant4/LEMuSR/src/LEMuSRSteppingAction.cc +++ b/geant4/LEMuSR/src/LEMuSRSteppingAction.cc @@ -66,6 +66,7 @@ void LEMuSRSteppingAction::UserSteppingAction(const G4Step* aStep) { LoopKiller(aStep); + // G4cout<<"LEMuSR STEPPING ACTION INVOKED!"< Draw(polyline); } + + // ParticleInfo(aStep); + } @@ -105,23 +109,14 @@ void LEMuSRSteppingAction::UserSteppingAction(const G4Step* aStep) void LEMuSRSteppingAction::LoopKiller(const G4Step *aStep) { // loop killa - if(aStep->GetTrack()->GetCurrentStepNumber()>2500) + if(aStep->GetTrack()->GetCurrentStepNumber()>1500) { aStep->GetTrack()->SetTrackStatus(fStopAndKill); - G4cout<<"killed "< 1500 "<GetTrack()->GetDefinition()->GetParticleName()=="e-" - ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="gamma" - ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="nu_e" - ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="anti_nu_e" - ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="nu_mu" - ||aStep->GetTrack()->GetDefinition()->GetParticleName()=="anti_nu_mu") - { - aStep->GetTrack()->SetTrackStatus(fStopAndKill); - } - */ - } @@ -246,8 +241,8 @@ void LEMuSRSteppingAction:: ParticleInfo(const G4Step *aStep) G4double spin= aStep->GetTrack()->GetDynamicParticle()->GetDefinition()->GetPDGSpin(); // spin in units of 1 //b - G4ThreeVector position = aStep->GetPostStepPoint()->GetPosition(); // position - G4ThreeVector momentum = aStep->GetPostStepPoint()->GetMomentum(); // momentum + G4ThreeVector position = aStep->GetTrack()->GetPosition(); // position + G4ThreeVector momentum = aStep->GetTrack()->GetMomentumDirection(); // momentum //c G4double tof = aStep->GetTrack()->GetLocalTime(); // time since track creation G4double globaltime = aStep->GetTrack()->GetGlobalTime();// time since the event in which the track belongs is created. @@ -263,8 +258,8 @@ void LEMuSRSteppingAction:: ParticleInfo(const G4Step *aStep) // G4cout << "Parent particles are : " << aStep->GetTrack()->GetCreatorProcess()->GetProcessName() <<" ;\n "; } - - G4cout << "particle name : " << p_name <<" ;\n " + G4cout.precision(6); + G4cout << "particle name : " << p_name <<" ;\n " << "volume name : " << v_name <<" ;\n " << "next volume name : " << nv_name <<" ;\n " << "spin : " << spin <<" ;\n " diff --git a/geant4/LEMuSR/src/LEMuSRTrack.cc b/geant4/LEMuSR/src/LEMuSRTrack.cc deleted file mode 100644 index f5c9647..0000000 --- a/geant4/LEMuSR/src/LEMuSRTrack.cc +++ /dev/null @@ -1,9 +0,0 @@ -#include "LEMuSRTrack.hh" - - -void LEMuSRTrack::SetDynamicParticle(G4DynamicParticle particle) -{ - G4DynamicParticle *theParticle; - // theParticle = this->GetDynamicParticle(); - *theParticle = particle; -} diff --git a/geant4/LEMuSR/src/LEMuSRdummydets.cc b/geant4/LEMuSR/src/LEMuSRdummydets.cc index b3127b7..c06f66b 100644 --- a/geant4/LEMuSR/src/LEMuSRdummydets.cc +++ b/geant4/LEMuSR/src/LEMuSRdummydets.cc @@ -96,10 +96,12 @@ void LEMuSRDetectorConstruction::lemuAsym() G4Sphere *Asym_sphr = new G4Sphere("sphAsymR",25*cm,25.000001*cm,0.*deg,360*deg,0*deg,90*deg); // logical volumes - lv_Asym = new G4LogicalVolume(Asym_tube,G4Material::GetMaterial("vacuum"),"lv_Asym",0,0,0); + G4String mat = "vacuum"; - lv_AsymL = new G4LogicalVolume(Asym_sphl,G4Material::GetMaterial("vacuum"),"lv_AsymL",0,0,0); - lv_AsymR = new G4LogicalVolume(Asym_sphr,G4Material::GetMaterial("vacuum"),"lv_AsymR",0,0,0); + lv_Asym = new G4LogicalVolume(Asym_tube,G4Material::GetMaterial(mat),"lv_Asym",0,0,0); + + lv_AsymL = new G4LogicalVolume(Asym_sphl,G4Material::GetMaterial(mat),"lv_AsymL",0,0,0); + lv_AsymR = new G4LogicalVolume(Asym_sphr,G4Material::GetMaterial(mat),"lv_AsymR",0,0,0); // G4Translate3D trl =G4ThreeVector(+50.*cm,0.*cm,0.*cm); //uncomment Asym_rotation = G4RotateY3D(90*deg);