Added cos(Theta) angular distribution

This commit is contained in:
nemu 2014-10-21 15:41:31 +02:00
parent 667a345144
commit 3440daafcc
5 changed files with 35 additions and 22 deletions

View File

@ -63,7 +63,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
void SetMomentumSmearing(G4double val) {pSigma=val;}
void SetMomentumBoundary(G4ThreeVector v){pMinAllowed=v[0]; pMaxAllowed=v[1];}
void SetTilt(G4ThreeVector v) {xangle0=v[0]; yangle0=v[1];}
void SetSigmaTilt(G4ThreeVector v) {xangleSigma=v[0]; yangleSigma=v[1];}
void SetSigmaTilt(G4ThreeVector v) {xangleSigma=v[0]; yangleSigma=v[1];zangleSigma=v[2];}
void SetPitch(G4double val) {pitch=val;}
void SetBeamDirection(G4ThreeVector vIniDir);
void SetInitialMuonPolariz(G4ThreeVector vIniPol);
@ -111,7 +111,7 @@ class musrPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
G4double xMaxSource0, yMaxSource0, zMaxSource0; //P.B. 15 Dec 2009
G4double xMaxSource, yMaxSource, zMaxSource; //P.B. 15 Dec 2009
G4double p0, pSigma, pMinAllowed, pMaxAllowed;
G4double xangle0, yangle0, xangleSigma, yangleSigma, pitch;
G4double xangle0, yangle0, xangleSigma, yangleSigma, zangleSigma,pitch;
G4bool UnpolarisedMuonBeam, TransversalyUnpolarisedMuonBeam;
G4double xPolarisIni, yPolarisIni, zPolarisIni;
G4double xDirection, yDirection, zDirection;

View File

@ -153,9 +153,9 @@
###################################################################################
######################### V I S U A L I S A T I O N ##############################
###################################################################################
#/vis/disable
/vis/disable
#/control/execute visFromToni.mac
/control/execute visDawn101.mac
#/control/execute visDawn101.mac
#/control/execute visVRML.mac
###################################################################################
######################### P A R T I C L E G U N #################################
@ -179,4 +179,4 @@
#/gps/ang/maxphi 2 deg
######################## B E A M O N #######################################
#/run/beamOn 1000
/run/beamOn 5
/run/beamOn 5

View File

@ -187,8 +187,8 @@
###################################################################################
######################### V I S U A L I S A T I O N ##############################
###################################################################################
/vis/disable
#/control/execute visDawn201.mac
#/vis/disable
/control/execute vis.mac
#/control/execute visVRML.mac
###################################################################################
######################### P A R T I C L E G U N #################################
@ -213,5 +213,5 @@
###################################################################################
######################## B E A M O N #########################################
###################################################################################
/run/beamOn 8000000
#/run/beamOn 10
#/run/beamOn 8000000
/run/beamOn 10

View File

@ -368,18 +368,18 @@ void musrPhysicsList::ConstructEM()
//cks 14.12.2011 else if (stringProcessName=="G4PenelopeIonisation") pManager->AddProcess(new G4PenelopeIonisation,nr1,nr2,nr3);
//cks 14.12.2011 else if (stringProcessName=="G4PenelopeBremsstrahlung") pManager->AddProcess(new G4PenelopeBremsstrahlung,nr1,nr2,nr3);
//cks 14.12.2011 else if (stringProcessName=="G4PenelopeAnnihilation") pManager->AddProcess(new G4PenelopeAnnihilation,nr1,nr2,nr3);
//cks 14.12.2011 else if (stringProcessName=="G4LowEnergyIonisation") pManager->AddProcess(new G4LowEnergyIonisation,nr1,nr2,nr3);
//cks 14.12.2011 else if (stringProcessName=="G4LowEnergyBremsstrahlung") pManager->AddProcess(new G4LowEnergyBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuMultipleScattering") pManager->AddProcess(new G4MuMultipleScattering,nr1,nr2,nr3);
// else if (stringProcessName=="G4LowEnergyIonisation") pManager->AddProcess(new G4LowEnergyIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4LowEnergyBremsstrahlung") pManager->AddProcess(new G4LowEnergyBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuMultipleScattering") pManager->AddProcess(new G4MuMultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4MuIonisation") pManager->AddProcess(new G4MuIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4MuBremsstrahlung") pManager->AddProcess(new G4MuBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4MuPairProduction") pManager->AddProcess(new G4MuPairProduction,nr1,nr2,nr3);
// cks 2009-06-08 G4StepLimiter and/or G4UserSpecialCuts are needed to activate the "G4UserLimits"
else if (stringProcessName=="G4StepLimiter") pManager->AddProcess(new G4StepLimiter,nr1,nr2,nr3);
else if (stringProcessName=="G4UserSpecialCuts") pManager->AddProcess(new G4UserSpecialCuts,nr1,nr2,nr3);
// else if (stringProcessName=="G4DecayWithSpin") pManager->AddProcess(new G4DecayWithSpin,nr1,nr2,nr3);
// else if (stringProcessName=="G4hIonisation") pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4hLowEnergyIonisation") pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4DecayWithSpin") pManager->AddProcess(new G4DecayWithSpin,nr1,nr2,nr3);
// else if (stringProcessName=="G4hIonisation") pManager->AddProcess(new G4hIonisation,nr1,nr2,nr3);
// else if (stringProcessName=="G4hLowEnergyIonisation") pManager->AddProcess(new G4hLowEnergyIonisation,nr1,nr2,nr3);
else if (stringProcessName=="musrMuFormation") pManager->AddProcess(new musrMuFormation,nr1,nr2,nr3);
else if (stringProcessName=="G4Scintillation") {
G4Scintillation* myScintillation = new G4Scintillation();

View File

@ -260,8 +260,26 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
if (pitch!=0) {yangle += - pitch * (y-y0); }
} // end of the part specific for the muons generated by random rather then from TURTLE
//------- Add new intial angle
G4double sinXangle, sinYangle, phi;
G4double px, py, pz;
if (zangleSigma<0) {
sinXangle=sqrt(G4UniformRand());
phi=2*CLHEP::pi*G4UniformRand();
// Calculate the final momentum
px = p*sinXangle*cos(phi);
py = p*sinXangle*sin(phi);
pz = std::sqrt(p*p - px*px - py*py);
}
if (zangleSigma> 0) {sinXangle=sin(xangle); sinYangle=sin(yangle);
px = p*sinXangle;
py = p*sinYangle;
pz = std::sqrt(p*p - px*px - py*py);
// printf("px, py, pz = %f, %f, %f\n", px, py, pz);
}
// Calculate particle (muon) starting time
G4double ParticleTime; //P.B. 13 May 2009
if (tSigma>0) {ParticleTime = G4RandGauss::shoot(t0,tSigma);} // Gaussian distribution P.B. 13 May 2009
@ -269,12 +287,7 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
else {ParticleTime = t0;} // Point-like P.B. 13 May 2009
// Calculate the final momentum
G4double px, py, pz;
px = p*sin(xangle);
py = p*sin(yangle);
pz = std::sqrt(p*p - px*px - py*py);
// Now rotate so that beam propagates along /gun/direction