193 lines
7.3 KiB
C++
193 lines
7.3 KiB
C++
// Geant4 simulation for MuSR
|
|
// AUTHOR: Toni SHIROKA, Paul Scherrer Institut, PSI
|
|
// DATE : 2008-05
|
|
//
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
#include "lem4EventAction.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 "lem4MagneticField.hh"
|
|
#include "lem4TabulatedField3D.hh"
|
|
#include "lem4TabulatedField2D.hh"
|
|
#include "lem4RootOutput.hh"
|
|
#include "lem4ErrorMessage.hh"
|
|
#include "lem4SteppingAction.hh"
|
|
#include "F04GlobalField.hh"
|
|
//#include "F04ElementField.hh"
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
G4bool lem4EventAction::setRandomNrSeedAccordingEventNr=0;
|
|
G4bool lem4EventAction::setRandomNrSeedFromFile=0;
|
|
G4int lem4EventAction::nHowOftenToPrintEvent=10000;
|
|
//vector<int> * lem4EventAction::RandomNrInitialisers=NULL;
|
|
|
|
//long lem4EventAction::myEventNr=0;
|
|
|
|
lem4EventAction::lem4EventAction() {
|
|
pointer=this;
|
|
fieldValueStart=0;
|
|
pointerToSeedVector = new vector<int>;
|
|
timeDependentField=false;
|
|
lastFieldValue=-10000*tesla;
|
|
pointerToMusrUniformField=NULL;
|
|
pointerToTabulatedField3D=NULL;
|
|
pointerToTabulatedField2D=NULL;
|
|
latestEventNr=-1;
|
|
}
|
|
lem4EventAction* lem4EventAction::pointer=0;
|
|
lem4EventAction* lem4EventAction::GetInstance() {
|
|
return pointer;
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
lem4EventAction::~lem4EventAction()
|
|
{
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void lem4EventAction::BeginOfEventAction(const G4Event* evt) {
|
|
// test error
|
|
// lem4ErrorMessage::GetInstance()->lem4Error(SERIOUS,"test error",true);
|
|
//
|
|
// G4cout<<"lem4EventAction::BeginOfEventAction: KAMIL"<<G4endl;
|
|
lem4SteppingAction::GetInstance()->DoAtTheBeginningOfEvent();
|
|
|
|
long thisEventNr = (long) (evt->GetEventID());
|
|
|
|
if (F04GlobalField::Exists()) {
|
|
// if (F04GlobalField::getObject() -> DoesAnyFieldValueNeedsToBeChanged(thisEventNr)) {
|
|
// G4cout<<"We should check each Element Field Object whether its field needs to be changed:"<<G4endl;
|
|
// Loop over all Element Fields and tell them to change their nominal field value, if that was foreseen.
|
|
|
|
// F04ElementField
|
|
// }
|
|
F04GlobalField::getObject() -> CheckWhetherAnyNominalFieldValueNeedsToBeChanged(thisEventNr);
|
|
}
|
|
|
|
latestEventNr = thisEventNr;
|
|
G4double actualFieldValue;
|
|
if (timeDependentField) {
|
|
// actualFieldValue=fieldAtEnd-fieldAtBeginning
|
|
G4int i=int(double(thisEventNr)/double(maxEventNr)*fieldNrOfSteps); // i ... nr of actual step in the field
|
|
actualFieldValue=fieldValueStart+fieldStep*i;
|
|
if (actualFieldValue!=lastFieldValue) {
|
|
lastFieldValue=actualFieldValue;
|
|
// G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
|
|
// // fieldMgr->SetFieldValue(actualFieldValue); //did not work
|
|
// const G4Field* pointerToField = NULL;
|
|
// pointerToField=fieldMgr->GetDetectorField();
|
|
// pointerToField->GetFieldValue(0,0,0,0,
|
|
if (pointerToMusrUniformField) {
|
|
pointerToMusrUniformField->SetFieldValue(actualFieldValue);
|
|
// pointerToField->SetFieldValue(actualFieldValue); does not work
|
|
G4cout<<"Event "<<thisEventNr<<": Uniform field value changed to ";
|
|
G4cout<<(actualFieldValue/tesla)<<"Tesla."<<G4endl;
|
|
}
|
|
else if (pointerToTabulatedField3D) {
|
|
pointerToTabulatedField3D->SetFieldValue(actualFieldValue);
|
|
G4cout<<"Event "<<thisEventNr<<": Tabulated Field (3D) value changed to ";
|
|
G4cout<<(actualFieldValue/tesla)<<"Tesla."<<G4endl;
|
|
}
|
|
else if (pointerToTabulatedField2D) {
|
|
pointerToTabulatedField2D->SetFieldValue(actualFieldValue);
|
|
G4cout<<"Event "<<thisEventNr<<": Tabulated Field (2D) value changed to ";
|
|
G4cout<<(actualFieldValue/tesla)<<"Tesla."<<G4endl;
|
|
}
|
|
}
|
|
}
|
|
else {actualFieldValue=fieldValueStart;}
|
|
lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance();
|
|
myRootOutput->SetFieldValue(actualFieldValue);
|
|
|
|
// if (lem4DetectorMessenger::setRandomNrSeedAccordingEventNr) {
|
|
if (setRandomNrSeedFromFile) {
|
|
// G4cout<<"RandomNrInitialisers.size()="<<RandomNrInitialisers->size()<<G4endl;
|
|
if (thisEventNr < (long) pointerToSeedVector->size()) {
|
|
G4cout <<"lem4EventAction.cc: seed will be set to="<< pointerToSeedVector->at(thisEventNr)<<G4endl;
|
|
CLHEP::HepRandom::setTheSeed(pointerToSeedVector->at(thisEventNr));
|
|
}
|
|
}
|
|
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(thisEventNr);
|
|
// G4cout <<" getTheSeed="<<CLHEP::HepRandom::getTheSeed()<< G4endl;
|
|
CLHEP::RandGauss::setFlag(false);
|
|
}
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
void lem4EventAction::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 and thisEventNr%nHowOftenToPrintEvent == 0) {
|
|
time_t curr=time(0);
|
|
//char * ctime(const time_t * tp);
|
|
G4cout << ">>> Event " << evt->GetEventID() <<" at "<< 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);
|
|
}
|
|
}
|
|
}
|
|
|
|
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
|
|
|
|
vector<int> * lem4EventAction::pointerToSeedVector=NULL;
|
|
vector<int> * lem4EventAction::GetPointerToSeedVector() {
|
|
return pointerToSeedVector;
|
|
}
|
|
|
|
|
|
void lem4EventAction::SetTimeDependentField(G4bool setFieldToBeTimeDependend, G4double initialField,
|
|
G4double finalField, G4int nrOfSteps) {
|
|
timeDependentField = setFieldToBeTimeDependend;
|
|
fieldValueStart = initialField;
|
|
fieldValueEnd = finalField;
|
|
fieldNrOfSteps = nrOfSteps;
|
|
fieldStep = (finalField-initialField)/(nrOfSteps-1);
|
|
G4cout<<"---&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&---"<<G4endl;
|
|
}
|