LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Kamil Sedlak 2009-05-18

Browse Source 
This is the first version of the muSR simulation code (musrSim)
based on the merged codes of Kamil Sedlak and Toni Shiroka.
It should be a running version of the simulation code, however 
it has not been very well tested, therefore it will probably
need some further development.
This commit is contained in:
sedlak
2009-05-18 09:59:52 +00:00
commit fcd5eea567
66 changed files with 13320 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

200
src/F04ElementField.cc Normal file
View File

@@ -0,0 +1,200 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
#include "G4GeometryManager.hh"
#include "F04ElementField.hh"
#include "F04GlobalField.hh"
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
G4Navigator* F04ElementField::aNavigator;
F04ElementField::F04ElementField(G4ThreeVector c, G4LogicalVolume* lv)
{
elementFieldName="NAME_NOT_DEFINED";
center = c;
minX = minY = minZ = -DBL_MAX;
maxX = maxY = maxZ = DBL_MAX;
// G4cout<<"Kamil: F04GlobalField: addElementField() will be called: "<<G4endl;
F04GlobalField::getObject()->addElementField(this);
color = "1,1,1";
userLimits = new G4UserLimits();
lvolume = lv;
lvolume->SetVisAttributes(getVisAttribute(color));
maxStep = 1*mm;
userLimits->SetMaxAllowedStep(maxStep);
userLimits->SetUserMaxTrackLength(500.*m);
userLimits->SetUserMaxTime(10*ms);
// userLimits->SetUserMinEkine(0.1*MeV);
// userLimits->SetUserMinRange(1*mm);
lvolume->SetUserLimits(userLimits);
}
void F04ElementField::construct()
{
G4Navigator* theNavigator =
G4TransportationManager::GetTransportationManager()->
GetNavigatorForTracking();
if (!aNavigator) {
aNavigator = new G4Navigator();
if ( theNavigator->GetWorldVolume() )
aNavigator->SetWorldVolume(theNavigator->GetWorldVolume());
}
G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
if (!geomManager->IsGeometryClosed()) {
geomManager->OpenGeometry();
geomManager->CloseGeometry(true);
}
G4cout<<"F04ElementField: center="<<center.x()/mm<<", "<<center.y()/mm<<", "<<center.z()/mm<<" mm"<<G4endl;
aNavigator->LocateGlobalPointAndSetup(center,0,false);
G4TouchableHistoryHandle fTouchable = aNavigator->
CreateTouchableHistoryHandle();
G4int depth = fTouchable->GetHistoryDepth();
for (G4int i = 0; i<depth; ++i) {
if(fTouchable->GetVolume()->GetLogicalVolume() == lvolume)break;
fTouchable->MoveUpHistory();
}
// Check that the point of origin of the field matches with the logical volume, to which it is assigned:
G4String volumeName=lvolume->GetName();
if (fTouchable->GetVolume()->GetLogicalVolume() != lvolume) {
char eMessage[200];
sprintf(eMessage,"F04ElementField.cc::construct(): Centre (point of origin) of the field outside the assigned logical volume \"%s\".",
volumeName.c_str());
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
// G4cout<<"+!+!+! global2local VOLUME NAME: "<<fTouchable->GetVolume()->GetName()<<G4endl;
// Set global2local transform. The centre of the transformation is set to the centre of the volume, not
// to the point "center", as one could naively expect. This is corrected a few lines later.
global2local = fTouchable->GetHistory()->GetTopTransform();
// Print out the point of the origin of the field in the local coordinate system of the logical volume:
G4ThreeVector local_center = global2local.TransformPoint(center);
G4cout<<"\t==> "<<elementFieldName<<" in vol.=\""<<volumeName<<"\", center(local coord. system): "
<<local_center.x()/mm<<", "<<local_center.y()/mm<<", "<<local_center.z()/mm<<" mm."<<G4endl;
// Now move the centre of the transformation such that it coincides with the point "center":
global2local*=G4AffineTransform(-local_center);
// After moving the centre of the transformation, the point "local_center" should be set to 0., 0., 0. in
// the following print-out message:
// local_center = global2local.TransformPoint(center);
// G4cout<<"\t==> "<<elementFieldName<<" (volume=\""<<volumeName<<"\", center of the field in local coordinate syst 2: "
// <<local_center.x()/mm<<", "<<local_center.y()/mm<<", "<<local_center.z()/mm<<" mm.)"<<G4endl;
// set global bounding box
G4double local[4], global[4];
G4ThreeVector globalPosition;
local[3] = 0.0;
for (int i=0; i<2; ++i) {
local[0] = (i==0 ? -1.0 : 1.0) * getWidth()/2.;
for (int j=0; j<2; ++j) {
local[1] = (j==0 ? -1.0 : 1.0) * getHeight()/2.;
for (int k=0; k<2; ++k) {
local[2] = (k==0 ? -1.0 : 1.0) * getLength()/2.;
G4ThreeVector localPosition(local[0],local[1],local[2]);
globalPosition =
global2local.Inverse().TransformPoint(localPosition);
global[0] = globalPosition.x();
global[1] = globalPosition.y();
global[2] = globalPosition.z();
setGlobalPoint(global);
}
}
}
}
G4VisAttributes* F04ElementField::getVisAttribute(G4String color)
{
G4VisAttributes* p = NULL;
if(color.size() > 0 &&
(isdigit(color.c_str()[0]) || color.c_str()[0] == '.')) {
G4double red=0.0, green=0.0, blue=0.0;
if (sscanf(color.c_str(),"%lf,%lf,%lf",&red,&green,&blue) == 3) {
p = new G4VisAttributes(true,G4Color(red,green,blue));
} else {
G4cout << " Invalid color " << color << G4endl;
}
}
if (!p) p = new G4VisAttributes(G4VisAttributes::Invisible);
p->SetDaughtersInvisible(false);
return p;
}
void F04ElementField::SetEventNrDependentField(G4double initialField, G4double finalField, G4int nrOfSteps) {
G4double eventNrStep = float(musrParameters::nrOfEventsToBeGenerated)/(nrOfSteps);
G4double fieldStep = (finalField-initialField)/(nrOfSteps-1);
// G4cout<<"musrParameters::nrOfEventsToBeGenerated="<<musrParameters::nrOfEventsToBeGenerated<<G4endl;
// G4cout<<"fieldStep="<<fieldStep<<" eventNrStep="<<eventNrStep<<G4endl;
for (G4int i=1; i<=nrOfSteps; i++) {
G4int eventNumber = int(i*eventNrStep);
G4double field = initialField + i*fieldStep;
changeFieldInStepsMap[eventNumber]=field;
}
G4cout << "Setting field in steps for field "<<elementFieldName<<G4endl;
std::map<G4int,G4double>::iterator it;
for ( it=changeFieldInStepsMap.begin() ; it != changeFieldInStepsMap.end(); it++ ) {
G4cout << "Field will be changed at event "<< (*it).first << " to the value of " << (*it).second/tesla<<" T" << G4endl;
// G4double nominalFieldValue=it->second;
// it->SetNominalFieldValue(nominalFieldValue);
}
}
void F04ElementField::SetElementFieldValueIfNeeded(G4int eventNr) {
std::map<G4int,G4double>::iterator itr;
itr = changeFieldInStepsMap.find(eventNr);
if (itr==F04ElementField::changeFieldInStepsMap.end()) {
// eventNr was not found in the map ==> field is not going to change at this eventNr
}
else {
G4double newFieldValue = itr->second;
SetNominalFieldValue(newFieldValue);
// G4cout<<"Nominal Field changed for "<<G4endl;
}
}

147
src/F04FieldMessenger.cc Normal file
View File

@@ -0,0 +1,147 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
#include "F04FieldMessenger.hh"
#include "F04GlobalField.hh"
#include "G4UIdirectory.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithoutParameter.hh"
//////////////////////////////////////////////////////////////////////////////
F04FieldMessenger::F04FieldMessenger(F04GlobalField* pEMfield)
: fGlobalField(pEMfield)
{
detDir = new G4UIdirectory("/field/");
detDir->SetGuidance(" Field tracking control ");
fStepperCMD = new G4UIcmdWithAnInteger("/field/setStepperType",this);
fStepperCMD->SetGuidance("Select stepper type for field");
fStepperCMD->SetParameterName("choice",true);
fStepperCMD->SetDefaultValue(4);
fStepperCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fUpdateCMD = new G4UIcmdWithoutParameter("/field/update",this);
fUpdateCMD->SetGuidance("Update Field");
fUpdateCMD->SetGuidance("This command MUST be applied before \"beamOn\" ");
fUpdateCMD->SetGuidance("if you changed field settings.");
fUpdateCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fMinStepCMD = new G4UIcmdWithADoubleAndUnit("/field/setMinStep",this);
fMinStepCMD->SetGuidance("Define minimal step");
fMinStepCMD->SetParameterName("min step",false,false);
fMinStepCMD->SetDefaultUnit("mm");
fMinStepCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fDeltaChordCMD = new G4UIcmdWithADoubleAndUnit("/field/setDeltaChord",this);
fDeltaChordCMD->SetGuidance("Define delta chord");
fDeltaChordCMD->SetParameterName("delta chord",false,false);
fDeltaChordCMD->SetDefaultUnit("mm");
fDeltaChordCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fDeltaOneStepCMD =
new G4UIcmdWithADoubleAndUnit("/field/setDeltaOneStep",this);
fDeltaOneStepCMD->SetGuidance("Define delta one step");
fDeltaOneStepCMD->SetParameterName("delta one step",false,false);
fDeltaOneStepCMD->SetDefaultUnit("mm");
fDeltaOneStepCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fDeltaIntersectionCMD =
new G4UIcmdWithADoubleAndUnit("/field/setDeltaIntersection",this);
fDeltaIntersectionCMD->SetGuidance("Define delta intersection");
fDeltaIntersectionCMD->SetParameterName("delta intersection",false,false);
fDeltaIntersectionCMD->SetDefaultUnit("mm");
fDeltaIntersectionCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fEpsMinCMD = new G4UIcmdWithADoubleAndUnit("/field/setEpsMin",this);
fEpsMinCMD->SetGuidance("Define eps min");
fEpsMinCMD->SetParameterName("eps min",false,false);
fEpsMinCMD->SetDefaultUnit("mm");
fEpsMinCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
fEpsMaxCMD = new G4UIcmdWithADoubleAndUnit("/field/setEpsMax",this);
fEpsMaxCMD->SetGuidance("Define eps max");
fEpsMaxCMD->SetParameterName("eps max",false,false);
fEpsMaxCMD->SetDefaultUnit("mm");
fEpsMaxCMD->AvailableForStates(G4State_PreInit,G4State_Idle);
}
F04FieldMessenger::~F04FieldMessenger()
{
delete detDir;
delete fStepperCMD;
delete fMinStepCMD;
delete fDeltaChordCMD;
delete fDeltaOneStepCMD;
delete fDeltaIntersectionCMD;
delete fEpsMinCMD;
delete fEpsMaxCMD;
delete fUpdateCMD;
}
void F04FieldMessenger::SetNewValue( G4UIcommand* command, G4String newValue)
{
if( command == fStepperCMD )
{
fGlobalField->SetStepperType(fStepperCMD->GetNewIntValue(newValue));
}
if( command == fUpdateCMD )
{
fGlobalField->updateField();
}
if( command == fMinStepCMD )
{
fGlobalField->SetMinStep(fMinStepCMD->GetNewDoubleValue(newValue));
}
if( command == fDeltaChordCMD )
{
fGlobalField->SetDeltaChord(fDeltaChordCMD->GetNewDoubleValue(newValue));
}
if( command == fDeltaOneStepCMD )
{
fGlobalField->
SetDeltaOneStep(fDeltaOneStepCMD->GetNewDoubleValue(newValue));
}
if( command == fDeltaIntersectionCMD )
{
fGlobalField->
SetDeltaIntersection(fDeltaIntersectionCMD->GetNewDoubleValue(newValue));
}
if( command == fEpsMinCMD )
{
fGlobalField->SetEpsMin(fEpsMinCMD->GetNewDoubleValue(newValue));
}
if( command == fEpsMaxCMD )
{
fGlobalField->SetEpsMax(fEpsMaxCMD->GetNewDoubleValue(newValue));
}
}

330
src/F04GlobalField.cc Normal file
View File

@@ -0,0 +1,330 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
//
#include <time.h>
#include "Randomize.hh"
#include "G4TransportationManager.hh"
#include "G4ExplicitEuler.hh"
#include "G4ImplicitEuler.hh"
#include "G4SimpleRunge.hh"
#include "G4SimpleHeum.hh"
#include "G4ClassicalRK4.hh"
#include "G4CashKarpRKF45.hh"
#include "F04GlobalField.hh"
// #include "F04GlobalFieldOldPoints.hh" //cks F04GlobalFieldOldPoints were introduced in order to keep the last few
// values of the field, so that they do not have to be recalculated, if the
// field at recent point is requested. It turned out not to be useful
// (i.e. it did not speed-up the simulation)
#include "musrRootOutput.hh"
F04GlobalField* F04GlobalField::object = 0;
F04GlobalField::F04GlobalField() : G4ElectroMagneticField(),
minStep(0.01*mm), deltaChord(3.0*mm),
deltaOneStep(0.01*mm), deltaIntersection(0.1*mm),
epsMin(2.5e-7*mm), epsMax(0.05*mm),
fEquation(0), fFieldManager(0),
fFieldPropagator(0), fStepper(0), fChordFinder(0)
// F04GlobalField::F04GlobalField() : G4MagneticField(),
// minStep(0.01*mm), deltaChord(3.0*mm),
// deltaOneStep(0.01*mm), deltaIntersection(0.1*mm),
// epsMin(2.5e-7*mm), epsMax(0.05*mm),
// fEquation(0), fFieldManager(0),
// fFieldPropagator(0), fStepper(0), fChordFinder(0)
{
fFieldMessenger = new F04FieldMessenger(this);
fields = new FieldList();
fStepperType = 4 ; // ClassicalRK4 is default stepper
// set object
object = this;
updateField();
//cks myOldFieldPoints = F04GlobalFieldOldPoints();
}
F04GlobalField::~F04GlobalField()
{
// cks 2009_05_14
// For some reason, the "clear()" on the next line was causing the crash in some cases
// (probably when more "Elementary fields" were created). The problem appeared
// when running the 1050.mac file for the lem4. Fixed by commenting out the clear():
// clear();
delete fFieldMessenger;
if (fEquation) delete fEquation;
// if (fFieldManager) delete fFieldManager;
// if (fFieldPropagator) delete fFieldPropagator;
if (fStepper) delete fStepper;
if (fChordFinder) delete fChordFinder;
}
void F04GlobalField::updateField()
{
first = true;
nfp = 0;
fp = 0;
clear();
// Construct equ. of motion of particles through B fields
// fEquation = new G4Mag_EqRhs(this);
// Construct equ. of motion of particles through e.m. fields
// fEquation = new G4EqMagElectricField(this);
// Construct equ. of motion of particles including spin through B fields
// fEquation = new G4Mag_SpinEqRhs(this);
// Construct equ. of motion of particles including spin through e.m. fields
fEquation = new G4EqEMFieldWithSpin(this);
// Get transportation, field, and propagator managers
G4TransportationManager* fTransportManager =
G4TransportationManager::GetTransportationManager();
fFieldManager = GetGlobalFieldManager();
fFieldPropagator = fTransportManager->GetPropagatorInField();
// Need to SetFieldChangesEnergy to account for a time varying electric
// field (r.f. fields)
fFieldManager->SetFieldChangesEnergy(true);
// Set the field
fFieldManager->SetDetectorField(this);
// Choose a stepper for integration of the equation of motion
SetStepper();
// Create a cord finder providing the (global field, min step length,
// a pointer to the stepper)
fChordFinder = new G4ChordFinder((G4MagneticField*)this,minStep,fStepper);
// Set accuracy parameters
fChordFinder->SetDeltaChord( deltaChord );
fFieldManager->SetAccuraciesWithDeltaOneStep(deltaOneStep);
fFieldManager->SetDeltaIntersection(deltaIntersection);
fFieldPropagator->SetMinimumEpsilonStep(epsMin);
fFieldPropagator->SetMaximumEpsilonStep(epsMax);
// G4cout << "Accuracy Parameters:" <<
// " MinStep=" << minStep <<
// " DeltaChord=" << deltaChord <<
// " DeltaOneStep=" << deltaOneStep << G4endl;
// G4cout << " " <<
// " DeltaIntersection=" << deltaIntersection <<
// " EpsMin=" << epsMin <<
// " EpsMax=" << epsMax << G4endl;
fFieldManager->SetChordFinder(fChordFinder);
}
F04GlobalField* F04GlobalField::getObject()
{
if (!object) new F04GlobalField();
return object;
}
void F04GlobalField::SetStepper()
{
if(fStepper) delete fStepper;
switch ( fStepperType )
{
case 0:
// fStepper = new G4ExplicitEuler( fEquation, 8 ); // no spin tracking
fStepper = new G4ExplicitEuler( fEquation, 12 ); // with spin tracking
G4cout << "G4ExplicitEuler is called" << G4endl;
break;
case 1:
// fStepper = new G4ImplicitEuler( fEquation, 8 ); // no spin tracking
fStepper = new G4ImplicitEuler( fEquation, 12 ); // with spin tracking
G4cout << "G4ImplicitEuler is called" << G4endl;
break;
case 2:
// fStepper = new G4SimpleRunge( fEquation, 8 ); // no spin tracking
fStepper = new G4SimpleRunge( fEquation, 12 ); // with spin tracking
G4cout << "G4SimpleRunge is called" << G4endl;
break;
case 3:
// fStepper = new G4SimpleHeum( fEquation, 8 ); // no spin tracking
fStepper = new G4SimpleHeum( fEquation, 12 ); // with spin tracking
G4cout << "G4SimpleHeum is called" << G4endl;
break;
case 4:
// fStepper = new G4ClassicalRK4( fEquation, 8 ); // no spin tracking
fStepper = new G4ClassicalRK4( fEquation, 12 ); // with spin tracking
G4cout << "G4ClassicalRK4 (default) is called" << G4endl;
break;
case 5:
// fStepper = new G4CashKarpRKF45( fEquation, 8 ); // no spin tracking
fStepper = new G4CashKarpRKF45( fEquation, 12 ); // with spin tracking
G4cout << "G4CashKarpRKF45 is called" << G4endl;
break;
default: fStepper = 0;
}
}
G4FieldManager* F04GlobalField::GetGlobalFieldManager()
{
return G4TransportationManager::GetTransportationManager()
->GetFieldManager();
}
void F04GlobalField::GetFieldValue(const G4double* point, G4double* field) const
{
// NOTE: this routine dominates the CPU time for tracking.
// Using the simple array fp[] instead of fields[]
// directly sped it up
// G4cout<<"GlobalField: Field requested at point ("<<point[0]<<","<<point[1]<<","<<point[2]<<")"<<G4endl;
//cks Check whether the requested point is the same as in the previous call
// if (myOldFieldPoints.ThisPointWasCalculatedRecently(point,field)) {return;}
field[0] = field[1] = field[2] = field[3] = field[4] = field[5] = 0.0;
// protect against Geant4 bug that calls us with point[] NaN.
if(point[0] != point[0]) return;
// (can't use nfp or fp, as they may change)
if (first) ((F04GlobalField*)this)->setupArray(); // (cast away const)
for (int i=0; i<nfp; ++i) {
const F04ElementField* p = fp[i];
if (p->isInBoundingBox(point)) {
p->addFieldValue(point,field);
}
}
// cks NOT SURE WHETHER THE FOLLOWING WAS STILL NEEDED, HOWEVER REMOVED NOT
// TO DISTURB THE LOW ENERGY MUONS.
// Set some small field if field is almost zero (to avoid internal problems
// of Geant at almost zero magnetic fields).
// if (sqrt(field[0]*field[0]+field[1]*field[1]+field[2]*field[2])<0.00001*tesla) {
// field[2] = 0.00001*tesla;
// }
// cks myOldFieldPoints.StoreTheFieldPointForFutureReuse(point,field);
}
void F04GlobalField::clear()
{
if (fields) {
if (fields->size()>0) {
FieldList::iterator i;
//cks 2009_05_14 : The following line seems to cause problems (sometimes)
// See the comment in F04GlobalField::~F04GlobalField() for more details.
for (i=fields->begin(); i!=fields->end(); ++i) delete *i;
fields->clear();
}
}
if (fp) delete[] fp;
first = true;
nfp = 0;
fp = NULL;
}
void F04GlobalField::setupArray()
{
first = false;
nfp = fields->size();
fp = new const F04ElementField* [nfp+1]; // add 1 so it's never 0
for (int i=0; i<nfp; ++i) {
fp[i] = (*fields)[i];
//
// Find the event numbers, for which the field changes, for each Element Field.
// Then mark these event numbers in the Global Field, such that it can be efficiently
// find out during the run-time, at which event number the field may change.
std::map<G4int,G4double> localChangeFieldInStepsMap = fp[i] ->GetEventNrDependentField();
std::map<G4int,G4double>::iterator it;
for ( it=localChangeFieldInStepsMap.begin() ; it != localChangeFieldInStepsMap.end(); it++ ) {
G4int eventNr = it->first;
G4cout<<"globalChangeFieldInStepsMap set for event number "<<eventNr<<G4endl;
globalChangeFieldInStepsMap[eventNr]=true;
}
}
musrRootOutput::GetRootInstance()->SetNrFieldNomVal(nfp);
}
void F04GlobalField::CheckWhetherAnyNominalFieldValueNeedsToBeChanged(G4int eventNumber) {
if (globalChangeFieldInStepsMap[eventNumber]) {
// G4cout<<"We should check each Element Field Object whether its field needs to be changed:"<<G4endl;
G4int jjj=0;
FieldList::iterator i;
for (i=fields->begin(); i!=fields->end(); ++i) {
// Set the nominal field value for the given field, if that has been requested for the given field
(*i)->SetElementFieldValueIfNeeded(eventNumber);
// Get the nominal field value for the given field and store it in the Root output
G4double nomFieldValue = (*i)->GetNominalFieldValue();
musrRootOutput::GetRootInstance()->SetFieldNomVal(jjj,nomFieldValue);
jjj++;
}
}
}
// Print field value at all points requested by the user:
void F04GlobalField::PrintFieldAtRequestedPoints() const {
G4ThreeVector p;
G4double point[4];
G4double Bfi[6]={0,0,0,0,0,0};
for (unsigned int i=0; i < pointsAtWhichUserWantsToPrintFieldValue.size(); i++) {
p = pointsAtWhichUserWantsToPrintFieldValue[i];
point[0] = p.x();
point[1] = p.y();
point[2] = p.z();
point[3] = 0.;
object->GetFieldValue(point,Bfi);
// printf (" Magnetic Field at %f, %f, %f mm is B= %10.10f, %10.10f, %10.10f tesla.\n",
// point[0]/mm,point[1]/mm,point[2]/mm,Bfi[0]/tesla,Bfi[1]/tesla,Bfi[2]/tesla);
printf (" EM field value at (%.f, %.f, %.f) mm is: B = (%0.10g, %0.10g, %0.10g) T, E = (%0.10g, %0.10g, %0.10g) kV/mm\n",
point[0]/mm, point[1]/mm, point[2]/mm,
Bfi[0]/tesla, Bfi[1]/tesla, Bfi[2]/tesla,
Bfi[3]/(kilovolt/mm), Bfi[4]/(kilovolt/mm), Bfi[5]/(kilovolt/mm));
}
}

220
src/G4DecayWithSpin.cc Normal file
View File

@@ -0,0 +1,220 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History:
// 17 August 2004 P. Gumplinger, T. MacPhail
// ------------------------------------------------------------
//
#include "G4DecayWithSpin.hh"
#include "G4Step.hh"
#include "G4Track.hh"
#include "G4DecayTable.hh"
#include "G4MuonDecayChannelWithSpin.hh"
#include "G4Vector3D.hh"
#include "G4TransportationManager.hh"
#include "G4PropagatorInField.hh"
#include "G4FieldManager.hh"
#include "G4Field.hh"
#include "G4Transform3D.hh"
//cks
#include "musrErrorMessage.hh"
#include "G4RunManager.hh"
//#include "musrRootOutput.hh"
//csk
G4DecayWithSpin::G4DecayWithSpin(const G4String& processName):G4Decay(processName){}
//cks G4DecayWithSpin::G4DecayWithSpin(const G4String& processName):G4Decay(processName){SetVerboseLevel(3);}
G4DecayWithSpin::~G4DecayWithSpin(){}
G4VParticleChange* G4DecayWithSpin::DecayIt(const G4Track& aTrack, const G4Step& aStep)
{
// debugging: G4RunManager* fRunManager = G4RunManager::GetRunManager();
// debugging: G4int eventNr=fRunManager->GetCurrentEvent()->GetEventID();
// debugging: if ((eventNr>447715)&&(eventNr<45005)) SetVerboseLevel(3);
// debugging: if (eventNr>=45005) SetVerboseLevel(0);
// get particle
const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
// get parent_polarization
G4ThreeVector parent_polarization = aParticle->GetPolarization();
if(parent_polarization == G4ThreeVector(0,0,0))
{
// Generate random polarization direction
G4double cost = 1. - 2.*G4UniformRand();
G4double sint = std::sqrt((1.-cost)*(1.+cost));
G4double phi = twopi*G4UniformRand();
G4double sinp = std::sin(phi);
G4double cosp = std::cos(phi);
G4double px = sint*cosp;
G4double py = sint*sinp;
G4double pz = cost;
parent_polarization.setX(px);
parent_polarization.setY(py);
parent_polarization.setZ(pz);
}else{
G4FieldManager* fieldMgr = aStep.GetTrack()->GetVolume()->
GetLogicalVolume()->GetFieldManager();
if (!fieldMgr) {
G4TransportationManager *transportMgr =
G4TransportationManager::GetTransportationManager();
G4PropagatorInField* fFieldPropagator =
transportMgr->GetPropagatorInField();
if (fFieldPropagator) fieldMgr =
fFieldPropagator->GetCurrentFieldManager();
}
const G4Field* field = NULL;
if(fieldMgr)field = fieldMgr->GetDetectorField();
// if (field && !(fieldMgr->DoesFieldChangeEnergy())) {
if (field) {
G4double point[4];
point[0] = (aStep.GetPostStepPoint()->GetPosition())[0];
point[1] = (aStep.GetPostStepPoint()->GetPosition())[1];
point[2] = (aStep.GetPostStepPoint()->GetPosition())[2];
point[3] = aTrack.GetGlobalTime();
G4double fieldValue[6];
field -> GetFieldValue(point,fieldValue);
G4ThreeVector B(fieldValue[0],fieldValue[1],fieldValue[2]);
if ((B.mag2())>0) { // Call the spin precession only for non-zero mag. field
parent_polarization = Spin_Precession(aStep,B,fRemainderLifeTime);
// G4cout<<"G4DecayWithSpin: parent_polarization="<<parent_polarization<<", B="<<B[0]/tesla<<","<<B[1]/tesla<<","<<B[2]/tesla
// <<", fRemainderLifeTime="<<fRemainderLifeTime<<G4endl;
}
}
}
// decay table
G4DecayTable *decaytable = aParticleDef->GetDecayTable();
if (decaytable) {
G4MuonDecayChannelWithSpin *decaychannel;
decaychannel = (G4MuonDecayChannelWithSpin*)decaytable->SelectADecayChannel();
if (decaychannel) decaychannel->SetPolarization(parent_polarization);
}
G4ParticleChangeForDecay* pParticleChangeForDecay;
pParticleChangeForDecay = (G4ParticleChangeForDecay*)G4Decay::DecayIt(aTrack,aStep);
pParticleChangeForDecay->ProposePolarization(parent_polarization);
// pParticleChangeForDecay->DumpInfo();
return pParticleChangeForDecay;
}
G4ThreeVector G4DecayWithSpin::Spin_Precession( const G4Step& aStep,
G4ThreeVector B, G4double deltatime )
{
G4double Bnorm = std::sqrt(sqr(B[0]) + sqr(B[1]) +sqr(B[2]) );
G4double q = aStep.GetTrack()->GetDefinition()->GetPDGCharge();
G4double a = 1.165922e-3;
G4double s_omega = 8.5062e+7*rad/(s*kilogauss);
G4double omega = -(q*s_omega)*(1.+a) * Bnorm;
G4double rotationangle = deltatime * omega;
G4Transform3D SpinRotation = G4Rotate3D(rotationangle,B.unit());
G4Vector3D Spin = aStep.GetTrack() -> GetPolarization();
//cks 2008.10.09. Check that the polarisatin is normalised do 1
G4double spinMagnitude = Spin.mag();
if ((spinMagnitude>1.1)||(spinMagnitude<0.9)) {
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4int eventNr=fRunManager->GetCurrentEvent()->GetEventID();
musrErrorMessage::GetInstance()->musrError(SERIOUS,
"G4DecayWithSpin.cc::Spin_Precession: spin magnitude is large than 1.1 or smaller than 0.9! Strange!!!",true);
G4cout<<"G4DecayWithSpin.cc::Spin_Precession: "
<<"Event nr.:"<<eventNr
<<", spin = ("<<Spin.x()<<","<<Spin.y()<<","<<Spin.z()<<")."<<G4endl;
G4cout.flush();
}
Spin = Spin.unit();
//csk
G4Vector3D newSpin = SpinRotation * Spin;
//cks 2008.10.09. Check that the polarisatin is normalised do 1
G4double newSpinMagnitude = newSpin.mag();
if ((newSpinMagnitude>1.1)||(newSpinMagnitude<0.9)) {
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4int eventNr=fRunManager->GetCurrentEvent()->GetEventID();
musrErrorMessage::GetInstance()->musrError(SERIOUS,
"G4DecayWithSpin.cc::Spin_Precession: NEW spin magnitude is large than 1.1 or smaller than 0.9! Strange!!!",true);
G4cout<<"G4DecayWithSpin.cc::Spin_Precession: "
<<"Event nr.:"<<eventNr
<<", NEW spin = ("<<newSpin.x()<<","<<newSpin.y()<<","<<newSpin.z()<<")."<<G4endl;
G4cout.flush();
}
newSpin = newSpin.unit();
//csk
#ifdef G4VERBOSE
if (GetVerboseLevel()>2) {
G4double normspin = std::sqrt(Spin*Spin);
G4double normnewspin = std::sqrt(newSpin*newSpin);
//G4double cosalpha = Spin*newSpin/normspin/normnewspin;
//G4double alpha = std::acos(cosalpha);
G4cout << "AT REST::: PARAMETERS " << G4endl;
G4cout << "Initial spin : " << Spin << G4endl;
G4cout << "Delta time : " << deltatime << G4endl;
G4cout << "Rotation angle: " << rotationangle/rad << G4endl;
G4cout << "New spin : " << newSpin << G4endl;
G4cout << "Checked norms : " << normspin <<" " << normnewspin << G4endl;
}
#endif
return newSpin;
}

177
src/G4EqEMFieldWithSpin.cc Normal file
View File

@@ -0,0 +1,177 @@
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: G4EqEMFieldWithSpin.cc,v 1.1 2008/10/21 14:50:32 sedlak Exp $
// GEANT4 tag $Name: $
//
//
// This is the standard right-hand side for equation of motion.
//
// The only case another is required is when using a moving reference
// frame ... or extending the class to include additional Forces,
// eg an electric field
//
// 30.08.2007 Chris Gong, Peter Gumplinger
//
// -------------------------------------------------------------------
#include "G4EqEMFieldWithSpin.hh"
#include "G4ThreeVector.hh"
#include "globals.hh"
//cks
#include "musrErrorMessage.hh"
#include "G4RunManager.hh"
#include "musrRootOutput.hh"
//csk
G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField )
: G4EquationOfMotion( emField ) { anomaly = 1.165923e-3; }
void
G4EqEMFieldWithSpin::SetChargeMomentumMass(G4double particleCharge, // e+ units
G4double MomentumXc,
G4double particleMass)
{
fElectroMagCof = eplus*particleCharge*c_light ;
fMassCof = particleMass*particleMass ;
omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
ParticleCharge = particleCharge;
E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
beta = MomentumXc/E;
gamma = E/particleMass;
}
void
G4EqEMFieldWithSpin::EvaluateRhsGivenB(const G4double y[],
const G4double Field[],
G4double dydx[] ) const
{
// Components of y:
// 0-2 dr/ds,
// 3-5 dp/ds - momentum derivatives
G4double pSquared = y[3]*y[3] + y[4]*y[4] + y[5]*y[5] ;
G4double Energy = std::sqrt( pSquared + fMassCof );
G4double cof2 = Energy/c_light ;
G4double pModuleInverse = 1.0/std::sqrt(pSquared) ;
// G4double inverse_velocity = Energy * c_light * pModuleInverse;
G4double inverse_velocity = Energy * pModuleInverse / c_light;
G4double cof1 = fElectroMagCof*pModuleInverse ;
// G4double vDotE = y[3]*Field[3] + y[4]*Field[4] + y[5]*Field[5] ;
dydx[0] = y[3]*pModuleInverse ;
dydx[1] = y[4]*pModuleInverse ;
dydx[2] = y[5]*pModuleInverse ;
dydx[3] = cof1*(cof2*Field[3] + (y[4]*Field[2] - y[5]*Field[1])) ;
dydx[4] = cof1*(cof2*Field[4] + (y[5]*Field[0] - y[3]*Field[2])) ;
dydx[5] = cof1*(cof2*Field[5] + (y[3]*Field[1] - y[4]*Field[0])) ;
dydx[6] = dydx[8] = 0.;//not used
// Lab Time of flight
dydx[7] = inverse_velocity;
G4ThreeVector BField(Field[0],Field[1],Field[2]);
G4ThreeVector u(y[3], y[4], y[5]);
u *= pModuleInverse;
G4double udb = anomaly*beta*gamma/(1.+gamma) * (BField * u);
G4double ucb = (anomaly+1./gamma)/beta;
G4ThreeVector Spin(y[9],y[10],y[11]);
G4ThreeVector dSpin;
dSpin = ParticleCharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
// cks 2008.10.08
// Check that Spin is not extremely high. The problem is, that in some events
// (for extremely low moementum muons?) the MomentumXc is of the order of 1.e-07,
// the ucb is extremely high (e.g. ucb = 690618107.) and then dSpin becomes
// also extremely high. This problem cummulates and results in an unreasonably
// high spin. For some reason Geant then ends up in an infinite loop.
//
// G4double spinMagnitude = Spin.mag();
// if (spinMagnitude>1.01) dSpin /= spinMagnitude;
// if (spinMagnitude>10.) {
// dSpin=G4ThreeVector(0.,0.,0.);
// G4RunManager* fRunManager = G4RunManager::GetRunManager();
// G4int eventNr=fRunManager->GetCurrentEvent()->GetEventID();
// if (eventNr!=musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction) {
// musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction = eventNr;
// musrErrorMessage::GetInstance()->musrError(SERIOUS,
// "G4EqEMFieldWithSpin.cc:EvaluateRhsGivenB: spin is large than 10.! Strange!!! -> Deleting this event",true);
// G4cout<<" + G4EqEMFieldWithSpin.cc:EvaluateRhsGivenB: "
// <<"Event nr.:"<<eventNr<<", particle mass = "<<sqrt(fMassCof)/MeV<<" MeV"
// <<", spin = ("<<Spin.x()<<","<<Spin.y()<<","<<Spin.z()<<")."<<G4endl;
// G4cout.flush();
// musrRootOutput::GetRootInstance()->SetEventWeight(0);
// fRunManager->AbortEvent();
// }
// }
// if ((spinMagnitude<0.99)&&(spinMagnitude>0.1)&&( fabs(sqrt(fMassCof)/MeV-105.6584)<1) ) dSpin /= spinMagnitude;
// if ((spinMagnitude<0.2)&&( fabs(sqrt(fMassCof)/MeV-105.6584)<1) ) {
// dSpin=G4ThreeVector(0.,0.,0.);
// G4RunManager* fRunManager = G4RunManager::GetRunManager();
// G4int eventNr=fRunManager->GetCurrentEvent()->GetEventID();
// if (eventNr!=musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction) {
// musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction = eventNr;
// musrErrorMessage::GetInstance()->musrError(SERIOUS,
// "G4EqEMFieldWithSpin.cc:EvaluateRhsGivenB: spin for muon is smaller than 0.2! Strange!!! -> Deleting this event",true);
// G4cout<<" - G4EqEMFieldWithSpin.cc:EvaluateRhsGivenB: "
// <<"Event nr.:"<<eventNr<<", particle mass = "<<sqrt(fMassCof)/MeV<<" MeV"
// <<", spin = ("<<Spin.x()<<","<<Spin.y()<<","<<Spin.z()<<")."<<G4endl;
// G4cout.flush();
// musrRootOutput::GetRootInstance()->SetEventWeight(0);
// fRunManager->AbortEvent();
// }
// }
//
// csk
dydx[ 9] = dSpin.x();
dydx[10] = dSpin.y();
dydx[11] = dSpin.z();
return ;
}

202
src/MuDecayChannel.cc Normal file
View File

@@ -0,0 +1,202 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: MuDecayChannel.cc,v 1.17 2006/06/29 19:25:34 gunter Exp $
// GEANT4 tag $Name: geant4-09-00 $
//
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History: first implementation, based on object model of
// 30 May 1997 H.Kurashige
//
// Fix bug in calcuration of electron energy in DecayIt 28 Feb. 01 H.Kurashige
//2005
// M. Melissas ( melissas AT cppm.in2p3.fr)
// J. Brunner ( brunner AT cppm.in2p3.fr)
// Adding V-A fluxes for neutrinos using a new algortithm :
// ------------------------------------------------------------
#include "G4ParticleDefinition.hh"
#include "G4DecayProducts.hh"
#include "G4VDecayChannel.hh"
#include "MuDecayChannel.hh"
#include "Randomize.hh"
#include "G4LorentzVector.hh"
#include "G4LorentzRotation.hh"
#include "G4RotationMatrix.hh"
MuDecayChannel::MuDecayChannel(const G4String& theParentName,
G4double theBR)
:G4VDecayChannel("Muonium Decay",1)
{
// set names for daughter particles
if (theParentName == "Mu") {
SetBR(theBR);
SetParent("Mu");
SetNumberOfDaughters(3);
SetDaughter(0, "e+");
SetDaughter(1, "nu_e");
SetDaughter(2, "anti_nu_mu");
} else {
#ifdef G4VERBOSE
if (GetVerboseLevel()>0) {
G4cout << "MuDecayChannel:: constructor :";
G4cout << " parent particle is not muon but ";
G4cout << theParentName << G4endl;
}
#endif
}
}
MuDecayChannel::~MuDecayChannel()
{
}
G4DecayProducts *MuDecayChannel::DecayIt(G4double)
{
// this version neglects muon polarization,and electron mass
// assumes the pure V-A coupling
// the Neutrinos are correctly V-A.
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) G4cout << "MuDecayChannel::DecayIt ";
#endif
if (parent == 0) FillParent();
if (daughters == 0) FillDaughters();
// parent mass
G4double parentmass = parent->GetPDGMass();
//daughters'mass
G4double daughtermass[3];
G4double sumofdaughtermass = 0.0;
for (G4int index=0; index<3; index++){
daughtermass[index] = daughters[index]->GetPDGMass();
sumofdaughtermass += daughtermass[index];
}
//create parent G4DynamicParticle at rest
G4ThreeVector dummy;
G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
//create G4Decayproducts
G4DecayProducts *products = new G4DecayProducts(*parentparticle);
delete parentparticle;
// calculate daughter momentum
G4double daughtermomentum[3];
// calcurate electron energy
G4double xmax = (1.0+daughtermass[0]*daughtermass[0]/parentmass/parentmass);
G4double x;
G4double Ee,Ene;
G4double gam;
G4double EMax=parentmass/2-daughtermass[0];
//Generating Random Energy
do {
Ee=G4UniformRand();
do{
x=xmax*G4UniformRand();
gam=G4UniformRand();
}while (gam >x*(1.-x));
Ene=x;
} while ( Ene < (1.-Ee));
G4double Enm=(2.-Ee-Ene);
//initialisation of rotation parameters
G4double costheta,sintheta,rphi,rtheta,rpsi;
costheta= 1.-2./Ee-2./Ene+2./Ene/Ee;
sintheta=std::sqrt(1.-costheta*costheta);
rphi=twopi*G4UniformRand()*rad;
rtheta=(std::acos(2.*G4UniformRand()-1.));
rpsi=twopi*G4UniformRand()*rad;
G4RotationMatrix rot;
rot.set(rphi,rtheta,rpsi);
//electron 0
daughtermomentum[0]=std::sqrt(Ee*Ee*EMax*EMax+2.0*Ee*EMax * daughtermass[0]);
G4ThreeVector direction0(0.0,0.0,1.0);
direction0 *= rot;
G4DynamicParticle * daughterparticle = new G4DynamicParticle ( daughters[0], direction0 * daughtermomentum[0]);
products->PushProducts(daughterparticle);
//electronic neutrino 1
daughtermomentum[1]=std::sqrt(Ene*Ene*EMax*EMax+2.0*Ene*EMax * daughtermass[1]);
G4ThreeVector direction1(sintheta,0.0,costheta);
direction1 *= rot;
G4DynamicParticle * daughterparticle1 = new G4DynamicParticle ( daughters[1], direction1 * daughtermomentum[1]);
products->PushProducts(daughterparticle1);
//muonic neutrino 2
daughtermomentum[2]=std::sqrt(Enm*Enm*EMax*EMax +2.0*Enm*EMax*daughtermass[2]);
G4ThreeVector direction2(-Ene/Enm*sintheta,0,-Ee/Enm-Ene/Enm*costheta);
direction2 *= rot;
G4DynamicParticle * daughterparticle2 = new G4DynamicParticle ( daughters[2],
direction2 * daughtermomentum[2]);
products->PushProducts(daughterparticle2);
// output message
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) {
G4cout << "MuDecayChannel::DecayIt ";
G4cout << " create decay products in rest frame " <<G4endl;
products->DumpInfo();
}
#endif
return products;
}

271
src/MuDecayChannelWithSpin.cc Normal file
View File

@@ -0,0 +1,271 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
// ------------------------------------------------------------
// GEANT 4 class header file
//
// History:
// 17 August 2004 P.Gumplinger and T.MacPhail
// samples Michel spectrum including 1st order
// radiative corrections
// Reference: Florian Scheck "Muon Physics", in Physics Reports
// (Review Section of Physics Letters) 44, No. 4 (1978)
// 187-248. North-Holland Publishing Company, Amsterdam
// at page 210 cc.
//
// W.E. Fisher and F. Scheck, Nucl. Phys. B83 (1974) 25.
//
// ------------------------------------------------------------
//
#include "MuDecayChannelWithSpin.hh"
#include "Randomize.hh"
#include "G4DecayProducts.hh"
#include "G4LorentzVector.hh"
MuDecayChannelWithSpin::MuDecayChannelWithSpin(const G4String& theParentName,
G4double theBR)
: MuDecayChannel(theParentName,theBR)
{
}
MuDecayChannelWithSpin::~MuDecayChannelWithSpin()
{
}
G4DecayProducts *MuDecayChannelWithSpin::DecayIt(G4double)
{
// This version assumes V-A coupling with 1st order radiative correctons,
// the standard model Michel parameter values, but
// gives incorrect energy spectrum for neutrinos
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) G4cout << "MuDecayChannelWithSpin::DecayIt ";
#endif
if (parent == 0) FillParent();
if (daughters == 0) FillDaughters();
// parent mass
G4double parentmass = parent->GetPDGMass();
EMMU = parentmass;
//daughters'mass
G4double daughtermass[3];
G4double sumofdaughtermass = 0.0;
for (G4int index=0; index<3; index++){
daughtermass[index] = daughters[index]->GetPDGMass();
sumofdaughtermass += daughtermass[index];
}
EMASS = daughtermass[0];
//create parent G4DynamicParticle at rest
G4ThreeVector dummy;
G4DynamicParticle * parentparticle = new G4DynamicParticle( parent, dummy, 0.0);
//create G4Decayproducts
G4DecayProducts *products = new G4DecayProducts(*parentparticle);
delete parentparticle;
// calcurate electron energy
G4double michel_rho = 0.75; //Standard Model Michel rho
G4double michel_delta = 0.75; //Standard Model Michel delta
G4double michel_xsi = 1.00; //Standard Model Michel xsi
G4double michel_eta = 0.00; //Standard Model eta
G4double rndm, x, ctheta;
G4double FG;
G4double FG_max = 2.00;
G4double W_mue = (EMMU*EMMU+EMASS*EMASS)/(2.*EMMU);
G4double x0 = EMASS/W_mue;
G4double x0_squared = x0*x0;
// ***************************************************
// x0 <= x <= 1. and -1 <= y <= 1
//
// F(x,y) = f(x)*g(x,y); g(x,y) = 1.+g(x)*y
// ***************************************************
// ***** sampling F(x,y) directly (brute force) *****
do{
// Sample the positron energy by sampling from F
rndm = G4UniformRand();
x = x0 + rndm*(1.-x0);
G4double x_squared = x*x;
G4double F_IS, F_AS, G_IS, G_AS;
F_IS = 1./6.*(-2.*x_squared+3.*x-x0_squared);
F_AS = 1./6.*std::sqrt(x_squared-x0_squared)*(2.*x-2.+std::sqrt(1.-x0_squared));
G_IS = 2./9.*(michel_rho-0.75)*(4.*x_squared-3.*x-x0_squared);
G_IS = G_IS + michel_eta*(1.-x)*x0;
G_AS = 3.*(michel_xsi-1.)*(1.-x);
G_AS = G_AS+2.*(michel_xsi*michel_delta-0.75)*(4.*x-4.+std::sqrt(1.-x0_squared));
G_AS = 1./9.*std::sqrt(x_squared-x0_squared)*G_AS;
F_IS = F_IS + G_IS;
F_AS = F_AS + G_AS;
// *** Radiative Corrections ***
G4double R_IS = F_c(x,x0);
G4double F = 6.*F_IS + R_IS/std::sqrt(x_squared-x0_squared);
// *** Radiative Corrections ***
G4double R_AS = F_theta(x,x0);
rndm = G4UniformRand();
ctheta = 2.*rndm-1.;
G4double G = 6.*F_AS - R_AS/std::sqrt(x_squared-x0_squared);
FG = std::sqrt(x_squared-x0_squared)*F*(1.+(G/F)*ctheta);
if(FG>FG_max){
G4cout<<"***Problem in Muon Decay *** : FG > FG_max"<<G4endl;
FG_max = FG;
}
rndm = G4UniformRand();
}while(FG<rndm*FG_max);
G4double energy = x * W_mue;
rndm = G4UniformRand();
G4double phi = twopi * rndm;
if(energy < EMASS) energy = EMASS;
// calculate daughter momentum
G4double daughtermomentum[3];
daughtermomentum[0] = std::sqrt(energy*energy - EMASS*EMASS);
G4double stheta = std::sqrt(1.-ctheta*ctheta);
G4double cphi = std::cos(phi);
G4double sphi = std::sin(phi);
//Coordinates of the decay positron with respect to the muon spin
G4double px = stheta*cphi;
G4double py = stheta*sphi;
G4double pz = ctheta;
G4ThreeVector direction0(px,py,pz);
direction0.rotateUz(parent_polarization);
G4DynamicParticle * daughterparticle0
= new G4DynamicParticle( daughters[0], daughtermomentum[0]*direction0);
products->PushProducts(daughterparticle0);
// daughter 1 ,2 (neutrinos)
// create neutrinos in the C.M frame of two neutrinos
G4double energy2 = parentmass*(1.0 - x/2.0);
G4double vmass = std::sqrt((energy2-daughtermomentum[0])*(energy2+daughtermomentum[0]));
G4double beta = -1.0*daughtermomentum[0]/energy2;
G4double costhetan = 2.*G4UniformRand()-1.0;
G4double sinthetan = std::sqrt((1.0-costhetan)*(1.0+costhetan));
G4double phin = twopi*G4UniformRand()*rad;
G4double sinphin = std::sin(phin);
G4double cosphin = std::cos(phin);
G4ThreeVector direction1(sinthetan*cosphin,sinthetan*sinphin,costhetan);
G4DynamicParticle * daughterparticle1
= new G4DynamicParticle( daughters[1], direction1*(vmass/2.));
G4DynamicParticle * daughterparticle2
= new G4DynamicParticle( daughters[2], direction1*(-1.0*vmass/2.));
// boost to the muon rest frame
G4LorentzVector p4;
p4 = daughterparticle1->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle1->Set4Momentum(p4);
p4 = daughterparticle2->Get4Momentum();
p4.boost( direction0.x()*beta, direction0.y()*beta, direction0.z()*beta);
daughterparticle2->Set4Momentum(p4);
products->PushProducts(daughterparticle1);
products->PushProducts(daughterparticle2);
daughtermomentum[1] = daughterparticle1->GetTotalMomentum();
daughtermomentum[2] = daughterparticle2->GetTotalMomentum();
// output message
#ifdef G4VERBOSE
if (GetVerboseLevel()>1) {
G4cout << "MuDecayChannelWithSpin::DecayIt ";
G4cout << " create decay products in rest frame " <<G4endl;
products->DumpInfo();
}
#endif
return products;
}
G4double MuDecayChannelWithSpin::R_c(G4double x){
G4int n_max = (int)(100.*x);
if(n_max<10)n_max=10;
G4double L2 = 0.0;
for(G4int n=1; n<=n_max; n++){
L2 += std::pow(x,n)/(n*n);
}
G4double omega = std::log(EMMU/EMASS);
G4double r_c;
r_c = 2.*L2-(pi*pi/3.)-2.;
r_c = r_c + omega * (1.5+2.*std::log((1.-x)/x));
r_c = r_c - std::log(x)*(2.*std::log(x)-1.);
r_c = r_c + (3.*std::log(x)-1.-1./x)*std::log(1.-x);
return r_c;
}

1032
src/meyer.cc Normal file
View File

File diff suppressed because it is too large Load Diff

1173
src/musrDetectorConstruction.cc Normal file
View File

File diff suppressed because it is too large Load Diff

168
src/musrDetectorMessenger.cc Normal file
View File

@@ -0,0 +1,168 @@
#include "musrDetectorMessenger.hh"
#include "musrDetectorConstruction.hh"
#include "G4UIdirectory.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWithoutParameter.hh"
#include "G4UIcmdWith3Vector.hh"
#include "G4UIcmdWith3VectorAndUnit.hh"
#include "G4RunManager.hh" //cks included in order to be able to change run ID
#include "Randomize.hh" //cks included in order to initialise the random nr. generator by time
#include <time.h> //cks -----------------------------||-------------------------------
#include "musrEventAction.hh" // cks needed for setting the variable "nHowOftenToPrintEvent"
#include "musrPrimaryGeneratorAction.hh" // cks needed for the initialisation of the random nr. generator by event nr.
#include <vector>
#include "globals.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrDetectorMessenger::musrDetectorMessenger(musrDetectorConstruction* myDet)
:myDetector(myDet)
{
musrDir = new G4UIdirectory("/musr/");
musrDir->SetGuidance("UI commands specific to this example.");
Ignore1Cmd = new G4UIcmdWithAString("/musr/ignore",this);
Ignore1Cmd->SetGuidance("This command is ignored by the messenger, but used for the detector construction.");
Ignore1Cmd->AvailableForStates(G4State_PreInit,G4State_Idle);
Ignore2Cmd = new G4UIcmdWithAString("/musr/command",this);
Ignore2Cmd->SetGuidance("This command is ignored by the messenger, but used for the detector construction.");
Ignore2Cmd->AvailableForStates(G4State_PreInit,G4State_Idle);
runDir = new G4UIdirectory("/musr/run/");
runDir->SetGuidance("musr run control");
RunIDSetCmd = new G4UIcmdWithAnInteger("/musr/run/runID",this);
RunIDSetCmd->SetGuidance("Set the run number");
RunIDSetCmd->SetParameterName("something",false);
RunIDSetCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
RandomOptionCmd = new G4UIcmdWithAnInteger("/musr/run/randomOption",this);
RandomOptionCmd->SetGuidance("Specify the random number generator initialisation");
RandomOptionCmd->SetGuidance(" 0 ... no initialisation (default)");
RandomOptionCmd->SetGuidance(" 1 ... use actual computer time to initialise now");
RandomOptionCmd->SetGuidance(" 2 ... use event number to initialise at the beginning of each event");
RandomOptionCmd->SetGuidance(" 3 ... read in the random no. initial values for each event from a file");
RandomOptionCmd->SetGuidance(" 4 ... read in the random no. initial values for each event from the file kamil.rndm");
RandomOptionCmd->SetParameterName("randomOpt",false);
RandomOptionCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
HowOftenToPrintEventCmd = new G4UIcmdWithAnInteger("/musr/run/howOftenToPrintEvent",this);
HowOftenToPrintEventCmd->SetGuidance("Each n-th event will be notified. Set _n_ by this command.");
HowOftenToPrintEventCmd->SetParameterName("howOftenToPrintEv",false);
HowOftenToPrintEventCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
RndmEventToSaveSeedsCmd = new G4UIcmdWithAnInteger("/musr/run/rndmEventToSaveSeeds",this);
RndmEventToSaveSeedsCmd -> SetGuidance("Save seeds of the random number generators of the given event number.");
RndmEventToSaveSeedsCmd -> SetParameterName("rndmEventToSaveSe",false);
RndmEventToSaveSeedsCmd -> AvailableForStates(G4State_PreInit,G4State_Idle);
detDir = new G4UIdirectory("/musr/det/");
detDir->SetGuidance("detector control.");
UpdateCmd = new G4UIcmdWithoutParameter("/musr/det/update",this);
UpdateCmd->SetGuidance("Update calorimeter geometry.");
UpdateCmd->SetGuidance("This command MUST be applied before \"beamOn\" ");
UpdateCmd->SetGuidance("if you changed geometrical value(s).");
UpdateCmd->AvailableForStates(G4State_Idle);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrDetectorMessenger::~musrDetectorMessenger()
{
delete UpdateCmd;
delete detDir;
delete musrDir;
delete Ignore1Cmd;
delete Ignore2Cmd;
delete RunIDSetCmd;
delete RandomOptionCmd;
delete HowOftenToPrintEventCmd;
delete RndmEventToSaveSeedsCmd;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue) {
if( command == UpdateCmd )
{ myDetector->UpdateGeometry(); }
if( command == RunIDSetCmd )
{ (G4RunManager::GetRunManager())->SetRunIDCounter(RunIDSetCmd->GetNewIntValue(newValue));}
if( command == RandomOptionCmd )
{
G4int RandomOption=RandomOptionCmd->GetNewIntValue(newValue);
if (RandomOption == 1) {
// G4long seed=time(0); //returns time in seconds as an integer
// HepRandom::setTheSeed(seed);//changes the seed of the random engine
G4cout << "******************************************" << G4endl;
G4cout << "*** Random Seed set by the system time ***" << G4endl;
G4cout << "******************************************" << G4endl;
long seeds[2];
time_t systime = time(NULL);
seeds[0] = (long) systime;
seeds[1] = (long) (systime*G4UniformRand());
G4cout << "seed1: " << seeds[0] << "; seed2: " << seeds[1] << G4endl;
CLHEP::HepRandom::setTheSeeds(seeds);
CLHEP::HepRandom::showEngineStatus();
}
else if (RandomOption == 2) {
G4cout << "*******************************************" << G4endl;
G4cout << "*** Random Seed set by the event number ***" << G4endl;
G4cout << "*******************************************" << G4endl;
// musrEventAction::setRandomNrSeedAccordingEventNr=1;
musrPrimaryGeneratorAction::setRandomNrSeedAccordingEventNr=1;
// musrEventAction::setMyEventNr(70);
}
else if (RandomOption == 3) {
G4cout << "*******************************************" << G4endl;
G4cout << "*** Random Seed set from external file ***" << G4endl;
G4cout << "*******************************************" << G4endl;
// musrEventAction::setRandomNrSeedFromFile=1;
musrPrimaryGeneratorAction::setRandomNrSeedFromFile=1;
std::ifstream indata;
int num;
indata.open("randomNum.dat"); // opens the file
if(!indata) { // file couldn't be opened
G4cout << "Error: file could not be opened" << G4endl;
exit(1);
}
std::vector<int> * seedVector = musrPrimaryGeneratorAction::GetPointerToSeedVector();
indata >> num;
while ( !indata.eof() ) { // keep reading until end-of-file
G4cout << "The next number is " << num << G4endl;
seedVector->push_back(num);
indata >> num; // sets EOF flag if no value found
}
indata.close();
G4cout << "End-of-file reached.." << seedVector->size()<<G4endl;
}
else if (RandomOption == 4) {
G4cout << "*********************************************" << G4endl;
G4cout << "*** Random Seed set from kamil.rndm file ***" << G4endl;
G4cout << "*********************************************" << G4endl;
musrPrimaryGeneratorAction::setRandomNrSeedFromFile_RNDM=1;
}
}
if ( command == HowOftenToPrintEventCmd )
{
G4int n = HowOftenToPrintEventCmd->GetNewIntValue(newValue);
musrEventAction::nHowOftenToPrintEvent=n;
}
if ( command == RndmEventToSaveSeedsCmd )
{
G4int n = RndmEventToSaveSeedsCmd->GetNewIntValue(newValue);
musrPrimaryGeneratorAction::nRndmEventToSaveSeeds=n;
}
}

62
src/musrErrorMessage.cc Normal file
View File

@@ -0,0 +1,62 @@
#include "musrErrorMessage.hh"
musrErrorMessage::musrErrorMessage():nErrors(1)
{
pointerToErrors=this;
severityWord[INFO]="INFO";
severityWord[WARNING]="WARNING";
severityWord[SERIOUS]="SERIOUS";
severityWord[FATAL]="FATAL";
}
musrErrorMessage::~musrErrorMessage() {}
musrErrorMessage* musrErrorMessage::pointerToErrors=NULL;
musrErrorMessage* musrErrorMessage::GetInstance() {
return pointerToErrors;
}
void musrErrorMessage::musrError(SEVERITY severity, G4String message, G4bool silent) {
std::map<G4String,ErrorStruct>::iterator it;
it = ErrorMapping.find(message);
if (it == ErrorMapping.end()) { // The error message is called for the first time
ErrorStruct actualErrorMessage;
actualErrorMessage.mesSeverity = severity;
actualErrorMessage.nTimes = 1;
ErrorMapping[message]=actualErrorMessage;
G4cout<<"!!!"<<severityWord[severity]<<"!!! "<<message<<" (First time occurrence)"<<G4endl;
nErrors++;
}
else { // The error message is called for more than the first time
(*it).second.nTimes++;
}
// Print out the error message if required
if ((!silent)||(severity==FATAL)) {
if ((*it).second.nTimes>1) {
G4cout<<"!!!"<<severityWord[severity]<<"!!! "<<message
<<" ("<<(*it).second.nTimes<<" occurences)"<<G4endl;
}
}
if (severity==FATAL) {
G4cout<<"S T O P F O R C E D!"<<G4endl;
exit(1);
}
}
void musrErrorMessage::PrintErrorSummary() {
std::map<G4String,ErrorStruct>::iterator it;
G4cout<<"------ ERROR SUMMARY: ----------------------------------------------------------------"<<G4endl;
for (G4int i=0; i<4; i++) {
for ( it=ErrorMapping.begin() ; it != ErrorMapping.end(); it++ ) {
if ((*it).second.mesSeverity==i) {
G4cout<<severityWord[(*it).second.mesSeverity]<<" ("<<(*it).second.nTimes<<" times):"
<< (*it).first <<G4endl;
}
}
}
G4cout<<"-----------------------------------------------------------------------------------------"<<G4endl;
}

92
src/musrEventAction.cc Normal file
View File

@@ -0,0 +1,92 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrEventAction.hh"
#include "G4Event.hh"
#include "G4EventManager.hh"
#include "G4TrajectoryContainer.hh"
#include "G4Trajectory.hh"
#include "G4VVisManager.hh"
#include "G4ios.hh"
#include "G4TransportationManager.hh"
#include "G4FieldManager.hh"
#include "musrRootOutput.hh"
#include "musrErrorMessage.hh"
#include "musrSteppingAction.hh"
#include "F04GlobalField.hh"
#include "G4RunManager.hh" // needed just in order to save random number generator seed
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int musrEventAction::nHowOftenToPrintEvent=10000;
G4double musrEventAction::maximumRunTimeAllowed=85000;
musrEventAction::musrEventAction() {
pointer=this;
}
musrEventAction* musrEventAction::pointer=0;
musrEventAction* musrEventAction::GetInstance() {
return pointer;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrEventAction::~musrEventAction()
{
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrEventAction::BeginOfEventAction(const G4Event* evt) {
// G4cout<<"musrEventAction::BeginOfEventAction:"<<G4endl;
long thisEventNr = (long) (evt->GetEventID());
// if (thisEventNr == 44654) {trackingManager->SetVerboseLevel(2);}
musrSteppingAction::GetInstance()->DoAtTheBeginningOfEvent();
if (F04GlobalField::Exists()) {
F04GlobalField::getObject() -> CheckWhetherAnyNominalFieldValueNeedsToBeChanged(thisEventNr);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrEventAction::EndOfEventAction(const G4Event* evt) {
// cout << ":." << flush;
long thisEventNr = (long) evt->GetEventID();
// get number of stored trajectories
G4TrajectoryContainer* trajectoryContainer = evt->GetTrajectoryContainer();
G4int n_trajectories = 0;
if (trajectoryContainer) n_trajectories = trajectoryContainer->entries();
// G4cout << ">>> Event " << evt->GetEventID() << G4endl;
// periodic printing
//
// if (thisEventNr != 0 and thisEventNr%10000 == 0) {
if (thisEventNr == 0) timeOfRunStart=time(0);
if ((time(0)-timeOfRunStart)>maximumRunTimeAllowed) {
// Stop the execution of the run - the run took already too long time
char eMessage[200];
sprintf(eMessage,"musrEventAction::EndOfEventAction(): Run execution exceeded the allowed maximum time (maximum = %f sec) ==> RUN STOPPED",maximumRunTimeAllowed);
musrErrorMessage::GetInstance()->musrError(WARNING,eMessage,false);
G4RunManager* fRunManager = G4RunManager::GetRunManager();
fRunManager->AbortRun(true);
}
if (thisEventNr != 0 and thisEventNr%nHowOftenToPrintEvent == 0) {
time_t curr=time(0);
//char * ctime(const time_t * tp);
G4cout << ">>> Event " << evt->GetEventID() <<". Running already for "<<curr-timeOfRunStart<<" seconds. Present time: "<< ctime(&curr);
G4cout.flush();
// G4cout << " seed set to "<< CLHEP::HepRandom::getTheSeed();//<< G4endl;
}
// extract the trajectories and draw them
if (G4VVisManager::GetConcreteInstance()) {
for (G4int i=0; i<n_trajectories; i++)
{ G4Trajectory* trj = (G4Trajectory*)
((*(evt->GetTrajectoryContainer()))[i]);
trj->DrawTrajectory(1000);
}
}
}

111
src/musrMuFormation.cc Normal file
View File

@@ -0,0 +1,111 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium Formation according to yield.cc function (through GetYields method).
// Id : musrMuFormation.cc, v 1.4
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#include "musrMuFormation.hh"
using namespace std;
musrMuFormation::musrMuFormation(const G4String& name, G4ProcessType aType)
: G4VDiscreteProcess(name, aType){}
musrMuFormation::~musrMuFormation(){}
G4VParticleChange* musrMuFormation::PostStepDoIt(const G4Track& trackData,
const G4Step& aStep)
{ // Initialize ParticleChange (by setting all its members equal to
// the corresponding members in G4Track)
fParticleChange.Initialize(trackData);
G4Track theNewTrack;
if(CheckCondition(aStep))
{
GetDatas(&aStep);
G4Step theStep;
PrepareSecondary( trackData);
fParticleChange.AddSecondary(aSecondary);
fParticleChange.ProposeTrackStatus(fStopAndKill) ;
}
else
{
fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ;
}
return &fParticleChange;
}
G4bool musrMuFormation::CheckCondition(const G4Step& aStep)
{ // Decide when to call the MuFormation process - i.e. for muons going through the C foil.
G4bool condition=false;
p_name = aStep.GetTrack()->GetDefinition()->GetParticleName(); // particle name
//if(p_name == "mu+"&&aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName()=="log_CFoil")
std::string logVolName = aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
if(p_name == "mu+" && ((logVolName=="log_coulombCFoil")||(logVolName=="log_CFoil")))
{
condition=true;
}
return condition;
}
G4double musrMuFormation::GetMeanFreePath(const G4Track&,
G4double,
G4ForceCondition* condition)
{
*condition = Forced;
return DBL_MAX;
}
void musrMuFormation::GetDatas(const G4Step* aStep)
{ // Particle generation according to yield table
particleTable=G4ParticleTable::GetParticleTable();
rnd=G4UniformRand();
G4double E = aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy()/keV;
Gonin.GetYields(E,105.658369*1000,yvector); // Energy [keV], muon mass [keV/c2], yield table
G4String p_new = "Mu";
// Positive muon
if(p_name=="mu+")
{
if(rnd<yvector[0])
{
particle = particleTable->FindParticle(p_name) ;
}
else
{
particle = particleTable->FindParticle(p_new);
}
// Set the new dynamic particle DP
DP = new G4DynamicParticle(particle,
aStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection(),
aStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy());
// IMPORTANT : COPY THOSE DATA TO GET THE SAME PARTICLE PROPERTIES!!!
// SHOULD BE KEPT WHEN BUILDING A PARTICLE CHANGE
DP->SetProperTime( aStep->GetTrack()->GetDynamicParticle()->GetProperTime());
DP->SetPolarization(aStep->GetTrack()->GetDynamicParticle()->GetPolarization().x(),
aStep->GetTrack()->GetDynamicParticle()->GetPolarization().y(),
aStep->GetTrack()->GetDynamicParticle()->GetPolarization().z());
DP->SetPreAssignedDecayProperTime(aStep->GetTrack()->GetDynamicParticle()->GetPreAssignedDecayProperTime());
}
}
void musrMuFormation::PrepareSecondary(const G4Track& track)
{
if(p_name=="mu+")
{
aSecondary = new G4Track(DP,track.GetGlobalTime(),track.GetPosition());
}
}

86
src/musrMuScatter.cc Normal file
View File

@@ -0,0 +1,86 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium "Scattering"
// Id : musrMuScatter.cc, v 1.4
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
// Notes : Simplified model for Mu scattering. Spin effects have been excluded.
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#include "musrMuScatter.hh"
using namespace std;
musrMuScatter::musrMuScatter(const G4String& name,
G4ProcessType aType)
: G4VDiscreteProcess(name, aType){}
musrMuScatter:: ~musrMuScatter(){}
/*! - At the end of the step, the current volume is checked and if Muonium is in a solid
material (except for the carbon foil where it is generated), it is stopped immediately. */
G4VParticleChange* musrMuScatter::PostStepDoIt(const G4Track& trackData,
const G4Step& aStep)
{
fParticleChange.Initialize(trackData);
//! Tao - Get time information */
itime = trackData.GetProperTime();
gtime = trackData.GetGlobalTime();
ftime = trackData.GetDynamicParticle()->GetPreAssignedDecayProperTime();
deltatime = ftime - itime;
fParticleChange.ProposeGlobalTime(deltatime + itime -gtime);
/*! - Set position, momentum, energy and time of the particle change. */
fParticleChange.ProposePosition(trackData.GetPosition());
fParticleChange.ProposeMomentumDirection(trackData.GetMomentumDirection());
fParticleChange.ProposeEnergy(trackData.GetKineticEnergy());
fParticleChange.ProposeGlobalTime(gtime);
fParticleChange.ProposeProperTime(itime);
fParticleChange.ProposeTrackStatus(trackData.GetTrackStatus()) ;
/*! - Verify the condition of applying the process: if Mu is in a material
different than vacuum and carbon foil, then stop it directly. */
if( CheckCondition(aStep))
{
fParticleChange.ProposePosition(trackData.GetStep()->GetPreStepPoint()->GetPosition());
fParticleChange.ProposeTrackStatus(fStopButAlive) ;
}
/*! - Return the changed particle object. */
return &fParticleChange;
}
/*! - Muonium will be stopped as soon as it enters a material different than vacuum or C foil. */
G4bool musrMuScatter::CheckCondition(const G4Step& aStep)
{
G4bool condition = false;
p_name = aStep.GetTrack()->GetDefinition()->GetParticleName(); // particle name
if(p_name == "Mu" && aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetName()!="log_CFoil" &&
aStep.GetTrack()->GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()!="G4_Galactic")
{
condition=true;
}
return condition;
}
G4double musrMuScatter::GetMeanFreePath(const G4Track&,
G4double,
G4ForceCondition* condition)
{
*condition = Forced;
return DBL_MAX;
}
void musrMuScatter::PrepareSecondary(const G4Track& track)
{
aSecondary = new G4Track(DP,track.GetDynamicParticle()->GetPreAssignedDecayProperTime(),track.GetPosition());
}

112
src/musrMuonium.cc Normal file
View File

@@ -0,0 +1,112 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//
// ********************************************************************
// * License and Disclaimer *
// * *
// * The Geant4 software is copyright of the Copyright Holders of *
// * the Geant4 Collaboration. It is provided under the terms and *
// * conditions of the Geant4 Software License, included in the file *
// * LICENSE and available at http://cern.ch/geant4/license . These *
// * include a list of copyright holders. *
// * *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work make any representation or warranty, express or implied, *
// * regarding this software system or assume any liability for its *
// * use. Please see the license in the file LICENSE and URL above *
// * for the full disclaimer and the limitation of liability. *
// * *
// * This code implementation is the result of the scientific and *
// * technical work of the GEANT4 collaboration. *
// * By using, copying, modifying or distributing the software (or *
// * any work based on the software) you agree to acknowledge its *
// * use in resulting scientific publications, and indicate your *
// * acceptance of all terms of the Geant4 Software license. *
// ********************************************************************
//
//
// $Id: musrMuonium.cc,v 1.13 2007/03/15 06:53:58 kurasige Exp $
// GEANT4 tag $Name: geant4-09-00 $
//
//
// ----------------------------------------------------------------------
// GEANT 4 class implementation file
//
// History: first implementation, based on object model of
// 4th April 1996, G. Cosmo
// **********************************************************************
// New implementation as an utility class M. Asai, 26 July 2004
// ----------------------------------------------------------------------
#include "musrMuonium.hh"
#include "G4ParticleTable.hh"
#include "MuDecayChannel.hh"
#include "G4DecayTable.hh"
// ######################################################################
// ### MUONIUM ###
// ######################################################################
musrMuonium* musrMuonium::theInstance = 0;
musrMuonium* musrMuonium::Definition()
{
if (theInstance !=0) return theInstance;
const G4String name = "Mu";
// search in particle table]
G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* anInstance = pTable->FindParticle(name);
if (anInstance ==0)
{
// create particle
//
// Arguments for constructor are as follows
// name mass width charge
// 2*spin parity C-conjugation
// 2*Isospin 2*Isospin3 G-parity
// type lepton number baryon number PDG encoding
// stable lifetime decay table
// shortlived subType anti_encoding
anInstance = new G4ParticleDefinition(
name, 0.1056584*GeV, 2.99591e-16*MeV, 0.*eplus,
1, 0, 0,
0, 0, 0,
"lepton", -1, 0, -1313,
false, 2197.03*ns, NULL,
false, "mu"
);
// Bohr magnetron of Muonium - T. Shiroka
// The magnetic moment of Mu is the sum of those of mu+ and e- with
// the respective gyromagnetic ratio anomalies as coefficients
G4double muBmu = 0.5*eplus*hbar_Planck/(0.10565840*GeV/c_squared);
G4double muBel = -0.5*eplus*hbar_Planck/(0.51099906*MeV/c_squared);
G4double muB = 1.0011659208*muBmu + 1.0011596521859*muBel;
anInstance->SetPDGMagneticMoment( muB );
//create Decay Table
G4DecayTable* table = new G4DecayTable();
// create a decay channel
G4VDecayChannel* mode = new MuDecayChannel("Mu",1.00);
table->Insert(mode);
anInstance->SetDecayTable(table);
}
theInstance = reinterpret_cast<musrMuonium*>(anInstance);
return theInstance;
}
musrMuonium* musrMuonium::MuoniumDefinition()
{
return Definition();
}
musrMuonium* musrMuonium::Muonium()
{
return Definition();
}

91
src/musrParameters.cc Normal file
View File

@@ -0,0 +1,91 @@
#include "musrParameters.hh"
// #include "musrErrorMessage.hh" - musrErrorMessage class can not be used inside "musrParameters" constructor, because
// musrErrorMessage is crated later!
musrParameters::musrParameters(G4String steeringFileName)
{
pointerToParameters = this;
boolG4RegionRequested = false;
mySteeringFileName = steeringFileName;
// Read in the parameters, which have to be known before the detector construction is run
// (and therefore the parameters can not be read in in the musrDetectorConstruction.cc class).
FILE *fSteeringFile=fopen(steeringFileName.c_str(),"r");
if (fSteeringFile==NULL) {
G4cout<<"musrParameters::musrParameters: steeringFileName=\""<<steeringFileName
<<"\" not opened for some reason."<<G4endl;
G4cout << "S T O P F O R C E D" << G4endl;
exit(1);
}
G4cout<<"musrParameters::musrParameters: steeringFileName=\""<<steeringFileName<<"\" opened."<<G4endl;
char line[501];
while (!feof(fSteeringFile)) {
fgets(line,500,fSteeringFile);
if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
char tmpString0[100]="Unset";
sscanf(&line[0],"%s",tmpString0);
// First find out how many events will be generated (needs to be known at an early stage, if the
// field is to be set in steps):
if (strcmp(tmpString0,"/run/beamOn")==0) {
int nev;
sscanf(&line[0],"%*s %d", &nev);
musrParameters::nrOfEventsToBeGenerated = nev;
}
if (strncmp(tmpString0,"/gps/",5)==0) {
musrParameters::boolG4GeneralParticleSource = true;
G4cout<<"\n========================================================================"<<G4endl;
G4cout<<"musrParameters.cc: GPS (General Particle Source) requested in the macro."<<G4endl;
G4cout<<" GPS will be used instead of the primary generator action."<<G4endl;
}
// Now find some private parameters that need to be initialised at early stage
if ( (strcmp(tmpString0,"/musr/ignore")!=0)&&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
char tmpString1[100]="Unset", tmpString2[100]="Unset", tmpString3[100]="Unset";
sscanf(&line[0],"%*s %s %s %s",tmpString1,tmpString2,tmpString3);
// if (strcmp(tmpString1,"G4GeneralParticleSource")==0){
// if (strcmp(tmpString2,"true")==0){ musrParameters::boolG4GeneralParticleSource = true; }
// }
if (strcmp(tmpString1,"G4OpticalPhotons")==0){
if (strcmp(tmpString2,"true")==0){ musrParameters::boolG4OpticalPhotons = true; }
}
if (strcmp(tmpString1,"region")==0) {
boolG4RegionRequested = true;
}
if ( (strcmp(tmpString1,"construct")==0) && boolG4RegionRequested) {
G4cout<<"musrParameters.cc: User requests to construct a detector volume "<<tmpString3<<G4endl;
G4cout<<" after a previous G4Region definition."<<G4endl;
G4cout<<" Perhaps not a crutial problem, but the execution will be stopped"<<G4endl;
G4cout<<" anyway just to be on the safe side. Please correct "<<steeringFileName<<" file"<<G4endl;
G4cout<<" such that G4Region is called only after all detector volumes have been created."<<G4endl;
G4cout<<" S T O P F O R C E D!"<<G4endl;
exit(1);
}
}
}
fclose(fSteeringFile);
}
musrParameters::~musrParameters() {}
musrParameters* musrParameters::pointerToParameters=NULL;
musrParameters* musrParameters::GetInstance() {
return pointerToParameters;
}
G4String musrParameters::mySteeringFileName="Unset";
G4bool musrParameters::storeOnlyEventsWithHits=true;
G4int musrParameters::storeOnlyEventsWithHitInDetID=0;
G4double musrParameters::signalSeparationTime=100*nanosecond;
G4bool musrParameters::storeOnlyTheFirstTimeHit=false;
G4bool musrParameters::field_DecayWithSpin=false;
G4bool musrParameters::killAllPositrons=false;
G4bool musrParameters::killAllGammas=false;
G4bool musrParameters::killAllNeutrinos=true;
G4bool musrParameters::boolG4GeneralParticleSource=false;
G4bool musrParameters::boolG4OpticalPhotons=false;
//cks G4bool musrParameters::includeMuoniumProcesses =true; // TS
//G4bool musrParameters::boolG4GeneralParticleSource=true;
G4int musrParameters::nrOfEventsToBeGenerated=0;

444
src/musrPhysicsList.cc Normal file
View File

@@ -0,0 +1,444 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "globals.hh"
#include "G4ios.hh"
#include "musrPhysicsList.hh"
#include "G4VPhysicsConstructor.hh"
#include "G4ProcessManager.hh"
#include "G4ParticleTypes.hh"
#include "G4MuonDecayChannel.hh"
#include "G4DecayTable.hh"
//cks Added to have Geant default muon decay with spin
#include "G4MuonDecayChannelWithSpin.hh"
#include "G4MuonRadiativeDecayChannelWithSpin.hh"
#include "G4RadioactiveDecay.hh"
#include "G4IonConstructor.hh"
//TS Classes which account for Muonium as "particle" and its spin
#include "musrMuonium.hh"
#include "MuDecayChannel.hh"
#include "MuDecayChannelWithSpin.hh"
//
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrPhysicsList::musrPhysicsList(): G4VUserPhysicsList()
{
defaultCutValue = 0.1*mm;
SetVerboseLevel(0);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrPhysicsList::~musrPhysicsList()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructParticle()
{
ConstructBosons();
ConstructLeptons();
ConstructMesons();
ConstructBaryons();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructBosons()
{
// pseudo-particles
G4Geantino::GeantinoDefinition();
G4ChargedGeantino::ChargedGeantinoDefinition();
// gamma
G4Gamma::GammaDefinition();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructLeptons()
{
// leptons
// e+/-
G4Electron::ElectronDefinition();
G4Positron::PositronDefinition();
// mu+/-
G4MuonPlus::MuonPlusDefinition();
G4MuonMinus::MuonMinusDefinition();
//cks
// G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
// MuonPlusDecayTable -> Insert(new musrMuonDecayChannel("mu+",1.00));
// G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
//csk
//
// Muonium - TS
musrMuonium::MuoniumDefinition();
//
// nu_e
G4NeutrinoE::NeutrinoEDefinition();
G4AntiNeutrinoE::AntiNeutrinoEDefinition();
// nu_mu
G4NeutrinoMu::NeutrinoMuDefinition();
G4AntiNeutrinoMu::AntiNeutrinoMuDefinition();
//cks: Trial to use Geant4 muon decay with spin
G4DecayTable* MuonPlusDecayTable = new G4DecayTable();
MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",1.00));
// MuonPlusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu+",0.986));
// MuonPlusDecayTable -> Insert(new G4MuonRadiativeDecayChannelWithSpin("mu+",0.014));
G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable);
//
// G4DecayTable* MuonMinusDecayTable = new G4DecayTable();
// MuonMinusDecayTable -> Insert(new G4MuonDecayChannelWithSpin("mu-",1.00));
// G4MuonMinus::MuonMinusDefinition() -> SetDecayTable(MuonMinusDecayTable);
//csk
//
//TS: Using the muonium decay with and without spin
G4DecayTable* MuoniumDecayTable = new G4DecayTable();
MuoniumDecayTable -> Insert(new MuDecayChannel("Mu",0.50));
MuoniumDecayTable -> Insert(new MuDecayChannelWithSpin("Mu",0.5));
musrMuonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable);
//MuoniumDecayTable ->DumpInfo(); // Info on muonium decay channels
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructMesons()
{
// mesons
// light mesons
G4PionPlus::PionPlusDefinition();
G4PionMinus::PionMinusDefinition();
G4PionZero::PionZeroDefinition();
G4Eta::EtaDefinition();
G4EtaPrime::EtaPrimeDefinition();
G4KaonPlus::KaonPlusDefinition();
G4KaonMinus::KaonMinusDefinition();
G4KaonZero::KaonZeroDefinition();
G4AntiKaonZero::AntiKaonZeroDefinition();
G4KaonZeroLong::KaonZeroLongDefinition();
G4KaonZeroShort::KaonZeroShortDefinition();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructBaryons()
{
// baryons
G4Proton::ProtonDefinition();
G4AntiProton::AntiProtonDefinition();
G4Neutron::NeutronDefinition();
G4AntiNeutron::AntiNeutronDefinition();
// ions
G4IonConstructor iConstructor;
iConstructor.ConstructParticle();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructProcess()
{
AddTransportation();
ConstructEM();
ConstructGeneral();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4ComptonScattering.hh"
#include "G4GammaConversion.hh"
#include "G4PhotoElectricEffect.hh"
#include "G4MultipleScattering.hh"
#include "G4eIonisation.hh"
#include "G4eBremsstrahlung.hh"
#include "G4eplusAnnihilation.hh"
#include "G4MuIonisation.hh"
#include "G4MuBremsstrahlung.hh"
#include "G4MuPairProduction.hh"
#include "G4hIonisation.hh"
#include "G4UserSpecialCuts.hh"
//#include "musrAtRestSpinRotation.hh"
// For low energy physics processes:
#include "G4LowEnergyCompton.hh"
//#include "G4LowEnergyPolarizedCompton.hh"
#include "G4LowEnergyGammaConversion.hh"
#include "G4LowEnergyPhotoElectric.hh"
#include "G4LowEnergyRayleigh.hh"
#include "G4LowEnergyBremsstrahlung.hh"
#include "G4LowEnergyIonisation.hh"
#include "G4hLowEnergyIonisation.hh"
// For Penelope processes:
#include "G4PenelopeCompton.hh"
#include "G4PenelopeGammaConversion.hh"
#include "G4PenelopePhotoElectric.hh"
#include "G4PenelopeRayleigh.hh"
#include "G4PenelopeIonisation.hh"
#include "G4PenelopeBremsstrahlung.hh"
#include "G4PenelopeAnnihilation.hh"
// For Coulomb scattering instead of multiple scattering
#include "G4CoulombScattering.hh"
#include "G4CoulombScatteringModel.hh"
// For Muonium formation in the Carbon foil
#include "musrMuFormation.hh" // includes the yield function Y = Y(E).
// For a simple Muonium "scattering" when Mu hits solid materials
#include "musrMuScatter.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ConstructEM()
{
// cks 2008.08.22. - Adding the possibility to define the processes from the steering file:
char charSteeringFileName[1000]; strcpy(charSteeringFileName,(musrParameters::mySteeringFileName).c_str());
FILE *fSteeringFile=fopen(charSteeringFileName,"r");
if (fSteeringFile==NULL) {
sprintf(eMessage,"musrPhysicsList::ConstructEM(): Failed to open macro file \"%s\" .",charSteeringFileName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
else {G4cout<<"musrPhysicsList: The Physics processes are being defined:"<<G4endl;}
char line[501];
while (!feof(fSteeringFile)) {
fgets(line,500,fSteeringFile);
if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
char tmpString0[100]="Unset", tmpString1[100]="Unset",tmpString2[100]="Unset";
sscanf(&line[0],"%s %s %s",tmpString0,tmpString1,tmpString2);
if ( (strcmp(tmpString0,"/musr/ignore")!=0)&&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
if (strcmp(tmpString1,"process")!=0) continue;
if ((strcmp(tmpString2,"addProcess")==0)||(strcmp(tmpString2,"addDiscreteProcess")==0)) {
char charParticleName[100], charProcessName[100];
sscanf(&line[0],"%*s %*s %s %s %s",tmpString2,charParticleName,charProcessName);
G4cout<<"musrPhysicsList: Defining process "<<charProcessName<<" for "<<charParticleName<<G4endl;
G4String stringProcessName = charProcessName;
G4String stringParticleName = charParticleName;
G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable() -> FindParticle(stringParticleName);
// G4cout<<"particleDefinition of "<<stringParticleName<<" = "<<particleDefinition<<G4endl;
if (particleDefinition==NULL) {
sprintf(eMessage,"musrPhysicsList: Partile \"%s\" not found in G4ParticleTable when trying to assign process \"%s\".",
charParticleName,charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
G4ProcessManager* pManager = particleDefinition->GetProcessManager();
if (strcmp(tmpString2,"addDiscreteProcess")==0) {
if (stringProcessName=="G4PhotoElectricEffect") pManager->AddDiscreteProcess(new G4PhotoElectricEffect);
else if (stringProcessName=="G4ComptonScattering") pManager->AddDiscreteProcess(new G4ComptonScattering);
else if (stringProcessName=="G4GammaConversion") pManager->AddDiscreteProcess(new G4GammaConversion);
else if (stringProcessName=="G4PenelopePhotoElectric") pManager->AddDiscreteProcess(new G4PenelopePhotoElectric);
else if (stringProcessName=="G4PenelopeCompton") pManager->AddDiscreteProcess(new G4PenelopeCompton);
else if (stringProcessName=="G4PenelopeGammaConversion") pManager->AddDiscreteProcess(new G4PenelopeGammaConversion);
else if (stringProcessName=="G4PenelopeRayleigh") pManager->AddDiscreteProcess(new G4PenelopeRayleigh);
else if (stringProcessName=="G4LowEnergyPhotoElectric") pManager->AddDiscreteProcess(new G4LowEnergyPhotoElectric);
else if (stringProcessName=="G4LowEnergyCompton") pManager->AddDiscreteProcess(new G4LowEnergyCompton);
else if (stringProcessName=="G4LowEnergyGammaConversion") pManager->AddDiscreteProcess(new G4LowEnergyGammaConversion);
else if (stringProcessName=="G4LowEnergyRayleigh") pManager->AddDiscreteProcess(new G4LowEnergyRayleigh);
else if (stringProcessName=="G4CoulombScattering") pManager->AddDiscreteProcess(new G4CoulombScattering);
else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addDiscreteProcess. It can be easily added.",
charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
}
else if (strcmp(tmpString2,"addProcess")==0) {
G4int nr1, nr2, nr3;
char charRegion1[100]="", charRegion2[100]="", charRegion3[100]="", charControlString[10]="";
sscanf(&line[0],"%*s %*s %*s %*s %*s %d %d %d %s %s %s %s",&nr1,&nr2,&nr3,charRegion1,charRegion2,charRegion3,charControlString);
if (stringProcessName=="G4MultipleScattering") pManager->AddProcess(new G4MultipleScattering,nr1,nr2,nr3);
else if (stringProcessName=="G4eIonisation") pManager->AddProcess(new G4eIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4eBremsstrahlung") pManager->AddProcess(new G4eBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4eplusAnnihilation") pManager->AddProcess(new G4eplusAnnihilation,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeIonisation") pManager->AddProcess(new G4PenelopeIonisation,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeBremsstrahlung") pManager->AddProcess(new G4PenelopeBremsstrahlung,nr1,nr2,nr3);
else if (stringProcessName=="G4PenelopeAnnihilation") pManager->AddProcess(new G4PenelopeAnnihilation,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=="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);
// 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);
// cks: musrMuScatter could be uncommented here, but testing is needed, because Toni has some strange comments
// in his original "musrPhysicsList.cc about implementing musrMuScatter.
// else if (stringProcessName=="musrMuScatter") pManager->AddProcess(new musrMuScatter,nr1,nr2,nr3);
else if (stringProcessName=="MultipleAndCoulombScattering") {
G4MultipleScattering* multScat = new G4MultipleScattering();
// G4CoulombScattering* coulScat = new G4CoulombScattering();
G4CoulombScatteringModel* coulScatModel = new G4CoulombScatteringModel();
if (strcmp(charRegion1,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion1,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion1<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
// multScat->AddEmModel(0,multScat,regionForCoulomb);
}
if (strcmp(charRegion2,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion2,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion2<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
}
if (strcmp(charRegion3,"")!=0) {
G4Region* regionForCoulomb = FindG4Region(charRegion3,line);
G4cout<<" Adding Coulomb scattering model to multiple scattering model for region "<<charRegion3<<G4endl;
multScat->AddEmModel(0,coulScatModel,regionForCoulomb);
}
if (strcmp(charControlString,"")!=0) {
G4cout<<"More than 3 regions requested for Coulomb Scattering, but presently only up to 3 such regions are supported."<<G4endl;
G4cout<<"Please extend the number of supported regions in musrPhysicsList.cc to higher number."<<G4endl;
G4cout<<"The extention of the code to larger number of regions is not very difficult."<<G4endl;
G4cout<<" S T O P F O R C E D"<<G4endl;
exit(1);
}
pManager->AddProcess(multScat,nr1,nr2,nr3);
}
else {
sprintf(eMessage,"musrPhysicsList: Process \"%s\" is not implemented in musrPhysicsList.cc for addProcess. It can be easily added.",
charProcessName);
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
}
}
else ReportProblemWithProcessDefinition(line);
}
}
fclose(fSteeringFile);
G4cout<<"\n\n\n\n"<<G4endl;
// csk 2008.08.22.
//del G4String myTypeOfProcesses = musrParameters::GetInstance()->GetMyTypeOfProcesses();
//del G4cout<<"musrPhysicsList::ConstructEM(): myTypeOfProcesses="<<myTypeOfProcesses<<G4endl;
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
G4String particleName = particle->GetParticleName();
if ((particleName == "gamma")||(particleName == "e-")||(particleName == "e+")) {
// do nothing
}
else if ((particleName=="mu+")||(particleName=="mu-")) { //muon
G4DecayWithSpin* theDecayProcess = new G4DecayWithSpin();
// theDecayProcess->SetVerboseLevel(2);
pmanager->AddProcess(theDecayProcess);
pmanager ->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep);
}
else if (particleName=="Mu") {
// TS:
// Muonium "scattering" Kamil: the following 3 lines could be replaced by reading the musrMuScatter
// process through the steering file
G4VProcess* aMuScatt = new musrMuScatter();
pmanager->AddProcess(aMuScatt);
pmanager->SetProcessOrdering(aMuScatt, idxPostStep, 1);
//
G4Decay* theDecayProcess = new G4Decay();
//musrDecayWithSpin* theDecayProcess = new musrDecayWithSpin();
pmanager->AddProcess(theDecayProcess);
pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest);
pmanager->SetProcessOrdering(theDecayProcess, idxPostStep);
}
else if ((!particle->IsShortLived()) &&
(particle->GetPDGCharge() != 0.0) &&
(particle->GetParticleName() != "chargedgeantino")) {
//all others charged particles except geantino
pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
// if (myTypeOfProcesses=="highenergy") {
// pmanager->AddProcess(new G4hIonisation, -1, 2,2);
// }
// else {
pmanager->AddProcess(new G4hLowEnergyIonisation, -1, 2,2);
// }
///pmanager->AddProcess(new G4UserSpecialCuts, -1,-1,3);
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Decay.hh"
void musrPhysicsList::ConstructGeneral() {
if (musrParameters::boolG4GeneralParticleSource) {
G4RadioactiveDecay* theRadioactiveDecay = new G4RadioactiveDecay();
G4GenericIon* ion = G4GenericIon::GenericIon();
theParticleIterator->reset();
while( (*theParticleIterator)() ){
G4ParticleDefinition* particle = theParticleIterator->value();
G4ProcessManager* pmanager = particle->GetProcessManager();
if (particle == ion) {
pmanager->AddProcess(theRadioactiveDecay, 0, -1, 3);
}
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "G4Region.hh"
#include "G4RegionStore.hh"
#include "G4ProductionCuts.hh"
void musrPhysicsList::SetCuts()
{
//G4VUserPhysicsList::SetCutsWithDefault method sets
//the default cut value for all particle types
//
SetCutsWithDefault();
if (verboseLevel>0) DumpCutValuesTable();
DumpCutValuesTable();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrPhysicsList::ReportProblemWithProcessDefinition(char myString[501]) {
G4cout<<"\nE R R O R in musrPhysicsList.cc: "
<<"Unknown keyword requested in the steering (*.mac) file :"<<G4endl;
G4cout<<" "<<myString<<G4endl;
G4cout<<"S T O P F O R C E D!"<<G4endl;
exit(1);
}
G4Region* musrPhysicsList::FindG4Region(G4String regionName, char* lineOfSteeringFile) {
G4Region* myRegion = G4RegionStore::GetInstance()->GetRegion(regionName,false);
if( myRegion != NULL ) { // G4Region found
return myRegion;
}
else { // G4Region not found
G4cout<<"musrPhysicsList: G4Region "<<regionName<<" not found."<<G4endl;
G4cout<<" The critical command line of the steering file is:"<<G4endl;
G4cout<<" "<<lineOfSteeringFile<<G4endl;
G4cout<<" S T O P F O R C E D"<<G4endl;
exit(1);
}
}

427
src/musrPrimaryGeneratorAction.cc Normal file
View File

@@ -0,0 +1,427 @@
#include "musrPrimaryGeneratorAction.hh"
#include "musrDetectorConstruction.hh"
#include "musrPrimaryGeneratorMessenger.hh"
#include "musrParameters.hh"
#include "G4Event.hh"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "Randomize.hh"
#include "G4ios.hh"
#include "G4UnitsTable.hh"
#include "globals.hh"
#include "G4Gamma.hh"
#include "G4ThreeVector.hh"
#include "G4RunManager.hh"
#include "time.h"
#include <iomanip>
#include "musrRootOutput.hh" //cks for storing some info in the Root output file
#include "musrErrorMessage.hh"
G4bool musrPrimaryGeneratorAction::setRandomNrSeedAccordingEventNr=0;
G4bool musrPrimaryGeneratorAction::setRandomNrSeedFromFile=0;
G4bool musrPrimaryGeneratorAction::setRandomNrSeedFromFile_RNDM=0;
G4int musrPrimaryGeneratorAction::nRndmEventToSaveSeeds=-2;
std::vector<int> * musrPrimaryGeneratorAction::pointerToSeedVector=NULL;
std::vector<int> * musrPrimaryGeneratorAction::GetPointerToSeedVector() {
return pointerToSeedVector;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
musrDetectorConstruction* musrDC)
:musrDetector(musrDC), x0(0), y0(0), z0(-10*cm), xSigma(0), ySigma(0), zSigma(0),
rMaxAllowed(1e10*mm), zMinAllowed(-1e10*mm), zMaxAllowed(1e10*mm),
p0(0), pSigma(0), pMinAllowed(0), pMaxAllowed(1e10*mm),
xangle0(0), yangle0(0), xangleSigma(0), yangleSigma(0), pitch(0),
UnpolarisedMuonBeam(false), TransversalyUnpolarisedMuonBeam(false), xPolarisIni(1.), yPolarisIni(0.), zPolarisIni(0.),
polarisFraction(1.),
muonDecayTimeMin(-1), muonDecayTimeMax(-1), muonMeanLife(2197.03*ns),
takeMuonsFromTurtleFile(false), z0_InitialTurtle(0),
numberOfGeneratedEvents(0),
turtleMomentumBite(false), turtleMomentumP0(0.), turtleSmearingFactor(0.)
//, firstCall(true)
{
//create a messenger for this class
gunMessenger = new musrPrimaryGeneratorMessenger(this);
// create a vector for storing the event numbers as seeds for the random number generator
pointerToSeedVector = new std::vector<int>;
// default particle kinematic
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* muonParticle= particleTable->FindParticle("mu+");
// cks Implement also alpha and proton particles for the simulation of Juan Pablo Urrego
alphaParticle= particleTable->FindParticle("alpha");
protonParticle= particleTable->FindParticle("proton");
// csk
G4int n_particle = 1;
if (musrParameters::boolG4GeneralParticleSource) {
G4cout<<"musrPrimaryGeneratorAction: G4GeneralParticleSource is going to be initialised"<<G4endl;
particleSource = new G4GeneralParticleSource ();
}
else {
G4cout<<"musrPrimaryGeneratorAction: G4ParticleGun is going to be initialised"<<G4endl;
particleGun = new G4ParticleGun(n_particle);
particleGun->SetParticleDefinition(muonParticle);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
musrPrimaryGeneratorAction::~musrPrimaryGeneratorAction()
{
if (musrParameters::boolG4GeneralParticleSource) {delete particleSource;}
else {delete particleGun;}
delete gunMessenger;
if (takeMuonsFromTurtleFile) {fclose(fTurtleFile);}
G4cout<<"musrPrimaryGeneratorAction: Number of Generated Events = "<<numberOfGeneratedEvents<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
// This function is called at the begining of event.
// Clear Root variables
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
myRootOutput->ClearAllRootVariables(); // Note that musrPrimaryGeneratorAction::GeneratePrimaries
// is called before the musrEventAction::BeginOfEventAction.
// Therefore "ClearAllRootVariables" is called already here
// (before the "SetInitialMuonParameters".
// Set or read the seeds of random number generator
boolPrintInfoAboutGeneratedParticles=false;
SetOrReadTheRandomNumberSeeds(anEvent->GetEventID());
// If radioactive source is used, use G4GeneralParticleSource :
if (musrParameters::boolG4GeneralParticleSource) {
particleSource->GeneratePrimaryVertex(anEvent);
return;
}
G4double x, y, z;
G4double p;
G4double xangle, yangle;
if (takeMuonsFromTurtleFile) {
char line[501];
G4int checkNrOfCounts=0;
do {
float xTmp, yTmp, xAngleTmp, yAngleTmp, pTmp;
float dummy1, dummy2;
int Ztmp=-1, Atmp=-1;
fgets(line,500,fTurtleFile);
if (feof(fTurtleFile)) {
rewind(fTurtleFile);
G4cout<<"End of TurtleFile, lets start from the beginning (numberOfGeneratedEvents = "<<numberOfGeneratedEvents<<")"<<G4endl;
fgets(line,500,fTurtleFile);
}
numberOfGeneratedEvents++;
sscanf(&line[0],"%g %g %g %g %g %g %g %d %d",&xTmp,&xAngleTmp,&yTmp,&yAngleTmp,&pTmp,&dummy1,&dummy2,&Ztmp,&Atmp);
if (boolPrintInfoAboutGeneratedParticles) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: Turtle input for this event: "
<<xTmp<<", "<<xAngleTmp<<" "<<yTmp<<" "<<yAngleTmp<<" "<< pTmp<<G4endl;
}
//cks Implement also alpha and proton particles for the simulation of Juan Pablo Urrego
if ((Ztmp==1)&&(Atmp==1)) {particleGun->SetParticleDefinition(protonParticle);}// G4cout<<"proton"<<G4endl;}
else if ((Ztmp==2)&&(Atmp==4)) {particleGun->SetParticleDefinition(alphaParticle);}// G4cout<<"alpha particle"<<G4endl;}
else if ((Ztmp==-1)&&(Atmp==-1)) {;}
else {
G4cout<<"musrPrimaryGeneratorAction: Unknown particle requested in the TURTLE input file: Z="
<<Ztmp<<", A="<<Atmp<<G4endl<<"S T O P F O R C E D" << G4endl;
G4cout<<xTmp<<", "<<xAngleTmp<<", "<<yTmp<<", "<<yAngleTmp<<", "<<pTmp<<", "<<dummy1<<", "<<dummy2<<", "<<Ztmp<<", "<<Atmp<<G4endl;
exit(1);
}
//csk
xangle = xAngleTmp*mrad;
yangle = yAngleTmp*mrad;
x = xTmp*cm + (z0-z0_InitialTurtle)*tan(xangle) ; // usually z0 is negative
y = yTmp*cm + (z0-z0_InitialTurtle)*tan(yangle) ; // z0_InitialTurtle is the z0 at whith the turtle file was generated.
p = pTmp*GeV;
// add some offset, if requested:
x = x + x0;
y = y + y0;
// add some beam tilt, if requested:
xangle = xangle + xangle0;
yangle = yangle + yangle0;
// add some beam pitch, if requested:
if (pitch!=0) {
xangle += - pitch * (x-x0);
yangle += - pitch * (y-y0);
}
// add/remove some momentum smearing, if requested
if (turtleMomentumBite) {
p = turtleMomentumP0 - (turtleMomentumP0-p)*turtleSmearingFactor;
}
checkNrOfCounts++;
if (checkNrOfCounts>1000) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: Too strict requirements on the r position!"<<G4endl;
}
} while( (x*x+y*y)>(rMaxAllowed*rMaxAllowed) );
z=z0;
// G4cout<<"x,y,z=("<<x/mm<<","<<y/mm<<","<<z/mm<<"), angles="<<xangle/mrad<<","<<yangle/mrad<<" p="<<p/MeV<<G4endl;
}
else { // Generate the starting position of the muon by random
// rMaxAllowed ... maximal radius, within which the muon can be generated
// x0, y0, z0 ... central point around which the muons are generated
// xSigma, ySigma, zSigma ... sigma of the (gaussian) distributions of the beam
// x, y, z ... actual initial position of the generated muon
G4int checkNrOfCounts=0;
numberOfGeneratedEvents++;
do {
if (xSigma>0) {x = G4RandGauss::shoot(x0,xSigma);} // Gaussian distribution
else if (xSigma<0) {x = x0 + xSigma*(G4UniformRand()*2.-1.);} // Uniform step distribution
else { x = x0;} // Point-like
if (ySigma>0) {y = G4RandGauss::shoot(y0,ySigma);}
else if (ySigma<0) {y = y0 + ySigma*(G4UniformRand()*2.-1.);}
else {y = y0;}
if (zSigma>0) {z = G4RandGauss::shoot(z0,zSigma);}
else if (zSigma<0) {z = z0 + zSigma*(G4UniformRand()*2.-1.);}
else {z = z0;}
checkNrOfCounts++;
if (checkNrOfCounts>1000) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: Too strict requirements on the r or z position!"<<G4endl;
}
} while( ((x*x+y*y)>(rMaxAllowed*rMaxAllowed))||(z>zMaxAllowed)||(z<zMinAllowed) );
// The generated muon has to stay
// within some well defined region,
// e.g. within the beampipe
// Now generate the momentum
checkNrOfCounts=0;
do {
if (pSigma>0) {p = G4RandGauss::shoot(p0,pSigma);}
else {p=p0;}
checkNrOfCounts++;
if (checkNrOfCounts>1000) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: Too strict requirements on the momentum!"<<G4endl;
}
} while ( (p>pMaxAllowed)||(p<pMinAllowed) );
// Add some initial angle (px and py component of the momentum)
if (xangleSigma>0) { xangle = G4RandGauss::shoot(xangle0,xangleSigma); }
else { xangle = xangle0; }
// Add the beam tilt, which depends on the distance from the beam centre.
// if (xSigma>0) {xangle += - pitch * (x-x0)/xSigma; }
if (pitch!=0) {xangle += - pitch * (x-x0); }
if (yangleSigma>0) { yangle = G4RandGauss::shoot(yangle0,yangleSigma); }
else { yangle = yangle0; }
// Add the beam tilt, which depends on the distance from the beam centre.
// if (ySigma>0) {yangle += - pitch * (y-y0)/ySigma; }
if (pitch!=0) {yangle += - pitch * (y-y0); }
} // end of the part specific for the muons generated by random rather then from TURTLE
// 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);
// Assign spin
G4double xpolaris=0, ypolaris=0, zpolaris=0;
if (UnpolarisedMuonBeam) {
// for genarating random numbers on the sphere see http://mathworld.wolfram.com/SpherePointPicking.html
G4double thetaTMP=pi/2;
if(!TransversalyUnpolarisedMuonBeam) thetaTMP = acos(2. * G4UniformRand()-1);
G4double phiTMP = 2. * pi * G4UniformRand();
xpolaris = std::sin(thetaTMP) * std::cos(phiTMP);;
ypolaris = std::sin(thetaTMP) * std::sin(phiTMP);
zpolaris = std::cos(thetaTMP);
}
else {
if (G4UniformRand()>((1.-polarisFraction)/2.)) {
xpolaris = xPolarisIni; ypolaris = yPolarisIni; zpolaris = zPolarisIni;
// G4cout<<"spin up"<<G4endl;
}
else {
xpolaris = -xPolarisIni; ypolaris = -yPolarisIni; zpolaris = -zPolarisIni;
// G4cout<<"spin down"<<G4endl;
}
}
particleGun->SetParticlePosition(G4ThreeVector(x,y,z));
G4double particle_mass = particleGun->GetParticleDefinition()->GetPDGMass();
G4double particleEnergy = std::sqrt(p*p+particle_mass*particle_mass)-particle_mass;
particleGun->SetParticleEnergy(particleEnergy);
particleGun->SetParticleMomentumDirection(G4ThreeVector(px,py,pz));
particleGun->SetParticlePolarization(G4ThreeVector(xpolaris,ypolaris,zpolaris));
particleGun->GeneratePrimaryVertex(anEvent);
// G4cout<<"musrPrimaryGeneratorAction: Parameters:"<<G4endl;
// G4cout<<" x0,y0,z0="<<x0/mm<<","<<y0/mm<<","<<z0/mm<<" Sigma="<<xSigma/mm<<","<<ySigma/mm<<","<<zSigma/mm<<G4endl;
// G4cout<<" rMaxAllowed="<<rMaxAllowed/mm<<" zMaxAllowed="<<zMaxAllowed/mm<<" zMinAllowed="<<zMinAllowed/mm<<G4endl;
// G4cout<<" p0="<<p0/MeV<<" pSigma="<<pSigma/MeV
// <<" pMinAllowed="<<pMinAllowed/MeV<<" pMaxAllowed=<<"<<pMaxAllowed/MeV<<G4endl;
// G4cout<<" angle0="<<xangle0/deg<<","<<yangle0/deg<<",nic"
// <<" Sigma="<<xangleSigma/deg<<","<<yangleSigma/deg<<",nic"<<G4endl;
// G4cout<<" pitch="<<pitch/deg<<G4endl;
//
// G4cout<<"musrPrimaryGeneratorAction: Generated muon:"<<G4endl;
// G4cout<<" x,y,z="<<x/mm<<","<<y/mm<<","<<z/mm<<" angle="<<xangle/deg<<","<< yangle/deg<<",nic"<<G4endl;
// G4cout<<" p="<<px/MeV<<","<<py/MeV<<","<<pz/MeV<<" E="<< (particleGun->GetParticleEnergy())/MeV<<G4endl;
// G4cout<<" polarisation="<<xpolaris<<","<<ypolaris<<","<<zpolaris<<G4endl;
// if requested by "/gun/decaytimelimits", set the decay time of the muon such that it is within
// the required time window. Otherwise the decay time is set internally by Geant.
if (muonDecayTimeMax>0.) {
// G4cout<<"muonDecayTimeMin="<<muonDecayTimeMin/ns<<" ns , muonDecayTimeMax="<<muonDecayTimeMax/ns
// <<" ns , muonMeanLife="<<muonMeanLife/ns<<" ns."<<G4endl;
// find the primary muon
G4PrimaryParticle* generatedMuon = anEvent->GetPrimaryVertex(0)->GetPrimary(0);
// G4double decayLowerLimit = 1-exp(-muonDecayTimeMin/muonMeanLife);
// G4double decayUpperLimit = 1-exp(-muonDecayTimeMax/muonMeanLife);
// G4double randomVal = G4UniformRand()*(decayUpperLimit-decayLowerLimit) + decayLowerLimit;
// G4double decaytime = -muonMeanLife*log(1-randomVal);
//
// The following code is numerically more stable compared to the commented lines above:
G4double expMin = exp(-muonDecayTimeMin/muonMeanLife);
G4double expMax = exp(-muonDecayTimeMax/muonMeanLife);
G4double decaytime = -muonMeanLife * log(G4UniformRand()*(expMax-expMin)+expMin);
// G4cout<<"decaytime="<<decaytime/ns<<"ns."<< G4endl;
generatedMuon->SetProperTime(decaytime);
}
// Save variables into ROOT output file:
myRootOutput->SetInitialMuonParameters(x,y,z,px,py,pz,xpolaris,ypolaris,zpolaris);
myRootOutput->StoreGeantParameter(7,float(numberOfGeneratedEvents));
if (boolPrintInfoAboutGeneratedParticles) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: x="<<x<<", y="<<y<<", z="<<z<<G4endl;
G4cout<<" px="<<px<<", py="<<py<<", pz="<<pz<<", xpolaris="<<xpolaris<<", ypolaris="<<ypolaris<<", zpolaris="<<zpolaris<<G4endl;
G4cout<<" numberOfGeneratedEvents="<<numberOfGeneratedEvents<<G4endl;
G4cout<<" ------------------------------------"<<G4endl;
}
}
//===============================================================================
void musrPrimaryGeneratorAction::SetInitialMuonPolariz(G4ThreeVector vIniPol)
{
G4double magnitude=vIniPol.mag();
if(magnitude<0.00000001) {
G4cout<< "Unpolarised initial muons"<<G4endl;
UnpolarisedMuonBeam=true;
if ((magnitude<0.0000000085)&&(magnitude>0.0000000075)) {
G4cout<< "Transversaly unpolarised initial muons"<<G4endl;
TransversalyUnpolarisedMuonBeam=true;
}
}
else {
xPolarisIni=vIniPol(0)/magnitude;
yPolarisIni=vIniPol(1)/magnitude;
zPolarisIni=vIniPol(2)/magnitude;
G4cout<< "Initial Muon Polarisation set to ("<<xPolarisIni<<","<<yPolarisIni<<","<<zPolarisIni<<")"<<G4endl;
}
}
//===============================================================================
void musrPrimaryGeneratorAction::SetMuonDecayTimeLimits(G4ThreeVector decayTimeLimits) {
muonDecayTimeMin = decayTimeLimits[0];
muonDecayTimeMax = decayTimeLimits[1];
muonMeanLife = decayTimeLimits[2];
// store the muon decay time parameters to the Root output
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
myRootOutput->StoreGeantParameter(2,muonDecayTimeMin/microsecond);
myRootOutput->StoreGeantParameter(3,muonDecayTimeMax/microsecond);
myRootOutput->StoreGeantParameter(4,muonMeanLife/microsecond);
}
//===============================================================================
void musrPrimaryGeneratorAction::SetTurtleInput(G4String turtleFileName) {
takeMuonsFromTurtleFile = true;
fTurtleFile = fopen(turtleFileName.c_str(),"r");
if (fTurtleFile==NULL) {
G4cout << "E R R O R : Failed to open TURTLE input file \"" << turtleFileName
<<"\"."<< G4endl;
G4cout << "S T O P F O R C E D" << G4endl;
exit(1);
}
else {G4cout << "Turtle input file \"" << turtleFileName <<"\" opened."<< G4endl;}
}
//===============================================================================
void musrPrimaryGeneratorAction::SetTurtleInputFileToEventNo(G4int lineNumberOfTurtleFile) {
if (fTurtleFile==NULL) {
G4cout << "musrPrimaryGeneratorAction::SetTurtleInputFileToEventNo:"
<<" TURTLE input file not found - line number can not be set."<<G4endl;
}
else {
char line[501];
for (Int_t i=0; i<lineNumberOfTurtleFile; i++) {
if (feof(fTurtleFile)) rewind(fTurtleFile);
fgets(line,500,fTurtleFile);
}
G4cout << "musrPrimaryGeneratorAction::SetTurtleInputFileToEventNo: Turtle input file will start at line no.:"
<< lineNumberOfTurtleFile <<G4endl;
}
}
//===============================================================================
void musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds(G4int eventID) {
if (eventID == nRndmEventToSaveSeeds) {
G4cout<<"musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds: S A V I N G R A N D O M N O. S E E D S"<<G4endl;
G4cout<<" (for even nr. "<<eventID<<")"<<G4endl;
G4RunManager::GetRunManager()->rndmSaveThisEvent();
boolPrintInfoAboutGeneratedParticles = true;
}
if (eventID == 0) {
if (setRandomNrSeedFromFile_RNDM) {
G4cout<<"musrPrimaryGeneratorAction::SetOrReadTheRandomNumberSeeds: Restoring random number seeds from file kamil.rndm"<<G4endl;
G4RunManager::GetRunManager()->RestoreRandomNumberStatus("kamil.rndm");
boolPrintInfoAboutGeneratedParticles = true;
}
}
if (setRandomNrSeedFromFile) {
// G4cout<<"RandomNrInitialisers.size()="<<RandomNrInitialisers->size()<<G4endl;
if (eventID < int(pointerToSeedVector->size())) {
G4cout <<"musrEventAction.cc: seed will be set to="<< pointerToSeedVector->at(eventID)<<G4endl;
CLHEP::HepRandom::setTheSeed(pointerToSeedVector->at(eventID));
CLHEP::RandGauss::setFlag(false);
boolPrintInfoAboutGeneratedParticles = true;
}
}
else if (setRandomNrSeedAccordingEventNr) {
// long seeds[2];
// seeds[0] = (long) 234567890+thisEventNr*117;
// seeds[1] = (long) 333222111+thisEventNr*173;
//
// // seeds[1] = (long) (evt->GetEventID());
// // seeds[0] = (long) 123456789; // This leads to a gap in the decay time histogram fro N=100000 events
// // seeds[1] = (long) 333222111+thisEventNr; // ----------------------------||------------------------------------
// thisEventNr++;
// CLHEP::HepRandom::setTheSeeds(seeds);
// // G4cout << "seed1: " << seeds[0] << "; seed2: " << seeds[1] << G4endl;
//
// G4cout <<" thisEventNr="<<thisEventNr;
CLHEP::HepRandom::setTheSeed(eventID);
// G4cout <<" getTheSeed="<<CLHEP::HepRandom::getTheSeed()<< G4endl;
CLHEP::RandGauss::setFlag(false);
}
}
void musrPrimaryGeneratorAction::SetKEnergy(G4double val) {
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4double mu_mass = particleTable->FindParticle("mu+")->GetPDGMass();
p0=std::sqrt(val*val + 2*mu_mass*val);
// G4cout<<"musrPrimaryGeneratorAction::SetKEnergy: Muon kinetic energy of "
// <<val<<" MeV requested ==> initial muon momentum set to "<<p0<<" MeV/c"<<G4endl;
}

177
src/musrPrimaryGeneratorMessenger.cc Normal file
View File

@@ -0,0 +1,177 @@
#include "musrPrimaryGeneratorMessenger.hh"
#include "musrPrimaryGeneratorAction.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
#include "G4UIcmdWithADouble.hh"
#include "G4UIcmdWithAnInteger.hh"
#include "G4UIcmdWith3Vector.hh"
#include "G4UIcmdWith3VectorAndUnit.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGeneratorAction* musrGun)
:musrAction(musrGun)
{
setvertexCmd = new G4UIcmdWith3VectorAndUnit("/gun/vertex",this);
setvertexCmd->SetGuidance(" Set x0, y0, z0 of the generated muons (with unit)");
setvertexCmd->SetParameterName("mes_x0","mes_y0","mes_z0",true,true);
setvertexCmd->SetDefaultUnit("mm");
setvertexSigmaCmd = new G4UIcmdWith3VectorAndUnit("/gun/vertexsigma",this);
setvertexSigmaCmd->SetGuidance(" Set xSigma, ySigma, ySigma of the generated muons (with unit)");
setvertexSigmaCmd->SetParameterName("mes_xSigma","mes_ySigma","mes_zSigma",true,true);
setvertexSigmaCmd->SetDefaultUnit("mm");
setvertexBoundaryCmd = new G4UIcmdWith3VectorAndUnit("/gun/vertexboundary",this);
setvertexBoundaryCmd->SetGuidance(" Set maximum allowed radius, zmin, zmax of the generated vertex (with unit)");
setvertexBoundaryCmd->SetParameterName("mes_rMaxAllowed","mes_zMinAllowed","mes_zMaxAllowed",true,true);
setvertexBoundaryCmd->SetDefaultUnit("mm");
setKEnergyCmd = new G4UIcmdWithADoubleAndUnit("/gun/kenergy",this);
setKEnergyCmd->SetGuidance(" Set kinetic energy of the generated muons (with unit)");
setKEnergyCmd->SetParameterName("mes_E0",true);
setKEnergyCmd->SetDefaultUnit("MeV");
setMomentumCmd = new G4UIcmdWithADoubleAndUnit("/gun/momentum",this);
setMomentumCmd->SetGuidance(" Set mean momentum of the generated muons (with unit)");
setMomentumCmd->SetParameterName("mes_p0",true);
setMomentumCmd->SetDefaultUnit("MeV");
setMomentumSmearingCmd = new G4UIcmdWithADoubleAndUnit("/gun/momentumsmearing",this);
setMomentumSmearingCmd->SetGuidance(" Set sigma of the momentum of the generated muons (with unit)");
setMomentumSmearingCmd->SetParameterName("mes_pSigma",true);
setMomentumSmearingCmd->SetDefaultUnit("MeV");
setMomentumBoundaryCmd = new G4UIcmdWith3VectorAndUnit("/gun/momentumboundary",this);
setMomentumBoundaryCmd->SetGuidance(" Set minimum and maximum momentum allowed (with unit, z component ignored)");
setMomentumBoundaryCmd->SetParameterName("mes_pMinAllowed","mes_pMaxAllowed","mes_dummy",true,true);
setMomentumBoundaryCmd->SetDefaultUnit("MeV");
setTiltAngleCmd = new G4UIcmdWith3VectorAndUnit("/gun/tilt",this);
setTiltAngleCmd->SetGuidance(" Set tilt angle of the generated muons (with unit, z component ignored)");
setTiltAngleCmd->SetParameterName("mes_xangle","mes_yangle","dummy",true,true);
setTiltAngleCmd->SetDefaultUnit("deg");
setSigmaTiltAngleCmd = new G4UIcmdWith3VectorAndUnit("/gun/tiltsigma",this);
setSigmaTiltAngleCmd->SetGuidance(" Set sigma of the tilt angle (with unit, z component ignored)");
setSigmaTiltAngleCmd->SetParameterName("mes_xangleSigma","mes_yangleSigma","dummy",true,true);
setSigmaTiltAngleCmd->SetDefaultUnit("deg");
setPitchCmd = new G4UIcmdWithADoubleAndUnit("/gun/pitch",this);
setPitchCmd->SetGuidance(" Set pitch angle of the generated muons (with unit)");
setPitchCmd->SetParameterName("mes_pitch",true);
setPitchCmd->SetDefaultUnit("deg");
// setMuonPolarizCmd = new G4UIcmdWithAnInteger("/gun/muonpolarization",this);
// setMuonPolarizCmd->SetGuidance(" Set initial mu polariz: 0=transverse, 1=longitudinal ");
// setMuonPolarizCmd->SetParameterName("IniPol",true);
// setMuonPolarizCmd->SetDefaultValue(0) ;
setMuonPolarizCmd = new G4UIcmdWith3Vector("/gun/muonPolarizVector",this);
setMuonPolarizCmd->SetGuidance("Set initial mu polarisation as a vector (without unit)");
setMuonPolarizCmd->SetGuidance(" The vector does not have to be normalised to 1");
setMuonPolarizCmd->SetParameterName("mes_polarisX","mes_polarisY","mes_polarisZ",true,true);
setMuonPolarizFractionCmd = new G4UIcmdWithADouble("/gun/muonPolarizFraction",this);
setMuonPolarizFractionCmd->SetGuidance(" Set the fraction of the muon polarisation (in the range of -1 to 1),");
setMuonPolarizFractionCmd->SetGuidance(" where fraction = (N_up_spin - N_down_spin) / (N_up_spin + N_down_spin)");
setMuonPolarizFractionCmd->SetParameterName("mes_polarisFraction",true);
setMuonDecayTimeCmd = new G4UIcmdWith3VectorAndUnit("/gun/decaytimelimits",this);
setMuonDecayTimeCmd->SetGuidance(" Set minimum and maximum decay time and the muon mean life ");
setMuonDecayTimeCmd->SetParameterName("decayMin","decayMax","decayTime",true,true);
setMuonDecayTimeCmd->SetDefaultUnit("ns");
setTurtleCmd = new G4UIcmdWithAString("/gun/turtlefilename",this);
setTurtleCmd->SetGuidance("Set the filename of the TURTLE input file.");
setTurtleCmd->SetGuidance("If this varialble is set, TURTLE input will be used for initial muons");
setTurtleCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
setTurtleZ0Cmd = new G4UIcmdWithADoubleAndUnit("/gun/turtleZ0position",this);
setTurtleZ0Cmd->SetGuidance("Set the z0, with which the TURTLE input file has been generated.");
setTurtleZ0Cmd->AvailableForStates(G4State_PreInit,G4State_Idle);
setTurtleZ0Cmd->SetParameterName("mes_z0Turtle",true);
setTurtleZ0Cmd->SetDefaultUnit("mm");
setTurtleMomentumBite = new G4UIcmdWith3Vector("/gun/turtleMomentumBite",this);
setTurtleMomentumBite->SetGuidance(" Modify smearing of the turtle momentum bite. The first value is the mean momentum in MeV/c,");
setTurtleMomentumBite->SetGuidance(" the second value is the smearing factor in per cent, by which the momentum bite,");
setTurtleMomentumBite->SetGuidance(" will be increased/decreased around the mean momemtum. 100 per cent correspond to no");
setTurtleMomentumBite->SetGuidance(" change, 0 per cent will create monoenergetic beam, 200 per cent will create a two times");
setTurtleMomentumBite->SetGuidance(" broader beam. The third parameter is dummy.");
setTurtleMomentumBite->SetParameterName("mes_turtleMomentumP0","mes_turtleSmearingFactor","mes_turtleSmearingDummy",true,true);
setTurtleEventNrCmd = new G4UIcmdWithAnInteger("/gun/turtleFirstEventNr",this);
setTurtleEventNrCmd->SetGuidance("Set the line number that should be taken as the first event from the turtle input file.");
setTurtleEventNrCmd->SetParameterName("mes_turtleFirstEvent",true);
setTurtleEventNrCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
setTurtleEventNrCmd->SetDefaultValue(0) ;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
{
delete setvertexCmd;
delete setvertexSigmaCmd;
delete setvertexBoundaryCmd;
delete setKEnergyCmd;
delete setMomentumCmd;
delete setMomentumSmearingCmd;
delete setMomentumBoundaryCmd;
delete setTiltAngleCmd;
delete setSigmaTiltAngleCmd;
delete setPitchCmd;
delete setMuonPolarizCmd;
delete setMuonPolarizFractionCmd;
delete setMuonDecayTimeCmd;
delete setTurtleCmd;
delete setTurtleZ0Cmd;
delete setTurtleMomentumBite;
delete setTurtleEventNrCmd;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
void musrPrimaryGeneratorMessenger::SetNewValue(G4UIcommand * command,G4String newValue)
{
if( command == setvertexCmd)
{ musrAction->Setvertex(setvertexCmd->GetNew3VectorValue(newValue));}
if( command == setvertexSigmaCmd)
{ musrAction->SetvertexSigma(setvertexSigmaCmd->GetNew3VectorValue(newValue));}
if( command == setvertexBoundaryCmd)
{ musrAction->SetvertexBoundary(setvertexBoundaryCmd->GetNew3VectorValue(newValue));}
if( command == setKEnergyCmd)
{ musrAction->SetKEnergy(setKEnergyCmd->GetNewDoubleValue(newValue));}
if( command == setMomentumCmd)
{ musrAction->SetMomentum(setMomentumCmd->GetNewDoubleValue(newValue));}
if( command == setMomentumSmearingCmd)
{ musrAction->SetMomentumSmearing(setMomentumSmearingCmd->GetNewDoubleValue(newValue));}
if( command == setMomentumBoundaryCmd)
{ musrAction->SetMomentumBoundary(setMomentumBoundaryCmd->GetNew3VectorValue(newValue));}
if( command == setTiltAngleCmd)
{ musrAction->SetTilt(setTiltAngleCmd->GetNew3VectorValue(newValue));}
if( command == setSigmaTiltAngleCmd)
{ musrAction->SetSigmaTilt(setSigmaTiltAngleCmd->GetNew3VectorValue(newValue));}
if( command == setPitchCmd)
{ musrAction->SetPitch(setPitchCmd->GetNewDoubleValue(newValue));}
if( command == setMuonPolarizCmd)
{ musrAction->SetInitialMuonPolariz(setMuonPolarizCmd->GetNew3VectorValue(newValue));}
if( command == setMuonPolarizFractionCmd)
{ musrAction->SetInitialPolarizFraction(setMuonPolarizFractionCmd->GetNewDoubleValue(newValue));}
if( command == setMuonDecayTimeCmd)
{ musrAction->SetMuonDecayTimeLimits(setMuonDecayTimeCmd->GetNew3VectorValue(newValue)); }
if( command == setTurtleCmd)
{ musrAction->SetTurtleInput(newValue); }
if( command == setTurtleZ0Cmd)
{ musrAction->SetTurtleZ0(setTurtleZ0Cmd->GetNewDoubleValue(newValue)); }
if( command == setTurtleMomentumBite)
{ musrAction-> SetTurtleMomentumBite(setTurtleMomentumBite->GetNew3VectorValue(newValue)); }
if( command == setTurtleEventNrCmd)
{ musrAction->SetTurtleInputFileToEventNo(setTurtleEventNrCmd->GetNewIntValue(newValue));}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

87
src/musrQuadrupole.cc Normal file
View File

@@ -0,0 +1,87 @@
#include "musrQuadrupole.hh"
#include "musrParameters.hh"
musrQuadrupole::musrQuadrupole(G4double halflengthVal, G4double fieldRadiusVal, G4double gradientVal, G4double fringeFactorVal, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume)
{
G4cout << "\n-----------------------------------------------------------"
<< "\n Quadrupole field "
<< "\n-----------------------------------------------------------"
<< G4endl;
gradient = gradientVal; //*(tesla/m);
fieldRadius = fieldRadiusVal; //*mm;
halflength = halflengthVal; //*mm;
enge = BLEngeFunction(ENGE_QUAD);
G4double fringeFactor = fringeFactorVal; // the default should be 1.
G4bool fringe = true; if (fringeFactor==0.) {fringe=false;}
fringeDepth = fringeFactor * fieldRadius * 2.0;
if (!fringe) {
enge.set(0,0,0,0,0,0);
fringeMaxZ = halflength;
}
else {
for(int i=0; i<1000; ++i) {
fringeMaxZ = i*fieldRadius/10.0 + halflength;
if(enge((fringeMaxZ-halflength)/fringeDepth) < FRINGE_ACCURACY) break;
}
}
G4cout << " Field gradient set to "<< gradient/(tesla/m) << " T/m"<< G4endl;
G4cout << " Field radius set to " << fieldRadius/mm << " mm"<< G4endl;
G4cout << " Field halflength set to " << halflength/mm << " mm"<<G4endl;
G4cout << " Field fringeDepth set to " << fringeDepth/mm << " mm"<<G4endl;
G4cout << " Field fringeMaxZ set to " << fringeMaxZ/mm << " mm"<<G4endl;
G4cout << "\n-----------------------------------------------------------" << G4endl;
}
void musrQuadrupole::addFieldValue(const G4double point[4], G4double *field ) const {
G4ThreeVector global(point[0],point[1],point[2]);
G4ThreeVector local;
local = global2local.TransformPoint(global);
G4double r = sqrt(local[0]*local[0]+local[1]*local[1]);
if (r > fieldRadius || fabs(local[2]) > fringeMaxZ) { return;}
// apply enge() to the scalar potential phi=-G0*x*y*enge(z);
// B is minus its gradient. Handle both edges properly.
G4double G0 = gradient;
G4double f,fp;
if (fringeDepth!=0) {
double fringeZ = (fabs(local[2])-halflength)/fringeDepth;
f = enge(fringeZ);
fp = enge.prime(fringeZ)/fringeDepth;
}
else {
f = ( fabs(local[2]) > halflength) ? 0:1;
fp = 0;
}
G4ThreeVector B(G0*f*local[1],G0*f*local[0],G0*fp*local[0]*local[1]);
if(local[2] < 0.0) B[2] = -B[2];
G4ThreeVector finalField(B[0],B[1],B[2]);
finalField = global2local.Inverse().TransformAxis(finalField);
field[0] += finalField.x();
field[1] += finalField.y();
field[2] += finalField.z();
// G4cout<<"musrQuadrupole.cc: field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;
}
G4double musrQuadrupole::GetNominalFieldValue() {
return gradient;
}
void musrQuadrupole::SetNominalFieldValue(G4double newFieldValue) {
// // Rescale the magnetic field for a new value of the magnetic field
gradient=newFieldValue;
G4cout<<"musrQuadrupole.cc: gradient changed to="<< gradient/(tesla/m)<<" T/m"<<G4endl;
}

446
src/musrRootOutput.cc Normal file
View File

@@ -0,0 +1,446 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrRootOutput.hh"
#include "G4RunManager.hh"
#include "G4Run.hh"
#include "musrErrorMessage.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrRootOutput::musrRootOutput() {
pointerToRoot=this;
boolIsAnySpecialSaveVolumeDefined=false;
nFieldNomVal=0;
ProcessIDMapping["DecayWithSpin"]=1;
ProcessIDMapping["eIoni"]=2;
ProcessIDMapping["eBrem"]=3;
ProcessIDMapping["annihil"]=4;
ProcessIDMapping["LowEnCompton"]=5;
ProcessIDMapping["LowEnConversion"]=6;
ProcessIDMapping["LowEnBrem"]=7;
ProcessIDMapping["LowEnergyIoni"]=8;
ProcessIDMapping["LowEnPhotoElec"]=9;
ProcessIDMapping["RadioactiveDecay"]=10;
ProcessIDMapping["muIoni"]=11;
ProcessIDMapping["MuFormation"]=12;
ProcessIDMapping["Decay"]=13;
ProcessIDMapping["initialParticle"]=100;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrRootOutput::~musrRootOutput() {}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrRootOutput* musrRootOutput::pointerToRoot=0;
musrRootOutput* musrRootOutput::GetRootInstance() {
return pointerToRoot;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool musrRootOutput::store_runID = true;
G4bool musrRootOutput::store_eventID = true;
G4bool musrRootOutput::store_weight = true;
G4bool musrRootOutput::store_BFieldAtDecay = true;
G4bool musrRootOutput::store_muIniPosX = true;
G4bool musrRootOutput::store_muIniPosY = true;
G4bool musrRootOutput::store_muIniPosZ = true;
G4bool musrRootOutput::store_muIniMomX = true;
G4bool musrRootOutput::store_muIniMomY = true;
G4bool musrRootOutput::store_muIniMomZ = true;
G4bool musrRootOutput::store_muIniPolX = true;
G4bool musrRootOutput::store_muIniPolY = true;
G4bool musrRootOutput::store_muIniPolZ = true;
G4bool musrRootOutput::store_muDecayDetID= true;
G4bool musrRootOutput::store_muDecayPosX = true;
G4bool musrRootOutput::store_muDecayPosY = true;
G4bool musrRootOutput::store_muDecayPosZ = true;
G4bool musrRootOutput::store_muDecayTime = true;
G4bool musrRootOutput::store_muDecayPolX = true;
G4bool musrRootOutput::store_muDecayPolY = true;
G4bool musrRootOutput::store_muDecayPolZ = true;
G4bool musrRootOutput::store_muTargetTime = false;
G4bool musrRootOutput::store_muTargetPolX = false;
G4bool musrRootOutput::store_muTargetPolY = false;
G4bool musrRootOutput::store_muTargetPolZ = false;
G4bool musrRootOutput::store_muM0Time = false;
G4bool musrRootOutput::store_muM0PolX = false;
G4bool musrRootOutput::store_muM0PolY = false;
G4bool musrRootOutput::store_muM0PolZ = false;
G4bool musrRootOutput::store_muM1Time = false;
G4bool musrRootOutput::store_muM1PolX = false;
G4bool musrRootOutput::store_muM1PolY = false;
G4bool musrRootOutput::store_muM1PolZ = false;
G4bool musrRootOutput::store_muM2Time = false;
G4bool musrRootOutput::store_muM2PolX = false;
G4bool musrRootOutput::store_muM2PolY = false;
G4bool musrRootOutput::store_muM2PolZ = false;
G4bool musrRootOutput::store_posIniMomX = true;
G4bool musrRootOutput::store_posIniMomY = true;
G4bool musrRootOutput::store_posIniMomZ = true;
G4bool musrRootOutput::store_det_ID = true;
G4bool musrRootOutput::store_det_edep = true;
G4bool musrRootOutput::store_det_edep_el = true;
G4bool musrRootOutput::store_det_edep_pos = true;
G4bool musrRootOutput::store_det_edep_gam = true;
G4bool musrRootOutput::store_det_edep_mup = true;
G4bool musrRootOutput::store_det_nsteps = true;
G4bool musrRootOutput::store_det_length = true;
G4bool musrRootOutput::store_det_start = true;
G4bool musrRootOutput::store_det_end = true;
G4bool musrRootOutput::store_det_x = true;
G4bool musrRootOutput::store_det_y = true;
G4bool musrRootOutput::store_det_z = true;
G4bool musrRootOutput::store_det_kine = true;
G4bool musrRootOutput::store_det_VrtxKine = true;
G4bool musrRootOutput::store_det_VrtxX = true;
G4bool musrRootOutput::store_det_VrtxY = true;
G4bool musrRootOutput::store_det_VrtxZ = true;
G4bool musrRootOutput::store_det_VrtxVolID = true;
G4bool musrRootOutput::store_det_VrtxProcID = true;
G4bool musrRootOutput::store_det_VrtxTrackID = true;
G4bool musrRootOutput::store_det_VrtxParticleID = true;
G4bool musrRootOutput::store_det_VvvKine = true;
G4bool musrRootOutput::store_det_VvvX = true;
G4bool musrRootOutput::store_det_VvvY = true;
G4bool musrRootOutput::store_det_VvvZ = true;
G4bool musrRootOutput::store_det_VvvVolID = true;
G4bool musrRootOutput::store_det_VvvProcID = true;
G4bool musrRootOutput::store_det_VvvTrackID = true;
G4bool musrRootOutput::store_det_VvvParticleID = true;
G4bool musrRootOutput::store_fieldNomVal = true;
G4bool musrRootOutput::store_fieldIntegralBx = false;
G4bool musrRootOutput::store_fieldIntegralBy = false;
G4bool musrRootOutput::store_fieldIntegralBz = false;
G4bool musrRootOutput::store_fieldIntegralBz1 = false;
G4bool musrRootOutput::store_fieldIntegralBz2 = false;
G4bool musrRootOutput::store_fieldIntegralBz3 = false;
G4int musrRootOutput::oldEventNumberInG4EqEMFieldWithSpinFunction=-1;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::BeginOfRunAction() {
G4cout << "musrRootOutput::BeginOfRunAction() Defining the Root tree and branches:"<<G4endl;
G4int tmpRunNr=(G4RunManager::GetRunManager())->GetCurrentRun()->GetRunID();
char RootOutputFileName[200];
sprintf(RootOutputFileName, "data/musr_%i.root", tmpRunNr);
rootFile=new TFile(RootOutputFileName,"recreate");
rootTree=new TTree("t1","a simple Tree with simple variables");
if (store_runID) {rootTree->Branch("runID",&runID_t,"runID/I");}
if (store_eventID) {rootTree->Branch("eventID",&eventID_t,"eventID/I");}
if (store_weight) {rootTree->Branch("weight",&weight_t,"weight/D");}
if (store_BFieldAtDecay) {rootTree->Branch("BFieldAtDecay",&B_t,"Bx/D:By:Bz:B3:B4:B5");}
if (store_muIniPosX) {rootTree->Branch("muIniPosX",&muIniPosX_t,"muIniPosX/D");}
if (store_muIniPosY) {rootTree->Branch("muIniPosY",&muIniPosY_t,"muIniPosY/D");}
if (store_muIniPosZ) {rootTree->Branch("muIniPosZ",&muIniPosZ_t,"muIniPosZ/D");}
if (store_muIniMomX) {rootTree->Branch("muIniMomX",&muIniMomX_t,"muIniMomX/D");}
if (store_muIniMomY) {rootTree->Branch("muIniMomY",&muIniMomY_t,"muIniMomY/D");}
if (store_muIniMomZ) {rootTree->Branch("muIniMomZ",&muIniMomZ_t,"muIniMomZ/D");}
if (store_muIniPolX) {rootTree->Branch("muIniPolX",&muIniPolX_t,"muIniPolX/D");}
if (store_muIniPolY) {rootTree->Branch("muIniPolY",&muIniPolY_t,"muIniPolY/D");}
if (store_muIniPolZ) {rootTree->Branch("muIniPolZ",&muIniPolZ_t,"muIniPolZ/D");}
if (store_muDecayDetID) {rootTree->Branch("muDecayDetID",&muDecayDetID_t,"muDecayDetID/I");}
if (store_muDecayPosX) {rootTree->Branch("muDecayPosX",&muDecayPosX_t,"muDecayPosX/D");}
if (store_muDecayPosY) {rootTree->Branch("muDecayPosY",&muDecayPosY_t,"muDecayPosY/D");}
if (store_muDecayPosZ) {rootTree->Branch("muDecayPosZ",&muDecayPosZ_t,"muDecayPosZ/D");}
if (store_muDecayTime) {rootTree->Branch("muDecayTime",&muDecayTime_t,"muDecayTime/D");}
if (store_muDecayPolX) {rootTree->Branch("muDecayPolX",&muDecayPolX_t,"muDecayPolX/D");}
if (store_muDecayPolY) {rootTree->Branch("muDecayPolY",&muDecayPolY_t,"muDecayPolY/D");}
if (store_muDecayPolZ) {rootTree->Branch("muDecayPolZ",&muDecayPolZ_t,"muDecayPolZ/D");}
if (store_muTargetTime) {rootTree->Branch("muTargetTime",&muTargetTime_t,"muTargetTime/D");}
if (store_muTargetPolX) {rootTree->Branch("muTargetPolX",&muTargetPolX_t,"muTargetPolX/D");}
if (store_muTargetPolY) {rootTree->Branch("muTargetPolY",&muTargetPolY_t,"muTargetPolY/D");}
if (store_muTargetPolZ) {rootTree->Branch("muTargetPolZ",&muTargetPolZ_t,"muTargetPolZ/D");}
if (store_muM0Time) {rootTree->Branch("muM0Time",&muM0Time_t,"muM0Time/D");}
if (store_muM0PolX) {rootTree->Branch("muM0PolX",&muM0PolX_t,"muM0PolX/D");}
if (store_muM0PolY) {rootTree->Branch("muM0PolY",&muM0PolY_t,"muM0PolY/D");}
if (store_muM0PolZ) {rootTree->Branch("muM0PolZ",&muM0PolZ_t,"muM0PolZ/D");}
if (store_muM1Time) {rootTree->Branch("muM1Time",&muM1Time_t,"muM1Time/D");}
if (store_muM1PolX) {rootTree->Branch("muM1PolX",&muM1PolX_t,"muM1PolX/D");}
if (store_muM1PolY) {rootTree->Branch("muM1PolY",&muM1PolY_t,"muM1PolY/D");}
if (store_muM1PolZ) {rootTree->Branch("muM1PolZ",&muM1PolZ_t,"muM1PolZ/D");}
if (store_muM2Time) {rootTree->Branch("muM2Time",&muM2Time_t,"muM2Time/D");}
if (store_muM2PolX) {rootTree->Branch("muM2PolX",&muM2PolX_t,"muM2PolX/D");}
if (store_muM2PolY) {rootTree->Branch("muM2PolY",&muM2PolY_t,"muM2PolY/D");}
if (store_muM2PolZ) {rootTree->Branch("muM2PolZ",&muM2PolZ_t,"muM2PolZ/D");}
if (store_posIniMomX) {rootTree->Branch("posIniMomX",&posIniMomx_t,"posIniMomX/D");}
if (store_posIniMomY) {rootTree->Branch("posIniMomY",&posIniMomy_t,"posIniMomY/D");}
if (store_posIniMomZ) {rootTree->Branch("posIniMomZ",&posIniMomz_t,"posIniMomZ/D");}
// if (store_globalTime) {rootTree->Branch("globalTime",&globalTime_t,"globalTime/D");}
// if (store_fieldValue) {rootTree->Branch("fieldValue",&fieldValue_t,"fieldValue/D");}
if (store_fieldNomVal) {
rootTree->Branch("nFieldNomVal",&nFieldNomVal,"nFieldNomVal/I");
rootTree->Branch("fieldNomVal",&fieldNomVal,"fieldNomVal[nFieldNomVal]/D");
}
if (store_fieldIntegralBx) {rootTree->Branch("BxIntegral",&BxIntegral_t,"BxIntegral/D");}
if (store_fieldIntegralBy) {rootTree->Branch("ByIntegral",&ByIntegral_t,"ByIntegral/D");}
if (store_fieldIntegralBz) {rootTree->Branch("BzIntegral",&BzIntegral_t,"BzIntegral/D");}
if (store_fieldIntegralBz1) {rootTree->Branch("BzIntegral1",&BzIntegral1_t,"BzIntegral1/D");}
if (store_fieldIntegralBz2) {rootTree->Branch("BzIntegral2",&BzIntegral2_t,"BzIntegral2/D");}
if (store_fieldIntegralBz3) {rootTree->Branch("BzIntegral3",&BzIntegral3_t,"BzIntegral3/D");}
rootTree->Branch("det_n",&det_n,"det_n/I");
if (store_det_ID) {rootTree->Branch("det_ID",&det_ID,"det_ID[det_n]/I");}
if (store_det_edep) {rootTree->Branch("det_edep",&det_edep,"det_edep[det_n]/D");}
if (store_det_edep_el) {rootTree->Branch("det_edep_el",&det_edep_el,"det_edep_el[det_n]/D");}
if (store_det_edep_pos) {rootTree->Branch("det_edep_pos",&det_edep_pos,"det_edep_pos[det_n]/D");}
if (store_det_edep_gam) {rootTree->Branch("det_edep_gam",&det_edep_gam,"det_edep_gam[det_n]/D");}
if (store_det_edep_mup) {rootTree->Branch("det_edep_mup",&det_edep_mup,"det_edep_mup[det_n]/D");}
if (store_det_nsteps) {rootTree->Branch("det_nsteps",&det_nsteps,"det_nsteps[det_n]/I");}
if (store_det_length) {rootTree->Branch("det_length",&det_length,"det_length[det_n]/D");}
if (store_det_start) {rootTree->Branch("det_time_start",&det_time_start,"det_time_start[det_n]/D");}
if (store_det_end) {rootTree->Branch("det_time_end",&det_time_end,"det_time_end[det_n]/D");}
if (store_det_x) {rootTree->Branch("det_x",&det_x,"det_x[det_n]/D");}
if (store_det_y) {rootTree->Branch("det_y",&det_y,"det_y[det_n]/D");}
if (store_det_z) {rootTree->Branch("det_z",&det_z,"det_z[det_n]/D");}
if (store_det_kine) {rootTree->Branch("det_kine",&det_kine,"det_kine[det_n]/D");}
if (store_det_VrtxKine) {rootTree->Branch("det_VrtxKine",&det_VrtxKine,"det_VrtxKine[det_n]/D");}
if (store_det_VrtxX) {rootTree->Branch("det_VrtxX",&det_VrtxX,"det_VrtxX[det_n]/D");}
if (store_det_VrtxY) {rootTree->Branch("det_VrtxY",&det_VrtxY,"det_VrtxY[det_n]/D");}
if (store_det_VrtxZ) {rootTree->Branch("det_VrtxZ",&det_VrtxZ,"det_VrtxZ[det_n]/D");}
if (store_det_VrtxVolID){rootTree->Branch("det_VrtxVolID",&det_VrtxVolID,"det_VrtxVolID[det_n]/I");}
if (store_det_VrtxProcID){rootTree->Branch("det_VrtxProcID",&det_VrtxProcID,"det_VrtxProcID[det_n]/I");}
if (store_det_VrtxTrackID){rootTree->Branch("det_VrtxTrackID",&det_VrtxTrackID,"det_VrtxTrackID[det_n]/I");}
if (store_det_VrtxParticleID){rootTree->Branch("det_VrtxParticleID",&det_VrtxParticleID,"det_VrtxParticleID[det_n]/I");}
if (store_det_VvvKine) {rootTree->Branch("det_VvvKine",&det_VvvKine,"det_VvvKine[det_n]/D");}
if (store_det_VvvX) {rootTree->Branch("det_VvvX",&det_VvvX,"det_VvvX[det_n]/D");}
if (store_det_VvvY) {rootTree->Branch("det_VvvY",&det_VvvY,"det_VvvY[det_n]/D");}
if (store_det_VvvZ) {rootTree->Branch("det_VvvZ",&det_VvvZ,"det_VvvZ[det_n]/D");}
if (store_det_VvvVolID){rootTree->Branch("det_VvvVolID",&det_VvvVolID,"det_VvvVolID[det_n]/I");}
if (store_det_VvvProcID){rootTree->Branch("det_VvvProcID",&det_VvvProcID,"det_VvvProcID[det_n]/I");}
if (store_det_VvvTrackID){rootTree->Branch("det_VvvTrackID",&det_VvvTrackID,"det_VvvTrackID[det_n]/I");}
if (store_det_VvvParticleID){rootTree->Branch("det_VvvParticleID",&det_VvvParticleID,"det_VvvParticleID[det_n]/I");}
if (boolIsAnySpecialSaveVolumeDefined) {
rootTree->Branch("save_n",&save_n,"save_n/I");
rootTree->Branch("save_detID",&save_detID,"save_detID[save_n]/I");
rootTree->Branch("save_particleID",&save_particleID,"save_particleID[save_n]/I");
rootTree->Branch("save_ke",&save_ke,"save_ke[save_n]/D");
rootTree->Branch("save_x",&save_x,"save_x[save_n]/D");
rootTree->Branch("save_y",&save_y,"save_y[save_n]/D");
rootTree->Branch("save_z",&save_z,"save_z[save_n]/D");
rootTree->Branch("save_px",&save_px,"save_px[save_n]/D");
rootTree->Branch("save_py",&save_py,"save_py[save_n]/D");
rootTree->Branch("save_pz",&save_pz,"save_pz[save_n]/D");
}
// htest1 = new TH1F("htest1","The debugging histogram 1",50,-4.,4.);
// htest2 = new TH1F("htest2","The debugging histogram 2",50,0.,3.142);
htest1 = new TH2F("htest1","x, y",50,-200.,200.,50,-200.,200.);
htest2 = new TH2F("htest2","R, z",50,0.,250.,50,-150.,150.);
htest3 = new TH1F("htest3","Energy in MeV",55,0.,55.);
htest4 = new TH1F("htest4","Radioactive electron kinetic energy",250,0.,2.5);
htest5 = new TH1F("htest5","The debugging histogram 5",50,-4.,4.);
htest6 = new TH1F("htest6","The debugging histogram 6",50,0.,3.142);
htest7 = new TH1F("htest7","The debugging histogram 7",50,-4.,4.);
htest8 = new TH1F("htest8","The debugging histogram 8",50,0.,3.142);
G4cout << "musrRootOutput::BeginOfRunAction() The Root tree and branches were defined."<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::EndOfRunAction() {
G4cout << "musrRootOutput::EndOfRunAction() - Writing out the Root tree:"<<G4endl;
rootTree->Write();
htest1->Write();
htest2->Write();
htest3->Write();
htest4->Write();
htest5->Write();
htest6->Write();
htest7->Write();
htest8->Write();
// Variables exported from Geant simulation to the Root output
// static const Int_t nGeantParamD=10;
TVectorD TVector_GeantParametersD(maxNGeantParameters);
for (Int_t i=0; i<maxNGeantParameters; i++) {
TVector_GeantParametersD[i]=GeantParametersD[i];
}
TVector_GeantParametersD.Write("geantParametersD");
rootFile->Close();
G4cout<<"musrRootOutput::EndOfRunAction() - Root tree written out."<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::FillEvent() {
htest5->Fill(atan2(posIniMomy_t,posIniMomx_t));
htest6->Fill(atan2(sqrt(posIniMomx_t*posIniMomx_t+posIniMomy_t*posIniMomy_t),posIniMomz_t));
if (weight_t>0.) {
rootTree->Fill();
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::ClearAllRootVariables() {
runID_t=-1000;
eventID_t=-1000;
weight_t=1.;
B_t[0]=-1000.;B_t[1]=-1000.;B_t[2]=-1000.;B_t[3]=-1000.;B_t[4]=-1000.;B_t[5]=-1000.;
muIniPosX_t=-1000; muIniPosY_t=-1000; muIniPosZ_t=-1000;
muIniMomX_t=-1000; muIniMomY_t=-1000; muIniMomZ_t=-1000;
muIniPolX_t=-1000; muIniPolY_t=-1000; muIniPolZ_t=-1000;
muDecayDetID_t=-1000;
muDecayPolX_t=-1000; muDecayPolY_t=-1000; muDecayPolZ_t=-1000;
muTargetTime_t=-1000; muTargetPolX_t=-1000; muTargetPolY_t=-1000; muTargetPolZ_t=-1000;
muM0Time_t=-1000; muM0PolX_t=-1000; muM0PolY_t=-1000; muM0PolZ_t=-1000;
muM1Time_t=-1000; muM1PolX_t=-1000; muM1PolY_t=-1000; muM1PolZ_t=-1000;
muM2Time_t=-1000; muM2PolX_t=-1000; muM2PolY_t=-1000; muM2PolZ_t=-1000;
muDecayPosX_t=-1000;muDecayPosY_t=-1000;muDecayPosZ_t=-1000;
muDecayTime_t=-1000;
posIniMomx_t=-1000;posIniMomy_t=-1000;posIniMomz_t=-1000;
BxIntegral_t = -1000; ByIntegral_t = -1000; BzIntegral_t = -1000;
BzIntegral1_t = -1000; BzIntegral2_t = -1000; BzIntegral3_t = -1000;
det_n=0;
save_n=0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::SetVolumeIDMapping(std::string logivol, int volumeID) {
// This function assigns a unique number to each sensitive detector name.
// The numbers are used in the root tree, as it is easier to work with numbers
// rather than with strings.
if (SensDetectorMapping[logivol]) {
char message[200];
sprintf(message,"musrRootOutput::SetVolumeIDMapping: Sensitive volume %s already assigned",logivol.c_str());
musrErrorMessage::GetInstance()->musrError(FATAL,message,false);
}
else{
SensDetectorMapping[logivol]=volumeID;
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int musrRootOutput::ConvertVolumeToID(std::string logivol) {
G4int volumeID = SensDetectorMapping[logivol];
if (volumeID==0) {
char message[200];
sprintf(message,"musrRootOutput::ConvertVolumeToID: No ID number assigned to sensitive volume %s .",logivol.c_str());
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,true);
}
return volumeID;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int musrRootOutput::ConvertProcessToID(std::string processName) {
G4int processID = ProcessIDMapping[processName];
if (processID==0) {
char message[200];
sprintf(message,"musrRootOutput::ConvertProcessToID: No ID number assigned to the process \"%s\" .",processName.c_str());
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
return processID;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke,
G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz) {
if (save_n>=save_nMax) {
char message[200];
sprintf(message,"musrRootOutput.cc::SetSaveDetectorInfo(): more \"save\" hits then allowed: save_nMax=%i",save_nMax);
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,true);
}
else {
save_detID[save_n]=ID;
save_particleID[save_n]=particleID;
save_ke[save_n]=ke/MeV;
save_x[save_n]=x/mm;
save_y[save_n]=y/mm;
save_z[save_n]=z/mm;
save_px[save_n]=px/MeV;
save_py[save_n]=py/MeV;
save_pz[save_n]=pz/MeV;
save_n++;
}
}
void musrRootOutput::SetFieldNomVal(G4int i, G4double value) {
if (i<maxNFieldnNominalValues) {
// cks the following will probably not be correct for electric field,
// because the units are tesla. Should be modified.
fieldNomVal[i]=value/tesla;
}
else {
char message[200];
sprintf(message,
"musrRootOutput.cc::SetFieldNomVal(): more electromagnetic fields then allowed: maxNFieldnNominalValues=%i",
maxNFieldnNominalValues);
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,true);
}
}
void musrRootOutput::SetDetectorInfo (G4int nDetectors, G4int ID, G4int particleID, G4double edep,
G4double edep_el, G4double edep_pos,
G4double edep_gam, G4double edep_mup,G4int nsteps, G4double length, G4double t1,
G4double t2, G4double x, G4double y, G4double z,
G4double ek, G4double ekVertex, G4double xVertex, G4double yVertex, G4double zVertex,
G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex)
{
if ((nDetectors<0)||(nDetectors>=(det_nMax-1))) {
char message[200];
sprintf(message,"musrRootOutput.cc::SetDetectorInfo: nDetectors %i is larger than det_nMax = %i",nDetectors,det_nMax);
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,false);
return;
}
else {
det_n=nDetectors+1;
det_ID[nDetectors]=ID;
det_edep[nDetectors]=edep/MeV;
det_edep_el[nDetectors]=edep_el/MeV;
det_edep_pos[nDetectors]=edep_pos/MeV;
det_edep_gam[nDetectors]=edep_gam/MeV;
det_edep_mup[nDetectors]=edep_mup/MeV;
det_nsteps[nDetectors]=nsteps;
det_length[nDetectors]=length/mm;
det_time_start[nDetectors]=t1/microsecond;
det_time_end[nDetectors]=t2/microsecond;
det_x[nDetectors]=x/mm;
det_y[nDetectors]=y/mm;
det_z[nDetectors]=z/mm;
det_kine[nDetectors]=ek/MeV;
det_VrtxKine[nDetectors]=ekVertex/MeV;
det_VrtxX[nDetectors]=xVertex/mm;
det_VrtxY[nDetectors]=yVertex/mm;
det_VrtxZ[nDetectors]=zVertex/mm;
det_VrtxVolID[nDetectors]=idVolVertex;
det_VrtxProcID[nDetectors]=idProcVertex;
det_VrtxTrackID[nDetectors]=idTrackVertex;
det_VrtxParticleID[nDetectors]=particleID;
}
}
void musrRootOutput::SetDetectorInfoVvv (G4int nDetectors,
G4double ekVertex, G4double xVertex, G4double yVertex, G4double zVertex,
G4int idVolVertex, G4int idProcVertex, G4int idTrackVertex, G4int particleID) {
if ((nDetectors<0)||(nDetectors>=(det_nMax-1))) {
char message[200];
sprintf(message,"musrRootOutput.cc::SetDetectorInfoVvv: nDetectors %i is larger than det_nMax = %i",nDetectors,det_nMax);
musrErrorMessage::GetInstance()->musrError(SERIOUS,message,false);
return;
}
else {
det_VvvKine[nDetectors]=ekVertex/MeV;
det_VvvX[nDetectors]=xVertex/mm;
det_VvvY[nDetectors]=yVertex/mm;
det_VvvZ[nDetectors]=zVertex/mm;
det_VvvVolID[nDetectors]=idVolVertex;
det_VvvProcID[nDetectors]=idProcVertex;
det_VvvTrackID[nDetectors]=idTrackVertex;
det_VvvParticleID[nDetectors]=particleID;
}
}

100
src/musrRunAction.cc Normal file
View File

@@ -0,0 +1,100 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
// Make G4Timer appear first!
#include "G4Timer.hh"
#include "musrRunAction.hh"
#include "musrEventAction.hh"
#include "musrSteppingAction.hh"
#include "G4Run.hh"
#include "musrErrorMessage.hh"
#include "F04GlobalField.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrRunAction::musrRunAction() {
timer = new G4Timer;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrRunAction::~musrRunAction() {
delete timer;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRunAction::BeginOfRunAction(const G4Run* aRun) {
timer->Start();
G4int run_id= aRun->GetRunID();
// if (run_id%100 == 0) {
G4cout << "### Run " << run_id << G4endl;
musrRootOutput::GetRootInstance()->StoreGeantParameter(6,run_id);
musrRootOutput::GetRootInstance()->BeginOfRunAction();
// Initiate global electromagnetic field (if it exists):
if (F04GlobalField::Exists()) {
FieldList* fields = F04GlobalField::getObject()->getFields();
if (fields) {
if (fields->size()>0) {
G4int jjj=0;
G4cout<<"\n------------ The following fields were defined: ----------------"<<G4endl;
FieldList::iterator i;
for (i=fields->begin(); i!=fields->end(); ++i) {
(*i)->construct();
// G4cout<<"\t==> "<<(*i)->GetElementFieldName()<<G4endl;
// Get the nominal field value for the given field and store it in the Root output
G4double nomFieldValue = (*i)->GetNominalFieldValue();
musrRootOutput::GetRootInstance()->SetFieldNomVal(jjj,nomFieldValue);
jjj++;
}
G4cout<<"-----------------------------------------------------------------"<<G4endl;
}
}
// Print out the field values at the points user requested to be printed out:
F04GlobalField::getObject()->PrintFieldAtRequestedPoints();
}
// initialise the "boolIsVvvInfoRequested" variable in the musrSteppingAction.cc and musrScintSD.cc
G4bool boolIsVvvInfoRequested = false;
if ((musrRootOutput::store_det_VvvKine)||(musrRootOutput::store_det_VvvX)||
(musrRootOutput::store_det_VvvY)||(musrRootOutput::store_det_VvvZ)||
(musrRootOutput::store_det_VvvVolID)||(musrRootOutput::store_det_VvvProcID)||
(musrRootOutput::store_det_VvvTrackID)||(musrRootOutput::store_det_VvvParticleID)) {
boolIsVvvInfoRequested = true;
}
musrSteppingAction::GetInstance()->SetVvvInfoRequested(boolIsVvvInfoRequested);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRunAction::EndOfRunAction(const G4Run* aRun) {
musrRootOutput::GetRootInstance()->StoreGeantParameter(5,aRun->GetNumberOfEvent());
musrRootOutput::GetRootInstance()->EndOfRunAction();
if (F04GlobalField::Exists()) {
F04GlobalField* myGlobalField = F04GlobalField::getObject(); //Was causing seg. fault when called here
if (myGlobalField!=NULL) {delete myGlobalField;}
}
musrErrorMessage::GetInstance()->PrintErrorSummary();
timer->Stop();
G4cout << "musrRunAction::EndOfRunAction:"<<G4endl;
G4cout << " Number of events = " << aRun->GetNumberOfEvent()<<G4endl;
// << " " << *timer << G4endl;
G4cout << " User elapsed time = "<<timer->GetUserElapsed()/3600<<"h = "
<<timer->GetUserElapsed()/60<<"min = "<<timer->GetUserElapsed()<<"s."<<G4endl;
G4cout << " Real elapsed time = "<<timer->GetRealElapsed()/3600<<"h = "
<<timer->GetRealElapsed()/60<<"min = "<<timer->GetRealElapsed()<<"s."<<G4endl;
G4cout << " System elapsed time = "<<timer->GetSystemElapsed()/3600<<"h = "
<<timer->GetSystemElapsed()/60<<"min = "<<timer->GetSystemElapsed()<<"s."<<G4endl;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

156
src/musrScintHit.cc Normal file
View File

@@ -0,0 +1,156 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrScintHit.hh"
#include "G4UnitsTable.hh"
#include "G4VVisManager.hh"
#include "G4Circle.hh"
#include "G4Colour.hh"
#include "G4VisAttributes.hh"
#include "G4ios.hh"
#include "G4MagneticField.hh"
#include "G4FieldManager.hh"
#include "G4TransportationManager.hh"
#include "globals.hh"
#include "G4Transform3D.hh"
#include "G4ProcessManager.hh"
#include "G4Track.hh"
#include "G4ThreeVector.hh"
#include "G4RunManager.hh"
#include "G4Run.hh"
#include <fstream>
#include <iostream>
#include <iomanip>
G4Allocator<musrScintHit> musrScintHitAllocator;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintHit::musrScintHit() {
}
G4int musrScintHit::ScintMultihit=0;
G4int musrScintHit::runIDoldScint=-1;
G4int musrScintHit::eventIDoldScint=-1;
G4int musrScintHit::NIS=0;
G4int musrScintHit::ScintChamberNbold=-1;
G4int musrScintHit::verboseLevel=0;
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintHit::~musrScintHit() {
//save the Tree header. The file will be automatically closed
//when going out of the function scope
// rootTree->Write();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintHit::musrScintHit(const musrScintHit& right)
: G4VHit()
{
particleName = right.particleName;
trackID = right.trackID;
edep = right.edep;
pre_pos = right.pre_pos;
pol = right.pol;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
const musrScintHit& musrScintHit::operator=(const musrScintHit& right)
{
particleName = right.particleName;
trackID = right.trackID;
edep = right.edep;
pre_pos = right.pre_pos;
pol = right.pol;
return *this;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4int musrScintHit::operator==(const musrScintHit& right) const
{
return (this==&right) ? 1 : 0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintHit::Draw()
{
G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance();
if(pVVisManager)
{
G4Circle circle(pre_pos);
circle.SetScreenSize(0.04);
circle.SetFillStyle(G4Circle::filled);
G4Colour colour(1.,0.,0.);
G4VisAttributes attribs(colour);
circle.SetVisAttributes(attribs);
pVVisManager->Draw(circle);
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintHit::Print()
{
if (verboseLevel>2) G4cout<<"VERBOSE 3: Kamil: musrScintHit::Print()"<<G4endl;
G4RunManager* fRunManager = G4RunManager::GetRunManager();
eventID = fRunManager->GetCurrentEvent()->GetEventID();
runID = fRunManager->GetCurrentRun()->GetRunID();
G4int ScintMultihitSwitch=0;
if (runID != runIDoldScint) {
NIS=0;
ScintMultihit = 0;
eventIDoldScint = -1;
ScintChamberNbold = -1;
}
//cks if (particleName== "e+" and (eventIDoldScint != eventID or ScintChamberNbold != IBchamberNb)) {
if (particleName== "e+") {
if (eventIDoldScint == eventID) {
ScintMultihit++;
ScintMultihitSwitch=1;
}
NIS++;
G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
point[0]=0.;
point[1]=0.;
point[2]=0.;
B[2]=0.0;
if(!fMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field
mfield = fMgr->GetDetectorField();
mfield->GetFieldValue(point,B);
B[0]=B[0]/tesla;
B[1]=B[1]/tesla;
B[2]=B[2]/tesla;
}
// G4cout << " Segment: " << IBchamberNb << G4endl;
// G4cout <<"Position " << pos.x()/cm<<" "<<pos.y()/cm <<" "<< pos.z()/cm <<G4endl;
// G4cout << "Field is "<< B[2]<<G4endl;
std::ofstream posfile1;
posfile1.open ("scint.dat", std::ios::out | std::ios::app);
posfile1 << runID << " " << eventID
<< " " << logicalVolume
<< " " << ScintMultihitSwitch
<<" " << edep
<< " " << fabs(B[2])
<<" "<< pre_pos.x()/cm<<" "<<pre_pos.y()/cm <<" "<< pre_pos.z()/cm
<< " " << globalTime/s
// << " " << IBchamberNb
// << " first=" << firstStepInVolume << " last=" << lastStepInVolume
<< G4endl;
posfile1.close();
}
eventIDoldScint=eventID;
runIDoldScint = runID;
// ScintChamberNbold = IBchamberNb;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

359
src/musrScintSD.cc Normal file
View File

@@ -0,0 +1,359 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrScintSD.hh"
#include "G4HCofThisEvent.hh"
#include "G4Step.hh"
#include "G4ThreeVector.hh"
#include "G4SDManager.hh"
#include "G4ios.hh"
#include <algorithm> // needed for the sort() function
#include "G4VProcess.hh" // needed for the degugging message of the process name
#include "G4RunManager.hh"
#include "G4Run.hh"
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
#include "musrSteppingAction.hh"
#include <vector>
//bool myREMOVEfunction (int i,int j) { return (i<j); }
//bool timeOrdering (musrScintHit hit1, musrScintHit hit2) {
// return (hit1.GetGlobalTime()<hit2.GetGlobalTime());
//}
//
//bool timeOrdering2 (std::map<int,double>::iterator i1, std::map<int,double>::iterator m2) {
// return ( (*i1).first()<(*i2).second() );
//}
//
//bool timeOrdering2 (std::pair<int,double> p1, std::pair<int,double> p2) {
// return ( p1.first()<p2.second() );
//}
//
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintSD::musrScintSD(G4String name)
:G4VSensitiveDetector(name)
{
G4String HCname;
collectionName.insert(HCname="scintCollection");
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrScintSD::~musrScintSD(){ }
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::Initialize(G4HCofThisEvent* HCE) {
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::Initialize\n";
scintCollection = new musrScintHitsCollection
(SensitiveDetectorName,collectionName[0]);
static G4int HCID = -1;
if(HCID<0) {
HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::HCID was <0\n, now HCID="<<HCID<<"\n";
}
HCE->AddHitsCollection( HCID, scintCollection );
myStoreOnlyEventsWithHits = musrParameters::storeOnlyEventsWithHits;
mySignalSeparationTime = musrParameters::signalSeparationTime;
myStoreOnlyTheFirstTimeHit= musrParameters::storeOnlyTheFirstTimeHit;
myStoreOnlyEventsWithHitInDetID = musrParameters::storeOnlyEventsWithHitInDetID;
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
boolIsVvvInfoRequested = myMusrSteppingAction->IsVvvInfoRequested();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
G4bool musrScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*)
{
if (verboseLevel>1) G4cout<<"VERBOSE 2: musrScintSD::ProcessHits\n";
G4double edep = aStep->GetTotalEnergyDeposit();
if(edep==0.) {
return false;
}
G4Track* aTrack = aStep->GetTrack();
G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName();
// If requested, store only the hit that happened first (usefull for some special studies, not for a serious simulation)
if (myStoreOnlyTheFirstTimeHit) {
G4int NbHits = scintCollection->entries();
if (NbHits>0) {
aTrack->SetTrackStatus(fStopAndKill);
return false;
}
}
musrScintHit* newHit = new musrScintHit();
newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName());
G4int particleID = aTrack->GetDefinition()->GetPDGEncoding();
newHit->SetParticleID (particleID);
newHit->SetEdep (edep);
newHit->SetPrePos (aStep->GetPreStepPoint()->GetPosition());
newHit->SetPostPos (aStep->GetPostStepPoint()->GetPosition());
newHit->SetPol (aTrack->GetPolarization());
G4LogicalVolume* hitLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
newHit->SetLogVolName(hitLogicalVolume->GetName());
newHit->SetGlobTime(aTrack->GetGlobalTime());
// Warning - aStep->IsFirstStepInVolume() only available in Geant version >= 4.8.2 !
// newHit->SetFirstStepInVolumeFlag(aStep->IsFirstStepInVolume());
// newHit->SetLastStepInVolumeFlag(aStep->IsLastStepInVolume());
newHit->SetKineticEnergy(aTrack->GetKineticEnergy());
// newHit->SetKineticEnergy(aTrack->GetKineticEnergy()+edep);
G4double vertexKineticEnergy = aTrack->GetVertexKineticEnergy();
newHit->SetVertexKineticEnergy(vertexKineticEnergy);
G4ThreeVector vertexPosition = aTrack->GetVertexPosition();
newHit->SetVertexPosition(vertexPosition);
const G4LogicalVolume* vertexLogicalVolume = aTrack->GetLogicalVolumeAtVertex();
G4String vertexLogicalVolumeName = vertexLogicalVolume->GetName();
newHit->SetLogicalVolumeAtVertex(vertexLogicalVolumeName);
G4String processName;
if ((aTrack->GetCreatorProcess())!=0) { processName=aTrack->GetCreatorProcess()->GetProcessName(); }
else {processName="initialParticle";} //if no process found, the track comes from the generated particle
newHit->SetCreatorProcessName(processName);
G4int trackID = aTrack->GetTrackID();
newHit->SetTrackID (trackID);
newHit->SetStepLength (aStep->GetStepLength());
scintCollection->insert( newHit );
// newHit->Print();
newHit->Draw();
return true;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrScintSD::EndOfEvent(G4HCofThisEvent*) {
if (verboseLevel>1) {
G4cout<<"VERBOSE 2: musrScintSD::EndOfEvent"<<G4endl;
G4int NbHits = scintCollection->entries();
G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
<< " hits in the scint chambers: " << G4endl;
}
// Positron_momentum_already_stored=0;
musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
G4RunManager* fRunManager = G4RunManager::GetRunManager();
myRootOutput->SetRunID(fRunManager->GetCurrentRun()->GetRunID());
myRootOutput->SetEventID(fRunManager->GetCurrentEvent()->GetEventID());
G4int NbHits = scintCollection->entries();
if (myStoreOnlyEventsWithHits) {
if (NbHits<=0) {
return;
}
else if (myStoreOnlyEventsWithHitInDetID!=0) {
for (G4int i=0; i<NbHits; i++) {
musrScintHit* aHit = (*scintCollection)[i];
G4String aHitVolumeName = aHit->GetLogVolName();
G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName);
if (aHitVolumeID==myStoreOnlyEventsWithHitInDetID) break; // hit in the requested detector was identified
if (i==(NbHits-1)) return; // no hit identified in the requested detector
}
}
}
// Sort out hits and fill them into root
if (NbHits>0) {
const G4int det_IDmax = musrRootOutput::det_nMax;
G4double det_edep[det_IDmax];
G4int det_nsteps[det_IDmax];
G4double det_length[det_IDmax];
G4int det_ID[det_IDmax];
G4double det_edep_el[det_IDmax];
G4double det_edep_pos[det_IDmax];
G4double det_edep_gam[det_IDmax];
G4double det_edep_mup[det_IDmax];
G4double det_time_start[det_IDmax];
G4double det_time_end[det_IDmax];
G4double det_x[det_IDmax];
G4double det_y[det_IDmax];
G4double det_z[det_IDmax];
G4double det_kine[det_IDmax];
G4double det_VrtxKine[det_IDmax];
G4double det_VrtxX[det_IDmax];
G4double det_VrtxY[det_IDmax];
G4double det_VrtxZ[det_IDmax];
G4int det_VrtxVolID[det_IDmax];
G4int det_VrtxProcID[det_IDmax];
G4int det_VrtxTrackID[det_IDmax];
G4int det_VrtxParticleID[det_IDmax];
G4int det_VvvTrackSign[det_IDmax];
// Sort hits according to the time. Using std::map is convenient, because it sorts
// its entries according to the key (the first variable of the pair).
std::multimap<G4double,G4int> myHitTimeMapping; // "map" replaced by "multimap" (needed for radioactive decay,
// in which case times are huge and due to the limited rounding
// precision become equal --> map ignores the same "keys",
// multimap does not.
std::map<G4double,G4int>::iterator it;
for (G4int i=0; i<NbHits; i++) {
musrScintHit* aHit = (*scintCollection)[i];
G4double tmptime=aHit->GetGlobalTime();
// G4cout<<"Hit nr "<<i<<" at time="<<tmptime<<" with edep="<<aHit->GetEdep()/MeV
// <<" detID="<<myRootOutput->ConvertVolumeToID(aHit->GetLogVolName())<< G4endl;
myHitTimeMapping.insert ( std::pair<G4double,G4int>(tmptime,i) );
}
// Loop over all hits (which are sorted according to their time):
G4int nSignals=0;
for (it=myHitTimeMapping.begin(); it!=myHitTimeMapping.end(); it++) {
// G4cout << "Key:" << it->first;
// G4cout << " Value:" << it->second << "\n";
G4int ii = it->second; // ii is the index of the hits, which is sorted according to time
musrScintHit* aHit = (*scintCollection)[ii];
G4String aHitVolumeName = aHit->GetLogVolName();
G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName);
G4double aHitTime = aHit->GetGlobalTime();
G4int aHitTrackID = aHit->GetTrackID();
// Loop over all already defined signals and check whether the hit falls into any of them
G4bool signalAssigned=false;
for (G4int j=0; j<nSignals; j++) {
if ( (aHitVolumeID==det_ID[j]) && ((aHitTime-det_time_end[j])<mySignalSeparationTime) ) {
signalAssigned=true;
det_edep[j] += aHit->GetEdep();
det_nsteps[j]++;
det_length[j] += aHit->GetStepLength();
det_time_end[j] = aHitTime;
G4String aParticleName = aHit->GetParticleName();
if (aParticleName=="e-") {
det_edep_el[j] += aHit->GetEdep();
} else if (aParticleName=="e+") {
det_edep_pos[j] += aHit->GetEdep();
} else if (aParticleName=="gamma") {
det_edep_gam[j] += aHit->GetEdep();
} else if (aParticleName=="mu+") {
det_edep_mup[j] += aHit->GetEdep();
} else {
char message[200];
sprintf(message,"musrScintSD.cc::EndOfEvent(): untreated particle \"%s\" deposited energy.",aParticleName.c_str());
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
// Check whether the signals consits of more then just one hit, in which case make the track ID negative:
if (abs(det_VrtxTrackID[j])!=aHitTrackID) {
det_VrtxTrackID[j]=-1*abs(det_VrtxTrackID[j]);
if (boolIsVvvInfoRequested) {
if (det_VvvTrackSign[j]==1) {
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
if (!(myMusrSteppingAction->AreTracksCommingFromSameParent(aHitTrackID,abs(det_VrtxTrackID[j]),aHitVolumeName))) {det_VvvTrackSign[j]=-1;}
}
}
}
break;
}
}
if (!signalAssigned) { // The hit does not belong to any existing signal --> create a new signal.
// Check, whether the maximum number of signals was not exceeded:
if ( nSignals >= (det_IDmax-1) ) {
char message[200];
sprintf(message,"musrScintSD.cc::EndOfEvent(): number of signals exceeds maximal allowed value.");
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
else {
det_edep[nSignals] = aHit->GetEdep();
det_nsteps[nSignals] = 1;
det_length[nSignals] = aHit->GetStepLength();
det_ID[nSignals] = aHitVolumeID;
det_time_start[nSignals] = aHitTime;
det_time_end[nSignals] = aHitTime;
det_edep_el[nSignals] = 0;
det_edep_pos[nSignals] = 0;
det_edep_gam[nSignals] = 0;
det_edep_mup[nSignals] = 0;
G4String aParticleName = aHit->GetParticleName();
if (aParticleName=="e-") {
det_edep_el[nSignals] += aHit->GetEdep();
} else if (aParticleName=="e+") {
det_edep_pos[nSignals] += aHit->GetEdep();
} else if (aParticleName=="gamma") {
det_edep_gam[nSignals] += aHit->GetEdep();
} else if (aParticleName=="mu+") {
det_edep_mup[nSignals] += aHit->GetEdep();
} else {
char message[200];
sprintf(message,"musrScintSD.cc::EndOfEvent(): UNTREATED PARTICLE \"%s\" deposited energy.",aParticleName.c_str());
musrErrorMessage::GetInstance()->musrError(WARNING,message,true);
}
G4ThreeVector prePos = aHit->GetPrePos();
det_x[nSignals]=prePos.x();
det_y[nSignals]=prePos.y();
det_z[nSignals]=prePos.z();
det_kine[nSignals] = aHit->GetKineticEnergy();
det_VrtxKine[nSignals] = aHit->GetVertexKineticEnergy();
G4ThreeVector VrtxPos = aHit->GetVertexPosition();
det_VrtxX[nSignals] = VrtxPos.x();
det_VrtxY[nSignals] = VrtxPos.y();
det_VrtxZ[nSignals] = VrtxPos.z();
G4String logicalVolumeAtVertex = aHit->GetLogicalVolumeAtVertex();
det_VrtxVolID[nSignals] = myRootOutput->ConvertVolumeToID(logicalVolumeAtVertex);
G4String creatorProcessName = aHit->GetCreatorProcessName();
det_VrtxProcID[nSignals] = myRootOutput->ConvertProcessToID(creatorProcessName);
det_VrtxTrackID[nSignals] = aHit->GetTrackID();
det_VrtxParticleID[nSignals] = aHit->GetParticleID();
det_VvvTrackSign[nSignals] = 1;
nSignals++;
}
} // end of "if (!signalAssigned)"
} // end of the for loop over the hits
// Sort the signals according to the energy (in decreasing order)
std::map<G4double,G4int> mySignalMapping;
std::map<G4double,G4int>::iterator itt;
for (G4int i=0; i<nSignals; i++) {
mySignalMapping.insert ( std::pair<G4double,G4int>(-det_edep[i],i) );
}
// Write out the signals (sorted according to energy) to the musrRootOutput class:
G4int jj=-1;
for (itt=mySignalMapping.begin(); itt!=mySignalMapping.end(); itt++) {
jj++;
G4int ii = itt->second;
myRootOutput->SetDetectorInfo(jj,det_ID[ii],det_VrtxParticleID[ii],det_edep[ii],
det_edep_el[ii],det_edep_pos[ii],
det_edep_gam[ii],det_edep_mup[ii],det_nsteps[ii],det_length[ii],
det_time_start[ii],det_time_end[ii],det_x[ii],det_y[ii],det_z[ii],
det_kine[ii],det_VrtxKine[ii],det_VrtxX[ii],det_VrtxY[ii],det_VrtxZ[ii],
det_VrtxVolID[ii],det_VrtxProcID[ii],det_VrtxTrackID[ii] );
if (boolIsVvvInfoRequested) {
G4int oldTrackID = abs(det_VrtxTrackID[ii]);
musrSteppingAction* myMusrSteppingAction = musrSteppingAction::GetInstance();
G4int vvvParentTrackID= -1; G4int vvvPparticleID; G4double vvvKine; G4ThreeVector vvvPosition; G4String vvvLogVol; G4String vvvProcess;
G4int vvvLogVolID=-999;
G4bool oldTrackRetrievedOK;
do {
if (vvvParentTrackID>-1) {oldTrackID=vvvParentTrackID;}
oldTrackRetrievedOK = myMusrSteppingAction->GetInfoAboutOldTrack(oldTrackID,
vvvParentTrackID, vvvPparticleID, vvvKine, vvvPosition, vvvLogVol, vvvProcess);
if (oldTrackRetrievedOK) {
// G4cout<<"musrScintSD: Old Track seems to be achieved fine -> Lets save it to Root tree"<<G4endl;
vvvLogVolID=myRootOutput->ConvertVolumeToID(vvvLogVol);
}
else {
G4cout<<" oldTrackRetrievedOK is false"<<G4endl;
oldTrackID = -999;
vvvParentTrackID = -999;
vvvPparticleID = -999;
vvvKine = -999;
vvvPosition = G4ThreeVector(-999,-999,-999);
vvvLogVol = "undefined";
vvvProcess = "undefined";
}
} while (oldTrackRetrievedOK && (vvvLogVolID==det_ID[ii]));
if (oldTrackRetrievedOK) {
G4int vvvProcessID=myRootOutput->ConvertProcessToID(vvvProcess);
myRootOutput->SetDetectorInfoVvv(jj,vvvKine,vvvPosition.x(),vvvPosition.y(),vvvPosition.z(),
vvvLogVolID,vvvProcessID,oldTrackID*det_VvvTrackSign[ii],vvvPparticleID);
}
} //end "boolIsVvvInfoRequested"
}
} //end "if (NbHits>0)"
myRootOutput->FillEvent();
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

443
src/musrSteppingAction.cc Normal file
View File

@@ -0,0 +1,443 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrSteppingAction.hh"
#include "G4SteppingManager.hh"
#include "G4UnitsTable.hh"
#include "G4RunManager.hh" // needed for the event nr. comparison
#include "G4Run.hh" // ---------------||------------------
#include "G4MagneticField.hh" // needed for storing the magnetic field to the Root class
#include "G4FieldManager.hh" // ---------------||------------------
#include "G4TransportationManager.hh" // ---------------||------------------
#include "musrErrorMessage.hh"
#include "musrParameters.hh"
#include "F04GlobalField.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrSteppingAction::musrSteppingAction() {
pointer=this;
boolIsAnySpecialSaveVolumeDefined = false;
boolIsVvvInfoRequested = false;
boolMuonEventReweighting = false;
boolCalculateFieldIntegral = false;
myRootOutput = musrRootOutput::GetRootInstance();
if (myRootOutput == NULL) {
musrErrorMessage::GetInstance()->musrError(FATAL,
"musrSteppingAction::musrSteppingAction(): pointer to the musrRootOutput class not found! ==> EXECUTION STOPPED",true);
}
lastActualVolume="Unset";
}
musrSteppingAction::~musrSteppingAction() {
}
musrSteppingAction* musrSteppingAction::pointer=0;
musrSteppingAction* musrSteppingAction::GetInstance()
{
return pointer;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrSteppingAction::DoAtTheBeginningOfEvent() {
// G4cout<<"musrSteppingAction::DoAtTheBeginningOfEvent: at the beginning"<<G4endl;
radioactiveElectronAlreadySavedInThisEvent=false;
muAlreadyWasInTargetInThisEvent=false;
muAlreadyWasInM0InThisEvent=false;
muAlreadyWasInM1InThisEvent=false;
muAlreadyWasInM2InThisEvent=false;
myOldTracksMap.clear();
indexOfOldTrack = -1;
realTimeWhenThisEventStarted=time(0);
BxIntegral=0; ByIntegral=0; BzIntegral=0; BzIntegral1=0; BzIntegral2=0; BzIntegral3=0;
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
G4Track* aTrack = aStep->GetTrack();
// suspend the track if too many steps has already happened (relevant at high field)
if (aTrack->GetCurrentStepNumber()>100000) {
musrErrorMessage::GetInstance()->musrError(WARNING,
"musrSteppingAction: Current number of steps for the track > 100000 ==> TRACK KILLED",true);
G4double x=aStep->GetPostStepPoint()->GetPosition().x()/mm;
G4double y=aStep->GetPostStepPoint()->GetPosition().y()/mm;
G4double z=aStep->GetPostStepPoint()->GetPosition().z()/mm;
G4double E=aTrack->GetVertexKineticEnergy()/MeV;
myRootOutput->htest1->Fill(x,y);
myRootOutput->htest2->Fill(sqrt(x*x+y*y),z);
myRootOutput->htest3->Fill(E);
aTrack->SetTrackStatus(fStopAndKill);
}
// abort the event if it takes too long to finish (e.g. more than 60 seconds)
if ((time(0) - realTimeWhenThisEventStarted)>60) {
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4cout<<"musrSteppingAction: event "<<fRunManager->GetCurrentEvent()->GetEventID()
<<" aborted because calculation took already 60 seconds."<<G4endl;
musrErrorMessage::GetInstance()->musrError(WARNING,
"musrSteppingAction: event aborted because its calculation takes more than 60 seconds.",true);
//delete musrRootOutput* myRootOutput = musrRootOutput::GetRootInstance();
myRootOutput->SetEventWeight(0);
fRunManager->AbortEvent();
}
// Temporary fix to avoid crashes caused by particles with unphysically high energies
// (probably corrupted event?)
if ((aStep->GetPreStepPoint()->GetKineticEnergy()) > (1*GeV)) {
musrErrorMessage::GetInstance()->musrError(SERIOUS,
"musrSteppingAction: kinetic energy of a particle larger than 1GeV! STRANGE FOR muSR!",false);
G4RunManager* fRunManager = G4RunManager::GetRunManager();
G4cout<<" Event nr.:"<<fRunManager->GetCurrentEvent()->GetEventID()
<<", the particle \""<< aTrack->GetDynamicParticle()->GetDefinition()->GetParticleName()
<<"\" has energy of "<<(aStep->GetPreStepPoint()->GetKineticEnergy())/GeV<<" GeV."<<G4endl;
G4cout<<" Deleting the event!"<<G4endl;
G4cout.flush();
myRootOutput->SetEventWeight(0);
fRunManager->AbortEvent();
}
if (aTrack->GetDefinition()) {
G4ParticleDefinition* p_definition = aTrack->GetDynamicParticle()->GetDefinition();
G4String p_name = p_definition->GetParticleName();
// G4ProcessManager* p_manager = p_definition->GetProcessManager();
G4LogicalVolume* actualLogicalVolume = aTrack->GetVolume()->GetLogicalVolume();
G4String actualVolume = actualLogicalVolume->GetName();
// Delete track if the particle is in the "kill" volume.
// There is an example how to delete the track in example/novice/N04.
// It is done in a different way here, because the example/novice/N04 was not doing
// exactly what I wanted.
if((actualVolume(0,8)=="log_kill")||(actualVolume(0,8)=="log_Kill")) {
aTrack->SetTrackStatus(fStopAndKill); // suspend the track
}
if ((p_name=="nu_mu")||(p_name=="anti_nu_mu")||(p_name=="nu_e")||(p_name=="anti_nu_e")) {
//aTrack->SetTrackStatus(fStopAndKill); // suspend the tracks of neutrinos
}
// Save info about the old tracks, if the user wish to have Vvv info in the output Root Tree.
if (boolIsVvvInfoRequested) {
G4VPhysicalVolume* nextVolume = aTrack->GetNextVolume();
if (nextVolume!=NULL) {
if ((nextVolume->GetLogicalVolume()->GetSensitiveDetector()!=NULL)||(actualLogicalVolume->GetSensitiveDetector()!=NULL)) {
G4int trackID = aTrack->GetTrackID();
std::map<G4int,G4int>::iterator itr;
itr = myOldTracksMap.find(trackID);
if (itr==myOldTracksMap.end()) {
// track with this trackID has not been found in the map (has not been saved yet) ==> save it
indexOfOldTrack++;
myOldTracksMap.insert( std::pair<G4int,G4int>(trackID,indexOfOldTrack) );
if (indexOfOldTrack<maxNumberOfOldTracks) {
particleID_oldTrack[indexOfOldTrack] = aTrack->GetDefinition()->GetPDGEncoding();
parentTrackID_oldTrack[indexOfOldTrack] = aTrack->GetParentID();
vertexKine_oldTrack[indexOfOldTrack] = aTrack->GetVertexKineticEnergy();
vertexPosition_oldTrack[indexOfOldTrack] = aTrack->GetVertexPosition();
vertexLogVol_oldTrack[indexOfOldTrack] = aTrack->GetLogicalVolumeAtVertex()->GetName();
if ((aTrack->GetCreatorProcess())!=NULL) { vertexProcess_oldTrack[indexOfOldTrack] = aTrack->GetCreatorProcess()->GetProcessName();}
else { vertexProcess_oldTrack[indexOfOldTrack] = "initialParticle";}
}
else {
musrErrorMessage::GetInstance()->musrError(WARNING,
"musrSteppingAction: Maximum number of oldTracks reached ==> det_VvvXXX variables might be affected.",true);
}
}
}
}
}
// This are the data just for the radioactive decay (when using the radioactive source):
if ((musrParameters::boolG4GeneralParticleSource)) {
// &&(!radioactiveElectronAlreadySavedInThisEvent)) {
if (aTrack->GetTrackID() != 1 ){
if (aTrack->GetCreatorProcess()->GetProcessName() == "RadioactiveDecay") {
if (aTrack->GetDefinition()->GetParticleName()=="e-") {
if (aTrack->GetCurrentStepNumber()==1) {
G4double electron_kinetic_energy=aStep->GetPreStepPoint()->GetKineticEnergy();
myRootOutput->htest4->Fill(electron_kinetic_energy);
}
}
}
}
}
// Check if particle comes to the special volume
if (boolIsAnySpecialSaveVolumeDefined) {
// G4bool isFirstStepInVolume=aStep->IsFirstStepInVolume();
// This does not work!!! (aStep->IsFirstStepInVolume() is always zero.) I do not understand why!
G4bool isFirstStepInVolume=false;
if (actualVolume!=lastActualVolume) {
lastActualVolume=actualVolume;
isFirstStepInVolume=true;
}
if (isFirstStepInVolume) {
G4int tmpVolumeID=saveVolumeMapping[actualVolume];
if (tmpVolumeID!=0) {
G4int particle_id_save=p_definition->GetPDGEncoding();
G4double ke_save=aStep->GetPreStepPoint()->GetKineticEnergy();
G4double x_save=aStep->GetPreStepPoint()->GetPosition().x();
G4double y_save=aStep->GetPreStepPoint()->GetPosition().y();
G4double z_save=aStep->GetPreStepPoint()->GetPosition().z();
G4double px_save=aStep->GetPreStepPoint()->GetMomentum().x();
G4double py_save=aStep->GetPreStepPoint()->GetMomentum().y();
G4double pz_save=aStep->GetPreStepPoint()->GetMomentum().z();
myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,px_save,py_save,pz_save);
}
}
}
if ((p_name == "mu+")||(p_name == "mu-")||(p_name == "Mu")) {
// Store the information about the muon when it enters the target, M0, M1 or M2 for the fist time
// in a given event (i.e. the code has to be called just once during the event).
if ((actualVolume=="log_target")||(actualVolume=="log_Target")) {
if (!muAlreadyWasInTargetInThisEvent) {
muAlreadyWasInTargetInThisEvent=true;
myRootOutput->SetPolInTarget(aTrack->GetPolarization());
myRootOutput->SetTimeInTarget(aTrack->GetGlobalTime());
}
}
else if ((actualVolume=="log_M0")||(actualVolume=="log_m0")) {
if (!muAlreadyWasInM0InThisEvent) {
muAlreadyWasInM0InThisEvent=true;
myRootOutput->SetPolInM0(aTrack->GetPolarization());
myRootOutput->SetTimeInM0(aTrack->GetGlobalTime());
}
}
else if ((actualVolume=="log_M1")||(actualVolume=="log_m1")) {
if (!muAlreadyWasInM1InThisEvent) {
muAlreadyWasInM1InThisEvent=true;
myRootOutput->SetPolInM1(aTrack->GetPolarization());
myRootOutput->SetTimeInM1(aTrack->GetGlobalTime());
}
}
else if ((actualVolume=="log_M2")||(actualVolume=="log_m2")) {
if (!muAlreadyWasInM2InThisEvent) {
muAlreadyWasInM2InThisEvent=true;
myRootOutput->SetPolInM2(aTrack->GetPolarization());
myRootOutput->SetTimeInM2(aTrack->GetGlobalTime());
}
}
// Calculate the field integral along the muon path, if requested by the user.
// (2008.10.20. - idea of Robert is to calculate the distribution of the field integral for different muon paths
// to see what (delta_I/I) is still acceptable for hith field muSR. To calculate it properly, the user must
// force Geant to use very small step size.
if (boolCalculateFieldIntegral) {
if (F04GlobalField::Exists()) {
G4ThreeVector position_tmp = aStep->GetPostStepPoint()->GetPosition();
CoordinateForFieldIntegral[0] = position_tmp.x();
CoordinateForFieldIntegral[1] = position_tmp.y();
CoordinateForFieldIntegral[2] = position_tmp.z();
CoordinateForFieldIntegral[3] = aTrack->GetGlobalTime();
F04GlobalField::getObject() -> GetFieldValue(CoordinateForFieldIntegral,FieldForFieldIntegral);
G4double stepLength = aStep->GetStepLength();
// BxIntegral += stepLength;
// ByIntegral += stepLength;
// BzIntegral += stepLength;
BxIntegral += stepLength * FieldForFieldIntegral[0];
ByIntegral += stepLength * FieldForFieldIntegral[1];
BzIntegral += stepLength * FieldForFieldIntegral[2];
BzIntegral1 += ( position_tmp.z() - aStep->GetPreStepPoint()->GetPosition().z() ) * FieldForFieldIntegral[2];
BzIntegral2 += stepLength;
BzIntegral3 += ( position_tmp.z() - aStep->GetPreStepPoint()->GetPosition().z() );
// G4cout<<"BzIntegral="<<BzIntegral<<" stepLength="<<stepLength<<"FieldForFieldIntegral[2]="<<FieldForFieldIntegral[2]<<G4endl;
}
}
// Pick up process "DecayWithSpin":
const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep();
if (process!=NULL) {
G4String processName = process->GetProcessName();
if (processName=="DecayWithSpin") {
// std::cout<<"musrSteppingAction: DecayWithSpin"<<std::endl;
musrParameters::field_DecayWithSpin=true;
// Test whether the event reweighting is requeseted for this volume
// (i.e. user may request reweighting of events depending on the volume,
// in which the muon stops and decays).
if (boolMuonEventReweighting) {
G4int weight = volumeMuonWeightMapping[actualVolume];
if (weight!=0) {
G4double randomNumber = weight * G4UniformRand();
if (randomNumber < (weight-1.)) {
// G4cout<<"Event will be aborted"<<G4endl;
musrErrorMessage::GetInstance()->musrError(INFO,
"musrSteppingAction: event deleted because of the reweighting.",true);
G4RunManager* fRunManager = G4RunManager::GetRunManager();
weight=0;
fRunManager->AbortEvent();
}
// else {
// G4cout<<"Event will be reweighted"<<G4endl;
// }
myRootOutput->SetEventWeight(weight);
}
}
// Store the information about the decaying muon and store it in the Root tree
G4double timeOfDecay_tmp=aTrack->GetGlobalTime();
G4ThreeVector positionOfDecay_tmp = aStep->GetPostStepPoint()->GetPosition();
G4double BFieldAtOrigin[6] = {0.,0.,0.,0.,0.,0.};
G4double PointOfDecay[4] ={positionOfDecay_tmp.x(),positionOfDecay_tmp.y(),positionOfDecay_tmp.z(),timeOfDecay_tmp};
if (F04GlobalField::Exists()) {
F04GlobalField* myGlobalField = F04GlobalField::getObject();
myGlobalField->GetFieldValue(PointOfDecay,BFieldAtOrigin);
}
else {
G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
// G4cout<<"Debug 1"<<G4endl;
if (fMgr!=NULL) {
// G4cout<<"Debug 2"<<G4endl;
if(!fMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field
// G4cout<<"Debug 3"<<G4endl;
fMgr->GetDetectorField()->GetFieldValue(PointOfDecay,BFieldAtOrigin);
}
}
}
myRootOutput->SetDecayTime(timeOfDecay_tmp);
myRootOutput->SetDecayPolarisation(aTrack->GetPolarization());
myRootOutput->SetDecayPosition(positionOfDecay_tmp);
myRootOutput->SetDecayDetectorID(actualVolume);
myRootOutput->SetBField(BFieldAtOrigin);
if (boolCalculateFieldIntegral) {
myRootOutput->SetBFieldIntegral(BxIntegral,ByIntegral,BzIntegral,BzIntegral1,BzIntegral2,BzIntegral3);
}
// G4cout<<"================================================== BzIntegral = "<<BzIntegral<<G4endl;
// store the information about the emerging positron
G4TrackVector* secondary = fpSteppingManager->GetSecondary();
G4int n_secondaries= (*secondary).size();
for (G4int i=0; i<n_secondaries; i++) {
if ( ((*secondary)[i]->GetDefinition()->GetParticleName()) == "e+" ) {
myRootOutput->SetInitialPositronMomentum((*secondary)[i]->GetMomentum());
}
}
}
}
// G4ThreeVector position = aStep->GetPostStepPoint()->GetPosition();
// G4ThreeVector polarization=aTrack->GetDynamicParticle()->GetPolarization();
}
else { // particle is not muon
// Delete track if the particle is far away from the detector (i.e. in the "shield" volume).
// There is an example how to delete the track in example/novice/N04.
// It is done in a different way here, because the example/novice/N04 was not doing
// exactly what I wanted.
if ( ((musrParameters::killAllPositrons)&&(p_name == "e+")) ||
((musrParameters::killAllGammas)&&(p_name == "gamma")) ||
((musrParameters::killAllNeutrinos)&&((p_name == "nu_mu")||(p_name == "anti_nu_mu")||(p_name == "nu_e")||(p_name == "anti_nu_e"))) ){
aTrack->SetTrackStatus(fStopAndKill); // suspend the track
}
if((actualVolume(0,10)=="log_shield")||(actualVolume(0,10)=="log_Shield")) {
aTrack->SetTrackStatus(fStopAndKill); // suspend the track
}
}
}
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrSteppingAction::SetLogicalVolumeAsSpecialSaveVolume(G4String logicName, G4int volumeID) {
boolIsAnySpecialSaveVolumeDefined = true;
saveVolumeMapping[logicName]=volumeID;
}
void musrSteppingAction::SetVolumeForMuonEventReweighting(G4String logicName, G4int weight) {
boolMuonEventReweighting = true;
volumeMuonWeightMapping[logicName]=weight;
}
G4bool musrSteppingAction::GetInfoAboutOldTrack(G4int trackID, G4int& parentTrackID, G4int& particleID, G4double& vertexKine,
G4ThreeVector& vertexPosition, G4String& vertexLogVol, G4String& vertexProcess) {
// G4int ind = myOldTracksMap[trackID];
// G4cout<<"musrSteppingAction::GetInfoAboutOldTrack: trackID="<<trackID<<"\t myOldTracksMap[trackID]="<<myOldTracksMap[trackID]<<G4endl;
std::map<G4int,G4int>::iterator itr;
itr = myOldTracksMap.find(trackID);
if ( itr==myOldTracksMap.end() ) {
// if ((ind==0)||(ind>=maxNumberOfOldTracks)) {
char eMessage[200];
sprintf(eMessage,"musrSteppingAction::GetInfoAboutOldTrack: trackID not found in myOldTracksMap, det_VvvXXX variables might be affected");
musrErrorMessage::GetInstance()->musrError(WARNING,eMessage,false);
G4cout<<" Requested trackID="<<trackID<<G4endl;
// G4cout<<"Saved tracks:"<<G4endl;
// for (itr=myOldTracksMap.begin(); itr!=myOldTracksMap.end(); itr++) {
// G4cout<<"first="<<itr->first<<"\tsecond="<<itr->second<<G4endl;
// }
return false;
}
else {
G4int ind = itr->second;
if (ind>=maxNumberOfOldTracks) {
G4cout<<"musrSteppingAction::GetInfoAboutOldTrack: Problem! ind>maxNumberOfOldTracks! ("<<ind<<">"<<maxNumberOfOldTracks<<")"<<G4endl;
G4cout<<" itr->first = "<<itr->first<<", trackID = "<<trackID<<G4endl;
return false;
}
parentTrackID=parentTrackID_oldTrack[ind];
if (trackID==parentTrackID) {
G4cout<<"musrSteppingAction::GetInfoAboutOldTrack: Problem! trackID==parentTrackID! ("<<trackID<<"=="<<parentTrackID<<")"<<G4endl;
return false;
}
particleID=particleID_oldTrack[ind];
vertexKine=vertexKine_oldTrack[ind];
vertexPosition= vertexPosition_oldTrack[ind];
vertexLogVol=vertexLogVol_oldTrack[ind];
vertexProcess=vertexProcess_oldTrack[ind];
}
// G4cout<<"GetInfoAboutOldTrack: trackID="<<trackID<<"\t parentTrackID="<<parentTrackID<<"\t particleID="<<particleID;
// G4cout<<"\t vertexKine="<<vertexKine<<"\t vertexLogVol="<<vertexLogVol<<"\t vertexProcess="<<vertexProcess<<G4endl;
return true;
}
G4bool musrSteppingAction::AreTracksCommingFromSameParent(G4int trackID1, G4int trackID2, G4String volumeName){
// There are two tracks with different track IDs. This routine finds the parents of both of them,
// which were created outside logical volume "volumeID". If both tracks have the same parent, the
// functions returns "true".
std::map<G4int,G4int>::iterator itr;
G4int ind;
G4int track1;
G4int trID = trackID1;
do {
track1=trID;
itr = myOldTracksMap.find(trID);
if ( itr==myOldTracksMap.end() ) {
G4cout<<"musrSteppingAction::AreTracksCommingFromSameParent() Strange, trackID1 ="<<trackID1<<" not found"<<G4endl;
return false;
}
ind = itr->second;
trID = parentTrackID_oldTrack[ind];
} while (vertexLogVol_oldTrack[ind]==volumeName);
G4int track2;
trID = trackID2;
do {
track2=trID;
itr = myOldTracksMap.find(trID);
if ( itr==myOldTracksMap.end() ) {
G4cout<<"musrSteppingAction::AreTracksCommingFromSameParent() Strange, trackID2 ="<<trackID2<<" not found"<<G4endl;
return false;
}
ind = itr->second;
trID = parentTrackID_oldTrack[ind];
} while (vertexLogVol_oldTrack[ind]==volumeName);
if (track1==track2) {return true;}
//G4cout<<"track1="<<track1<<"\ttrack2="<<track2<<G4endl; return true;}
// G4cout<<"\t\t\t\ttrack1="<<track1<<"\ttrack2="<<track2<<G4endl;
return false;
}

151
src/musrSteppingVerbose.cc Normal file
View File

@@ -0,0 +1,151 @@
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
#include "musrSteppingVerbose.hh"
#include "G4SteppingManager.hh"
#include "G4UnitsTable.hh"
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrSteppingVerbose::musrSteppingVerbose()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
musrSteppingVerbose::~musrSteppingVerbose()
{}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrSteppingVerbose::StepInfo()
{
CopyState();
G4int prec = G4cout.precision(3);
if( verboseLevel >= 1 ){
if( verboseLevel >= 4 ) VerboseTrack();
if( verboseLevel >= 3 ){
G4cout << G4endl;
G4cout << std::setw( 5) << "#Step#" << " "
<< std::setw( 6) << "X" << " "
<< std::setw( 6) << "Y" << " "
<< std::setw( 6) << "Z" << " "
<< std::setw( 9) << "KineE" << " "
<< std::setw( 9) << "dEStep" << " "
<< std::setw(10) << "StepLeng"
<< std::setw(10) << "TrakLeng"
<< std::setw(10) << "Volume" << " "
<< std::setw(10) << "Process" << G4endl;
}
G4cout << std::setw(5) << fTrack->GetCurrentStepNumber() << " "
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().x(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy")
<< std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
<< std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length")
<< " ";
// if( fStepStatus != fWorldBoundary){
if( fTrack->GetNextVolume() != 0 ) {
G4cout << std::setw(10) << fTrack->GetVolume()->GetName();
} else {
G4cout << std::setw(10) << "OutOfWorld";
}
if(fStep->GetPostStepPoint()->GetProcessDefinedStep() != NULL){
G4cout << " "
<< std::setw(10) << fStep->GetPostStepPoint()->GetProcessDefinedStep()
->GetProcessName();
} else {
G4cout << " UserLimit";
}
G4cout << G4endl;
if( verboseLevel == 2 ){
G4int tN2ndariesTot = fN2ndariesAtRestDoIt +
fN2ndariesAlongStepDoIt +
fN2ndariesPostStepDoIt;
if(tN2ndariesTot>0){
G4cout << " :----- List of 2ndaries - "
<< "#SpawnInStep=" << std::setw(3) << tN2ndariesTot
<< "(Rest=" << std::setw(2) << fN2ndariesAtRestDoIt
<< ",Along=" << std::setw(2) << fN2ndariesAlongStepDoIt
<< ",Post=" << std::setw(2) << fN2ndariesPostStepDoIt
<< "), "
<< "#SpawnTotal=" << std::setw(3) << (*fSecondary).size()
<< " ---------------"
<< G4endl;
for(size_t lp1=(*fSecondary).size()-tN2ndariesTot;
lp1<(*fSecondary).size(); lp1++){
G4cout << " : "
<< std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().x(),"Length")
<< std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().y(),"Length")
<< std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetPosition().z(),"Length")
<< std::setw(6)
<< G4BestUnit((*fSecondary)[lp1]->GetKineticEnergy(),"Energy")
<< std::setw(10)
<< (*fSecondary)[lp1]->GetDefinition()->GetParticleName();
G4cout << G4endl;
}
G4cout << " :-----------------------------"
<< "----------------------------------"
<< "-- EndOf2ndaries Info ---------------"
<< G4endl;
}
}
}
G4cout.precision(prec);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrSteppingVerbose::TrackingStarted()
{
CopyState();
G4int prec = G4cout.precision(3);
if( verboseLevel > 0 ){
G4cout << std::setw( 5) << "Step#" << " "
<< std::setw( 6) << "X" << " "
<< std::setw( 6) << "Y" << " "
<< std::setw( 6) << "Z" << " "
<< std::setw( 9) << "KineE" << " "
<< std::setw( 9) << "dEStep" << " "
<< std::setw(10) << "StepLeng"
<< std::setw(10) << "TrakLeng"
<< std::setw(10) << "Volume" << " "
<< std::setw(10) << "Process" << G4endl;
G4cout << std::setw(5) << fTrack->GetCurrentStepNumber() << " "
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().x(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().y(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetPosition().z(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetKineticEnergy(),"Energy")
<< std::setw(6) << G4BestUnit(fStep->GetTotalEnergyDeposit(),"Energy")
<< std::setw(6) << G4BestUnit(fStep->GetStepLength(),"Length")
<< std::setw(6) << G4BestUnit(fTrack->GetTrackLength(),"Length")
<< " ";
if(fTrack->GetNextVolume()){
G4cout << std::setw(10) << fTrack->GetVolume()->GetName();
} else {
G4cout << std::setw(10) << "OutOfWorld";
}
G4cout << " initStep" << G4endl;
}
G4cout.precision(prec);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

596
src/musrTabulatedElementField.cc Normal file
View File

@@ -0,0 +1,596 @@
#include "musrTabulatedElementField.hh"
#include "musrParameters.hh"
#include "musrErrorMessage.hh"
#include "G4UnitsTable.hh"
musrTabulatedElementField::musrTabulatedElementField( const char* filename, const char* fldTableType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume),
ffieldValue(fieldValue)
{
// The following posibilities of the format of the field map are distinguieshed:
// 3DBOpera = 3D ... 3D field , magnetic, Opera format (Kamil-like)
// 3DE ... 3D field , electric, Toni-like (3DE WAS TESTED)
// 3DB ... 3D field , magnetic, Toni-like
// 2DBOpera = 2D .... 2D field, magnetic, Opera format (Kamil-like), X and Z components (2D WAS TESTED)
// 2DBOperaXY = 2D_OperaXY .... 2D field, magnetic, Opera format with (Kamil-like), X and Y components
// 2DE ... 2D field , electric, Toni-like (2DE WAS TESTED)
// 2DB ... 2D field , magnetic, Toni-like
G4double lenUnit = 1*m; // length unit of the field map grid coordinates
G4double fieldNormalisation = 1.; // Normalisation factor by which the the field map has to be multiplied
// in order to get 1T (in the case of magnetic field)
// "lenUnit" and "fieldNormalisation" are needed only if not specified
// inside the fieldmap file;
strcpy(fieldTableType,fldTableType);
fldType = 'B';
fUnit = "T";
fieUnit = tesla;
if (fieldTableType[2]=='E') {
fldType = 'E';
fUnit = "kV/mm";
fieUnit= kilovolt/mm;
}
fldDim = (fieldTableType[0]=='3') ? 3:2;
if (fldDim==2) ny=1;
G4cout << "\n-----------------------------------------------------------" << G4endl;
G4cout << " Field (of type "<<fieldTableType<<") set to "<< fieldValue/fieUnit << " "<< fUnit << G4endl;
G4cout << "\n ---> " "Reading the field grid from " << filename << " ... " << G4endl;
std::ifstream file( filename ); // Open the file for reading.
if (!(file.is_open())) {
G4cout << "Fieldmap file \""<< filename << "\" not opened (found) !!!"<<G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrTabulatedElementField: Field map file not found !",false);
}
char buffer[256];
G4bool boolMinimaAndMaximaDefinedInTheFile = false;
if (fldType=='E') G4cout<<" Electric field ";
if (fldType=='B') G4cout<<" Magnetic field ";
if ((strcmp(fieldTableType,"3DE")==0)||(strcmp(fieldTableType,"3DB")==0)) {
// File is in Toni Shiroka format:
// Read the number of arguments and decide filetype - 3 or 6 columns
G4cout << "3D, field-map file format by Toni Shiroka" << G4endl;
char lenUnitFromFile[50];
double fieldNormalisationFromFile;
file.getline(buffer,256);
int n_arg = sscanf (buffer,"%d %d %d %s %lf %lf %lf %lf %lf %lf %lf",
&nx, &ny, &nz, lenUnitFromFile, &fieldNormalisationFromFile,
&minimumx, &maximumx, &minimumy, &maximumy, &minimumz, &maximumz);
lenUnit = G4UnitDefinition::GetValueOf(lenUnitFromFile);
fieldNormalisation = fieldNormalisationFromFile;
if (n_arg==11) {
// The length unit and norm. factor have to be manually added to the field-map file!
G4cout << " ---> Assumed order (3 col.): "<<fldType<<"x, "<<fldType<<"y, "<<fldType<<"z"<<G4endl;
boolMinimaAndMaximaDefinedInTheFile = true;
minimumx = minimumx * lenUnit;
minimumy = minimumy * lenUnit;
minimumz = minimumz * lenUnit;
maximumx = maximumx * lenUnit;
maximumy = maximumy * lenUnit;
maximumz = maximumz * lenUnit;
}
else {
G4cout << " ---> Assumed order (6 col.): x, y, z, "<<fldType<<"x, "<<fldType<<"y, "<<fldType<<"z"<<G4endl;
}
// Ignore header information. All lines whose first character
// is '%' are considered to be part of the header.
do {
file.ignore(256, '\n');
} while (file.peek() == '%');
}
else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)) {
// OPERA format of the input file:
// Read table dimensions
lenUnit = 1*m;
fieldNormalisation = 1.;
G4cout << "3D, field-map file format from OPERA (Kamil)" << G4endl;
G4cout << " ---> Assumed order (7 col.): x, y, z, "<<fldType<<"x, "<<fldType<<"y, "<<fldType<<"z, Dummy"<<G4endl;
file.getline(buffer,256); // Skip the first empty line of the file
file >> nx >> ny >> nz; // Note dodgy order
// Ignore other header information
// The first line whose second character is '0' is considered to
// be the last line of the header.
do {
file.getline(buffer,256);
} while ( buffer[1]!='0');
}
else if ((strcmp(fieldTableType,"2DE")==0)||(strcmp(fieldTableType,"2DB")==0)||(strcmp(fieldTableType,"2DEf")==0)) {
// File is in Toni Shiroka format:
G4cout << "2D, field-map file format by Toni Shiroka" << G4endl;
G4cout << " ---> Assumed order (4 col.): r, z, "<<fldType<<"r, "<<fldType<<"z"<<G4endl;
char lenUnitFromFile[50];
file >> nx >> nz >> lenUnitFromFile >> fieldNormalisation;
lenUnit = G4UnitDefinition::GetValueOf(lenUnitFromFile);
// Ignore header information. All lines whose first character
// is '%' are considered to be part of the header.
do {
file.ignore(256, '\n');
} while (file.peek() == '%');
}
else if ((strcmp(fieldTableType,"2D_OperaXY")==0)||(strcmp(fieldTableType,"2DBOperaXY")==0)) {
int nDummy;
lenUnit = 1*cm;
fieldNormalisation = 0.00001;
G4cout << "2D, field-map file format from OPERA with X,Y components (Kamil)" << G4endl;
G4cout << " ---> Assumed order (6 col.): r, z, dummy, "<<fldType<<"r, "<<fldType<<"z, Dummy"<<G4endl;
file >> nx >> nz >> nDummy;
do {
file.getline(buffer,256);
} while ( buffer[1]!='0');
}
else if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)) {
int nDummy;
lenUnit = 1*cm;
fieldNormalisation = 0.00001;
G4cout << "2D, field-map file format from OPERA with X,Z components (Kamil)" << G4endl;
G4cout << " ---> Assumed order (6 col.): r, dummy, z, "<<fldType<<"r, "<<fldType<<"z, Dummy"<<G4endl;
file >> nx >> nDummy >> nz;
// G4cout << nx <<" "<< nDummy <<" "<< nz<<G4endl;
do {
file.getline(buffer,256);
} while ( buffer[1]!='0');
}
else {
G4cout << " musrTabulatedElementField::musrTabulatedElementField: Unknown field required!"
<< " ("<<fieldTableType<<")" << G4endl;
G4cout << " =====> S T O P " << G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrTabulatedElementField: Unknown field required!",false);
}
// SET UP STORAGE SPACE FOR THE TABLE
int ix, iy, iz;
if (fldDim==3) {
G4cout << " The grid consists of [" << nx << " x " << ny << " x " << nz << "] x, y, z values" << G4endl;
G4cout << " Field map length unit = " << lenUnit/mm<<" mm"<< G4endl;
G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl;
// Set up storage space for the table
xField.resize( nx );
yField.resize( nx );
zField.resize( nx );
for (ix=0; ix<nx; ix++) {
xField[ix].resize(ny);
yField[ix].resize(ny);
zField[ix].resize(ny);
for (iy=0; iy<ny; iy++) {
xField[ix][iy].resize(nz);
yField[ix][iy].resize(nz);
zField[ix][iy].resize(nz);
}
}
}
else if (fldDim==2) {
G4cout << " The grid consists of [" << nx << " x " << nz << "] R, z values" << G4endl;
G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl;
G4cout << " Field map length unit = " << lenUnit << G4endl;
if ((nx<2)||(nz<2)) {
char eMessage[200];
sprintf(eMessage,"musrTabulatedElementField(): Strange Field table! nx=%i, nz=%i",nx,nz);
musrErrorMessage::GetInstance()->musrError(WARNING,eMessage,false);
}
// Set up storage space for the table
xField2D.resize( nx );
zField2D.resize( nx );
for (ix=0; ix<nx; ix++) {
xField2D[ix].resize(nz);
zField2D[ix].resize(nz);
}
}
// Read in the data
double xval,yval,zval,bx,by,bz;
double permeability; // Not used in this example.
for (ix=0; ix<nx; ix++) {
for (iy=0; iy<ny; iy++) {
for (iz=0; iz<nz; iz++) {
if ((strcmp(fieldTableType,"3DE")==0)||(strcmp(fieldTableType,"3DB")==0)) {
if (boolMinimaAndMaximaDefinedInTheFile) {
file >> bx >> by >> bz; // Read ONLY field values
}
else {
file >> xval >> yval >> zval >> bx >> by >> bz;
}
}
else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)) {
file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
}
else if ((strcmp(fieldTableType,"2DE")==0)||(strcmp(fieldTableType,"2DB")==0)||(strcmp(fieldTableType,"2DEf")==0)) {
file >> xval >> zval >> bx >> bz;
}
else if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)) {
file >> xval >> yval >> zval >> bx >> bz >> permeability;
// G4cout<< xval <<" "<< yval <<" "<< zval <<" "<< bx <<" "<< bz <<G4endl;
}
else if (strcmp(fieldTableType,"2D_OperaXY")==0) {
file >> xval >> zval >> yval >> bx >> bz >> permeability;
}
else {
G4cout << " musrTabulatedElementField::musrTabulatedElementField: Undefined field required!"
<< " ("<<fieldTableType<<")" << G4endl;
G4cout << " =====> S T O P " << G4endl;
musrErrorMessage::GetInstance()->musrError(FATAL,"musrTabulatedElementField: Undefined field required!",false);
}
if (fldDim==3) { // 3D field
if ((!boolMinimaAndMaximaDefinedInTheFile) && ( ix==0 && iy==0 && iz==0 ) ) {
minimumx = xval * lenUnit;
minimumy = yval * lenUnit;
minimumz = zval * lenUnit;
}
xField[ix][iy][iz] = bx*fieldNormalisation;
yField[ix][iy][iz] = by*fieldNormalisation;
zField[ix][iy][iz] = bz*fieldNormalisation;
}
else { // 2D field
if ((!boolMinimaAndMaximaDefinedInTheFile) && ( ix==0 && iz==0 ) ) {
minimumx = xval * lenUnit;
minimumz = zval * lenUnit;
}
xField2D[ix][iz] = bx*fieldNormalisation;
zField2D[ix][iz] = bz*fieldNormalisation;
}
}
}
}
file.close();
if (!boolMinimaAndMaximaDefinedInTheFile) {
maximumx = xval * lenUnit;
maximumz = zval * lenUnit;
if (fldDim==3) maximumy = yval * lenUnit;
}
G4cout << " ---> ... reading of the field map finished." << G4endl;
if (fldDim==3) G4cout<<" ---> Min values of x,y,z: ";
else G4cout<<" ---> Min values of R,z: ";
G4cout << minimumx/cm << ", ";
if (fldDim==3) G4cout<< minimumy/cm << ", " ;
G4cout << minimumz/cm << " cm "<<G4endl;
if (fldDim==3) G4cout<<" ---> Max values of x,y,z: ";
else G4cout<<" ---> Max values of R,z: ";
G4cout << maximumx/cm << ", ";
if (fldDim==3) G4cout<< maximumy/cm << ", " ;
G4cout << maximumz/cm << " cm " << G4endl;
// Should really check that the limits are not the wrong way around.
G4bool reorderingDone = false;
if (maximumx < minimumx) {Invert("x"); reorderingDone=true;}
if (fldDim==3) {if (maximumy < minimumy) {Invert("y"); reorderingDone=true;} }
if (maximumz < minimumz) {Invert("z"); reorderingDone=true;}
if (reorderingDone) {
G4cout << "\n Reordering of the field grid was neccessary - after reordering:"<<G4endl;
if (fldDim==3) G4cout<<" ---> Min values of x,y,z: ";
else G4cout<<" ---> Min values of R,z: ";
G4cout << minimumx/cm << ", ";
if (fldDim==3) G4cout<< minimumy/cm << ", " ;
G4cout << minimumz/cm << " cm "<<G4endl;
if (fldDim==3) G4cout<<" ---> Max values of x,y,z: ";
else G4cout<<" ---> Max values of R,z: ";
G4cout << maximumx/cm << ", ";
if (fldDim==3) G4cout<< maximumy/cm << ", " ;
G4cout << maximumz/cm << " cm " << G4endl;
}
dx = maximumx - minimumx;
if (fldDim==3) dy = maximumy - minimumy;
dz = maximumz - minimumz;
if (fldDim==3) G4cout << "\n ---> Dif values x,y,z (range): ";
else G4cout << "\n ---> Dif values R,z (range): ";
G4cout << dx/cm << ", ";
if (fldDim==3) G4cout << dy/cm << ", ";
G4cout << dz/cm << " cm."<<G4endl;
G4cout << "-----------------------------------------------------------" << G4endl;
// Set maximum width, height and lenght of the field - very important! This
// dimensions are used to decide whether the "addFieldValue" method will be
// called at all for a given elementField.
if (fldDim==2) {
maximumWidth = 2*dx;
maximumHeight = 2*dx;
if ( (strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)||(strcmp(fieldTableType,"2D_OperaXY")) ) {
maximumLength = 2*dz;
}
else maximumLength = dz;
}
else {
maximumWidth = dx;
maximumHeight = dy;
maximumLength = dz;
}
}
void musrTabulatedElementField::addFieldValue(const G4double point[4],
G4double *field ) const
{
// G4cout<<"musrTabulatedElementField::addFieldValue"<<G4endl;
if (fldDim==2) addFieldValue2D(point,field);
else addFieldValue3D(point,field);
}
void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
G4double *field ) const
{
G4double B[3]; // Field value obtained from the field table
G4ThreeVector global(point[0],point[1],point[2]);
G4ThreeVector local;
local = global2local.TransformPoint(global);
double x = local.x();
double y = local.y();
double z = local.z();
// G4cout<<"Global points= "<<point[0]<<", "<<point[1]<<", "<<point[2]<<", Local point= "<<x<<", "<<y<<", "<<z<<G4endl;
// Check that the point is within the defined region
if ( x>minimumx && x<maximumx &&
y>minimumy && y<maximumy &&
z>minimumz && z<maximumz ) {
// Position of given point within region, normalized to the range
// [0,1]
double xfraction = (x - minimumx) / dx;
double yfraction = (y - minimumy) / dy;
double zfraction = (z - minimumz) / dz;
// Need addresses of these to pass to modf below.
// modf uses its second argument as an OUTPUT argument.
double xdindex, ydindex, zdindex;
// Position of the point within the cuboid defined by the
// nearest surrounding tabulated points
double xlocal = ( modf(xfraction*(nx-1), &xdindex));
double ylocal = ( modf(yfraction*(ny-1), &ydindex));
double zlocal = ( modf(zfraction*(nz-1), &zdindex));
// The indices of the nearest tabulated point whose coordinates
// are all less than those of the given point
int xindex = static_cast<int>(xdindex);
int yindex = static_cast<int>(ydindex);
int zindex = static_cast<int>(zdindex);
//cks The following check is necessary - even though xindex and zindex should never be out of range,
// it may happen (due to some rounding error ?). It is better to leave the check here.
if ((xindex<0)||(xindex>(nx-2))) {
std::cout<<"SERIOUS PROBLEM: xindex out of range! xindex="<<xindex<<" x="<<x<<" xfraction="<<xfraction<<std::endl;
if (xindex<0) xindex=0;
else xindex=nx-2;
}
if ((yindex<0)||(yindex>(ny-2))) {
std::cout<<"SERIOUS PROBLEM: yindex out of range! yindex="<<yindex<<" y="<<y<<" yfraction="<<yfraction<<std::endl;
if (yindex<0) yindex=0;
else yindex=ny-2;
}
if ((zindex<0)||(zindex>(nz-2))) {
std::cout<<"SERIOUS PROBLEM: zindex out of range! zindex="<<zindex<<" z="<<z<<" zfraction="<<zfraction<<std::endl;
if (zindex<0) zindex=0;
else zindex=nz-2;
}
// Full 3-dimensional version
B[0] =
xField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
xField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
xField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
xField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
xField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
xField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
xField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
xField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
B[1] =
yField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
yField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
yField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
yField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
yField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
yField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
yField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
yField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
B[2] =
zField[xindex ][yindex ][zindex ] * (1-xlocal) * (1-ylocal) * (1-zlocal) +
zField[xindex ][yindex ][zindex+1] * (1-xlocal) * (1-ylocal) * zlocal +
zField[xindex ][yindex+1][zindex ] * (1-xlocal) * ylocal * (1-zlocal) +
zField[xindex ][yindex+1][zindex+1] * (1-xlocal) * ylocal * zlocal +
zField[xindex+1][yindex ][zindex ] * xlocal * (1-ylocal) * (1-zlocal) +
zField[xindex+1][yindex ][zindex+1] * xlocal * (1-ylocal) * zlocal +
zField[xindex+1][yindex+1][zindex ] * xlocal * ylocal * (1-zlocal) +
zField[xindex+1][yindex+1][zindex+1] * xlocal * ylocal * zlocal ;
B[0] *= ffieldValue;
B[1] *= ffieldValue;
B[2] *= ffieldValue;
G4ThreeVector finalField(B[0],B[1],B[2]);
finalField = global2local.Inverse().TransformAxis(finalField);
if (fldType == 'E') {
field[3] += finalField.x();
field[4] += finalField.y();
field[5] += finalField.z();
}
else {
field[0] += finalField.x();
field[1] += finalField.y();
field[2] += finalField.z();
}
}
// G4cout<<"Kamil: Field: ("<<field[0]/tesla<<","<<field[1]/tesla<<","<<field[2]/tesla<<")"<<G4endl;
}
void musrTabulatedElementField::addFieldValue2D(const G4double point[4],
G4double *field ) const
{
G4double B[3]; // Field value obtained from the field table
G4ThreeVector global(point[0],point[1],point[2]);
G4ThreeVector local;
local = global2local.TransformPoint(global);
double x, z, z_sign;
if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)||
(strcmp(fieldTableType,"2D_OperaXY"))||(strcmp(fieldTableType,"2DEf")==0)) {
// Field is defined in just positive range of z; i.e. it is expected to be "symmetric"
// and the field for negative z is calculated from the positive z half.
x = sqrt(local.x()*local.x()+local.y()*local.y());
z = fabs(local.z());
z_sign = (local.z()>0) ? 1.:-1.;
}
else {
// Field is defined along the whole range of the z axis (i.e. asymmetric field is expected)
x = sqrt(local.x()*local.x()+local.y()*local.y());
z = local.z();
z_sign = 1;
}
// Check that the point is within the defined region
if ( x<maximumx && z<maximumz ) {
// if (evNr>evNrKriz) std::cout<<"bol som tu"<<std::endl;
// Position of given point within region, normalized to the range
// [0,1]
double xfraction = (x - minimumx) / dx;
double zfraction = (z - minimumz) / dz;
// Need addresses of these to pass to modf below.
// modf uses its second argument as an OUTPUT argument.
double xdindex, zdindex;
// Position of the point within the cuboid defined by the
// nearest surrounding tabulated points
double xlocal = ( modf(xfraction*(nx-1), &xdindex));
double zlocal = ( modf(zfraction*(nz-1), &zdindex));
// The indices of the nearest tabulated point whose coordinates
// are all less than those of the given point
int xindex = static_cast<int>(xdindex);
int zindex = static_cast<int>(zdindex);
//cks The following check is necessary - even though xindex and zindex should never be out of range,
// it may happen (due to some rounding error ?). It is better to leave the check here.
if ((xindex<0)||(xindex>(nx-2))) {
std::cout<<"SERIOUS PROBLEM: xindex out of range! xindex="<<xindex<<" x="<<x<<" xfraction="<<xfraction<<std::endl;
if (xindex<0) xindex=0;
else xindex=nx-2;
}
if ((zindex<0)||(zindex>(nz-2))) {
std::cout<<"SERIOUS PROBLEM: zindex out of range! zindex="<<zindex<<" z="<<z<<" zfraction="<<zfraction<<std::endl;
if (zindex<0) zindex=0;
else zindex=nz-2;
}
// G4cout<<"xField2D["<<xindex<<"]["<<zindex<<"]="<<xField2D[xindex ][zindex ]<<G4endl;
// G4cout<<"zField2D["<<xindex<<"]["<<zindex<<"]="<<zField2D[xindex ][zindex ]<<G4endl;
// Interpolate between the neighbouring points
double Bfield_R =
xField2D[xindex ][zindex ] * (1-xlocal) * (1-zlocal) +
xField2D[xindex ][zindex+1] * (1-xlocal) * zlocal +
xField2D[xindex+1][zindex ] * xlocal * (1-zlocal) +
xField2D[xindex+1][zindex+1] * xlocal * zlocal ;
B[0] = (x>0) ? Bfield_R * (local.x() /x) : 0.;
B[1] = (x>0) ? Bfield_R * (local.y() /x) : 0.;
B[2] =
zField2D[xindex ][zindex ] * (1-xlocal) * (1-zlocal) +
zField2D[xindex ][zindex+1] * (1-xlocal) * zlocal +
zField2D[xindex+1][zindex ] * xlocal * (1-zlocal) +
zField2D[xindex+1][zindex+1] * xlocal * zlocal ;
if (fldType == 'E') { // Electric field
B[0] *= ffieldValue;
B[1] *= ffieldValue;
B[2] *= ffieldValue * z_sign;
}
else { // Magnetic field
B[0] *= ffieldValue * z_sign;
B[1] *= ffieldValue * z_sign;
B[2] *= ffieldValue;
}
G4ThreeVector finalField(B[0],B[1],B[2]);
finalField = global2local.Inverse().TransformAxis(finalField);
if (fldType == 'E') {
field[3] += finalField.x();
field[4] += finalField.y();
field[5] += finalField.z();
}
else {
field[0] += finalField.x();
field[1] += finalField.y();
field[2] += finalField.z();
}
// G4cout<<"F= "<<field[0]<<" "<<field[1]<<" "<<field[2]<<" "<<field[3]<<" "<<field[4]<<" "<<field[5]<<G4endl;
}
}
G4double musrTabulatedElementField::GetNominalFieldValue() {
return ffieldValue;
}
void musrTabulatedElementField::SetNominalFieldValue(G4double newFieldValue) {
// // Rescale the magnetic field for a new value of the magnetic field
ffieldValue=newFieldValue;
G4cout<<"musrTabulatedElementField.cc: ffieldValue changed to="<< ffieldValue/fieUnit<<" "<<fUnit<<G4endl;
}
void musrTabulatedElementField::Invert(const char* indexToInvert) {
// This function inverts the indexes of the field table for a given axis (x or z).
// It should be called in the case when the x or z coordinate in the initial
// field table is ordered in the decreasing order.
std::vector< std::vector< std::vector< double > > > xFieldTemp(xField);
std::vector< std::vector< std::vector< double > > > yFieldTemp(yField);
std::vector< std::vector< std::vector< double > > > zFieldTemp(zField);
G4bool invertX=false;
G4bool invertY=false;
G4bool invertZ=false;
G4cout<<"Check that the musrTabulatedElementField::Invert() function works properly!"<<G4endl;
G4cout<<"It has not been tested yet!"<<G4endl;
if (strcmp(indexToInvert,"x")==0) {invertX=true; std::swap(maximumx,minimumx);}
if (strcmp(indexToInvert,"y")==0) {invertY=true; std::swap(maximumx,minimumx);}
if (strcmp(indexToInvert,"z")==0) {invertZ=true; std::swap(maximumz,minimumz);}
for (int ix=0; ix<nx; ix++) {
for (int iy=0; iy<ny; iy++) {
for (int iz=0; iz<nz; iz++) {
if (invertX) {
xField[ix][iy][iz] = xFieldTemp[nx-1-ix][iy][iz];
yField[ix][iy][iz] = yFieldTemp[nx-1-ix][iy][iz];
zField[ix][iy][iz] = zFieldTemp[nx-1-ix][iy][iz];
}
else if(invertY) {
xField[ix][iy][iz] = xFieldTemp[ix][ny-1-iy][iz];
yField[ix][iy][iz] = yFieldTemp[ix][ny-1-iy][iz];
zField[ix][iy][iz] = zFieldTemp[ix][ny-1-iy][iz];
}
else if(invertZ) {
xField[ix][iy][iz] = xFieldTemp[ix][iy][nz-1-iz];
yField[ix][iy][iz] = yFieldTemp[ix][iy][nz-1-iz];
zField[ix][iy][iz] = zFieldTemp[ix][iy][nz-1-iz];
}
}
}
}
}

91
src/musrUniformField.cc Normal file
View File

@@ -0,0 +1,91 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
#include "globals.hh"
#include "G4GeometryManager.hh"
#include "musrUniformField.hh"
musrUniformField::musrUniformField(G4double EMF[6], G4double half_X, G4double half_Y, G4double half_Z, G4LogicalVolume* lv, G4ThreeVector c)
: F04ElementField(c, lv) {
// EMF[6] ... constant vector of the field: 3 components magnetic + 3 components electric field
// half_X, half_Y, half_Z ... half dimenstions of the box, within which the field is defined
for (int i = 0; i < 6; i++){
EMfield[i] = EMF[i];
}
fieldLength = 2.*half_Z;
fieldWidth = 2.*half_X;
fieldHeight = 2.*half_Y;
G4cout << "\n-----------------------------------------------------------"
<< "\n Uniform electromagnetic field"
<< G4endl;
G4String volName = lv->GetName().substr(4);
G4cout << "\n ---> EM field in volume " << volName << " set to:" << G4endl;
printf (" B = (%0.3g, %0.3g, %0.3g) T, E = (%0.3g, %0.3g, %0.3g) kV/mm\n",
EMF[0]/tesla, EMF[1]/tesla, EMF[2]/tesla,
EMF[3]/(kilovolt/mm), EMF[4]/(kilovolt/mm), EMF[5]/(kilovolt/mm));
G4cout << "-----------------------------------------------------------" << G4endl;
}
void musrUniformField::addFieldValue(const G4double point[4],
G4double field[6]) const
{
G4ThreeVector global(point[0],point[1],point[2]);
G4ThreeVector local;
local = global2local.TransformPoint(global);
if (isOutside(local)) return;
G4ThreeVector B(EMfield[0],EMfield[1],EMfield[2]);
G4ThreeVector E(EMfield[3],EMfield[4],EMfield[5]);
B = global2local.Inverse().TransformAxis(B);
E = global2local.Inverse().TransformAxis(E);
field[0] += B[0];
field[1] += B[1];
field[2] += B[2];
field[3] += E[0];
field[4] += E[1];
field[5] += E[2];
//printf (" EM field components: B = (%0.3g, %0.3g, %0.3g) T, E = (%0.3g, %0.3g, %0.3g) kV/mm\n",
//field[0]/tesla, field[1]/tesla, field[2]/tesla,
//field[3]/(kilovolt/mm), field[4]/(kilovolt/mm), field[5]/(kilovolt/mm));
}
G4double musrUniformField::GetNominalFieldValue() {
G4double val = std::sqrt(EMfield[0]*EMfield[0]+EMfield[1]*EMfield[1]+EMfield[2]*EMfield[2]+
EMfield[3]*EMfield[3]+EMfield[4]*EMfield[4]+EMfield[5]*EMfield[5]);
return val;
}
void musrUniformField::SetNominalFieldValue(G4double newFieldValue){
EMfield[0] *= newFieldValue;
EMfield[1] *= newFieldValue;
EMfield[2] *= newFieldValue;
EMfield[3] *= newFieldValue;
EMfield[4] *= newFieldValue;
EMfield[5] *= newFieldValue;
}
G4bool musrUniformField::isOutside(G4ThreeVector& local) const
{
return (std::fabs(local.z()) > fieldLength/2.0 || std::fabs(local.x()) > fieldWidth/2.0 || std::fabs(local.y()) > fieldHeight/2.0);
}
G4bool musrUniformField::isWithin(G4ThreeVector& local) const
{
return (std::fabs(local.z()) < fieldLength/2.0 && std::fabs(local.x()) < fieldWidth/2.0 && std::fabs(local.y()) < fieldHeight/2.0);
}

107
src/yields.cc Normal file
View File

@@ -0,0 +1,107 @@
// Geant4 simulation for MuSR
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
// DATE : 2008-05
//
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
// Muonium yields as a function of initial mu+ energies.
// The method GetYields is used by MuFormation.
// Id : yields.cc, v 1.1
// Author: Taofiq PARAISO, T. Shiroka
// Date : 2007-12
// Notes : First implemented in Fortran by A. Hofer
// C++ conversion by T.K. Paraiso 04-2005
// Slight modifications by T. Shiroka
//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
#include "yields.hh"
#include <iomanip>
#include <fstream>
#include <iostream>
#include <stdlib.h>
/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
The Muonium Yield function as well as the parameters are taken from:
M. Gonin, R. Kallenbach, P. Bochsler: "Charge exchange of hydrogen atoms
in carbon foils at 0.4 - 120 keV", Rev.Sci.Instrum. 65 (3), March 1994
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
Yields:: Yields(){;}
Yields::~Yields(){;}
void Yields::GetYields(
double E, // kinetic energy in keV
double mass, // mass in keV/c**2
double yvector[3]) // pointer to the yields table
{
// Parameter NAMES for the muonium yield function
double a_zero, a_minus;
double k_Fermi, k_zero, k_minus;
double two_k_Fermi;
double k_Fermi_Quad, k_zero_Quad, k_minus_Quad;
double vc_minus, vc_plus, v_Bohr, v_rel;
// Parameter VALUES for the muonium yield function
a_zero = 0.953;
a_minus = 0.029;
k_Fermi = 1.178; // [v_Bohr]
k_Fermi_Quad = k_Fermi * k_Fermi;
two_k_Fermi = 2. * k_Fermi;
k_zero = 0.991*k_Fermi; // [v_Bohr]
k_zero_Quad = k_zero * k_zero;
k_minus = 0.989*k_Fermi; // [v_Bohr]
k_minus_Quad = k_minus * k_minus;
vc_minus = 0.284;
vc_plus = 0.193; // [v_Bohr]
v_Bohr = 7.2974E-3; // [c]
// std::cout<<"E = "<< E <<std::endl;
// Abort in case of negative energies ---------------------------
if (E < 0)
{
std::cout<< "Error in method ''Yields'':" <<std::endl;
std::cout<< "E = "<< E <<" < 0!" <<std::endl;
std::cout<< "-> ABORTED!" <<std::endl;
return;
}
//---------------------------------------------------------------
// Definition of variables (aux_n are some auxiliary variables)
// Calculate energy in (classical) terms of speed (in units of v_Bohr):
v_rel = sqrt(2.*E/mass)/ v_Bohr;
aux1 = v_rel*v_rel;
aux2 = two_k_Fermi*v_rel;
Q_zero = 1. + (k_zero_Quad - k_Fermi_Quad - aux1) / aux2;
Q_minus = 1. + (k_minus_Quad - k_Fermi_Quad - aux1) / aux2;
aux1 = a_zero * Q_zero;
aux2 = a_minus * Q_minus;
aux3 = (1.-Q_zero)*(1.-Q_minus);
D = aux1*(aux2 + (1.-Q_minus)) + aux3;
Yield_minus = aux1*aux2 / D;
Yield_plus = aux3 / D;
Yield_minus = Yield_minus* exp(-vc_minus/v_rel);
Yield_plus = Yield_plus * exp(-vc_plus /v_rel);
if(Yield_minus > exp(-vc_minus/v_rel)) Yield_minus=exp(-vc_minus/v_rel);
if(Yield_plus > exp(-vc_plus/v_rel)) Yield_plus=exp(-vc_plus/v_rel);
Yield_zero = 1. - (Yield_minus + Yield_plus);
yvector[0]=Yield_plus;
yvector[1]=Yield_zero;
yvector[2]=Yield_minus;
// std::cout<<"Y+ : "<< Yield_plus << std::endl;
// std::cout<<"Y0 : "<< Yield_zero << std::endl;
}