diff --git a/geant4/LEMuSR/src/F04ElementField.cc b/geant4/LEMuSR/src/F04ElementField.cc new file mode 100644 index 0000000..28e81d8 --- /dev/null +++ b/geant4/LEMuSR/src/F04ElementField.cc @@ -0,0 +1,207 @@ +// +// ******************************************************************** +// * 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 "lem4Parameters.hh" +#include "lem4ErrorMessage.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: "<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="<LocateGlobalPointAndSetup(center,0,false); + + G4TouchableHistoryHandle fTouchable = aNavigator->CreateTouchableHistoryHandle(); + + G4int depth = fTouchable->GetHistoryDepth(); + for (G4int i = 0; iGetVolume()->GetLogicalVolume() == lvolume)break; + fTouchable->MoveUpHistory(); + } + + // Check that the point of origin of the field matches that of 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(): Origin of the field outside the respective logical volume!! \"%s\".", + volumeName.c_str()); + /// NOTE: Cannot apply fields to non simply connected bodies, as e.g. rings, tori, etc. ! + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + + // G4cout<<"+!+!+! global2local VOLUME NAME: "<GetVolume()->GetName()<GetHistory()->GetTopTransform(); + G4ThreeVector local_center = global2local.TransformPoint(center); + + G4cout<<"F04ElementField: " << elementFieldName << " in \""< "< "< 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(lem4Parameters::nrOfEventsToBeGenerated)/(nrOfSteps); + G4double fieldStep = (finalField-initialField)/(nrOfSteps-1); + // G4cout<<"lem4Parameters::nrOfEventsToBeGenerated="<::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::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 "<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)); + } +} diff --git a/geant4/LEMuSR/src/F04GlobalField.cc b/geant4/LEMuSR/src/F04GlobalField.cc new file mode 100644 index 0000000..009c274 --- /dev/null +++ b/geant4/LEMuSR/src/F04GlobalField.cc @@ -0,0 +1,313 @@ +// +// ******************************************************************** +// * 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 + +#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 "lem4RootOutput.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(); +} + +F04GlobalField::~F04GlobalField() +{ + 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; + + G4cout << G4endl; + 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 + + 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; iisInBoundingBox(point)) { + p->addFieldValue(point,field); + } + } + +} + + // cks NOT SURE WHETHER THE FOLLOWING IS STILL NEEDED, HOWEVER ADDED HERE JUST + // FOR THE CASE: + // Set some small field if field is almost zero (to avoid internal problems of Geant). + ///if (sqrt(field[0]*field[0]+field[1]*field[1]+field[2]*field[2])<0.00001*tesla) { + /// field[2] = 0.00001*tesla; + ///} + // csk + + +void F04GlobalField::clear() +{ + if (fields) { + if (fields->size()>0) { + FieldList::iterator i; + 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 localChangeFieldInStepsMap = fp[i] ->GetEventNrDependentField(); + std::map::iterator it; + for ( it=localChangeFieldInStepsMap.begin() ; it != localChangeFieldInStepsMap.end(); it++ ) { + G4int eventNr = it->first; + G4cout<<"globalChangeFieldInStepsMap set for event number "<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:"<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(); + lem4RootOutput::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 EMfield[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,EMfield); + // [MV/m] = [10^3 kV / 10^3 mm] = [kV/mm]. But MV and mm are the BASE units in Geant4 => Fact. = 0.001 + ///printf (" EM field value at (%.f, %.f, %.f) mm is: B = (%.f, %.f, %.f) T, E = (%.10f, %.10f, %.10f) kV/mm\n", + printf (" EM field value at (%.f, %.f, %.f) mm is: B = (%0.3g, %0.3g, %0.3g) T, E = (%0.3g, %0.3g, %0.3g) kV/mm\n", + point[0]/mm, point[1]/mm, point[2]/mm, + EMfield[0]/tesla, EMfield[1]/tesla, EMfield[2]/tesla, + EMfield[3]/(kilovolt/mm), EMfield[4]/(kilovolt/mm), EMfield[5]/(kilovolt/mm)); + } +} diff --git a/geant4/LEMuSR/src/G4DecayWithSpin.cc b/geant4/LEMuSR/src/G4DecayWithSpin.cc new file mode 100644 index 0000000..eef7c49 --- /dev/null +++ b/geant4/LEMuSR/src/G4DecayWithSpin.cc @@ -0,0 +1,175 @@ +// +// ******************************************************************** +// * 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" + +G4DecayWithSpin::G4DecayWithSpin(const G4String& processName):G4Decay(processName){} + +G4DecayWithSpin::~G4DecayWithSpin(){} + +G4VParticleChange* G4DecayWithSpin::DecayIt(const G4Track& aTrack, const G4Step& aStep) +{ +// 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) { + // G4cout << "++++ LOCAL VERSION OF G4 DECAY WITH SPIN !!!! ++++"<< G4endl; + 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[3]; + 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); + } + + } + } + +// 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); + + 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(); + + G4Vector3D newSpin = SpinRotation * Spin; + +#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; + +} diff --git a/geant4/LEMuSR/src/MuDecayChannel.cc b/geant4/LEMuSR/src/MuDecayChannel.cc new file mode 100644 index 0000000..10cd200 --- /dev/null +++ b/geant4/LEMuSR/src/MuDecayChannel.cc @@ -0,0 +1,197 @@ +// +// ******************************************************************** +// * 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 " <DumpInfo(); + } +#endif + return products; +} + + + + + + diff --git a/geant4/LEMuSR/src/MuDecayChannelWithSpin.cc b/geant4/LEMuSR/src/MuDecayChannelWithSpin.cc new file mode 100644 index 0000000..10d08d3 --- /dev/null +++ b/geant4/LEMuSR/src/MuDecayChannelWithSpin.cc @@ -0,0 +1,266 @@ +// +// ******************************************************************** +// * 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"<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 " <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; +} diff --git a/geant4/LEMuSR/src/lem4Axial2DElField.cc b/geant4/LEMuSR/src/lem4Axial2DElField.cc new file mode 100644 index 0000000..e8db609 --- /dev/null +++ b/geant4/LEMuSR/src/lem4Axial2DElField.cc @@ -0,0 +1,341 @@ +#include "lem4Axial2DElField.hh" +#include "G4UnitsTable.hh" +//#include "lem4Parameters.hh" + +lem4Axial2DElField::lem4Axial2DElField(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), + ffieldValue(fieldValue) + +//lem4Axial2DElField::lem4Axial2DElField(const G4String filename, double fieldValue, double lenUnit, double fieldNormalisation, double offset) : ffieldValue(fieldValue/1000.0), zOffset(offset), invertR(false), invertZ(false) +{ + // double lenUnit = decimeter; = 100 // Base unit mm + // double fieldUnit = kV; = 0.001 // Base unit MV + // double fieldNorm = 0.00001 = 1E-5 // = 0.001/100 + G4cout << "\n-----------------------------------------------------------" + << "\n Electric field" + << "\n-----------------------------------------------------------" + << G4endl; + // G4cout << " kilovolt = "<< kilovolt << G4endl; // = 1e-3 MV - default G4 unit + ///G4cout << "\n Potential set to "<< fieldValue/1000.0/kilovolt << " kV"<< G4endl; + G4cout << "\n Potential set to "<< fieldValue/kilovolt << " kV"<< G4endl; + G4cout << "\n ---> " "Reading 2D E-field map from " << filename << " ... " << G4endl; + std::ifstream file(filename); // Open the file for reading + + // Start reading file content + char buffer[256]; + + // Read table dimensions first + file >> nr >> nz; // r - radial, z - longitudinal coordinate + + G4cout << " 2D field table dimensions (r-z): [" + << nr << ", " << nz << "] " + << G4endl; + + // Set up storage space for table + rField.resize(nr); + zField.resize(nr); + int ir, iz; + + for (ir = 0; ir < nr; ir++) { + rField[ir].resize(nz); + zField[ir].resize(nz); + } + + // Ignore header information. All lines whose first character + // is '%' are considered to be part of the header. + do { + file.getline(buffer, 256); + } while (buffer[0] != '%'); + + // Read in the data: [r, z, Er, Ez] + double rval, zval, er, ez; + for (ir = 0; ir < nr; ir++) { + for (iz = 0; iz < nz; iz++) { + file >> rval >> zval >> er >> ez; + if (ir == 0 && iz == 0) { + minr = rval * lenUnit; + minz = zval * lenUnit; + } + rField[ir][iz] = er*fieldNormalisation; + zField[ir][iz] = ez*fieldNormalisation; + } + } + G4cout << " Normalisation coeff.: " << fieldNormalisation << G4endl; + file.close(); + + maxr = rval * lenUnit; + maxz = zval * lenUnit; + + G4cout << "\n ---> ... done reading (assumed order: r, z, Er, Ez)." << G4endl; + G4cout << "\n ---> Min values r, z: " + << minr/cm << " " << minz/cm << " cm " + << "\n ---> Max values r, z: " + << maxr/cm << " " << maxz/cm << " cm " << G4endl; + + // Should really check that the limits are not the wrong way around. + if (maxr < minr) {Invert("x");} ///{std::swap(maxr, minr); invertR = true;} + if (maxz < minz) {Invert("z");} ///{std::swap(maxz, minz); invertZ = true;} + if (maxr < minr || maxz < minz){ ///(invertR == true || invertZ == true){ + G4cout << "\n After reordering:" + << "\n ---> Min values r, z: " + << minr/cm << " " << minz/cm << " cm " + << "\n ---> Max values r, z: " + << maxr/cm << " " << maxz/cm << " cm " << G4endl; + } + + dr = maxr - minr; + dz = maxz - minz; + + G4cout << " ---> Range of values: " + << dr/cm << " cm (in r) and " << dz/cm << " cm (in z)." + << "\n-----------------------------------------------------------\n" << G4endl; +} + + +///lem4Axial2DElField::~lem4Axial2DElField() /// Made virtual! +///{ +/// delete [] &rField ; +/// delete [] &zField ; +///} + + +void lem4Axial2DElField::addFieldValue(const G4double point[4], + G4double *field) const +{ + G4double B[3]; // Field value obtained from the field table + ///TS + ///G4cout<<"TTfileTT Global coord. ("<GetName()<0) ? 1.:-1.; + + // Check that the point is within the defined region + if ( r < maxr && z < maxz ) { + // Relative point position within region, normalized to [0,1]. + double rfraction = (r - minr) / dr; + double zfraction = (z - minz) / dz; + + // Need addresses of these to pass to modf below. + // modf uses its second argument as an OUTPUT argument. + double rdindex, zdindex; + + // Position of the point within the cuboid defined by the + // nearest surrounding tabulated points + double rlocal = ( modf(rfraction*(nr-1), &rdindex)); + 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 rindex = static_cast(rdindex); + int zindex = static_cast(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 ((rindex<0)||(rindex>(nr-2))) { + std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="<GetLatestEventNr(); + // G4int evNrKriz=6795; //457; //14250 + // if (evNr==evNrKriz) { + // std::cout<<"evNr="<0) ? 1.:-1.; + //if (evNr==evNrKriz) std::cout<<"r ="<evNrKriz) std::cout<<"bol som tu"<(rdindex); + int zindex = static_cast(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 ((rindex<0)||(rindex>(nr-2))) { + std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="<0) ? Efield_R * (point[0]/r) : 0.; + Efield[1] = (r>0) ? Efield_R * (point[1]/r) : 0.; + Efield[2] = + zField[rindex ][zindex ] * (1-rlocal) * (1-zlocal) + + zField[rindex ][zindex+1] * (1-rlocal) * zlocal + + zField[rindex+1][zindex ] * rlocal * (1-zlocal) + + zField[rindex+1][zindex+1] * rlocal * zlocal ; +// ATTENTION if dealing with MAGNETIC FIELDS reverse x, y and z signs!! + Efield[0] = Efield[0] * ffieldValue; // * z_sign; + Efield[1] = Efield[1] * ffieldValue; // * z_sign; + Efield[2] = Efield[2] * ffieldValue * z_sign; + } + + + + // Set some small field if field is almost zero (to avoid internal problems of Geant). + if (sqrt(Efield[0]*Efield[0]+Efield[1]*Efield[1]+Efield[2]*Efield[2])<1.0E-12*volt) { + // if (evNr>evNrKriz) std::cout<<"bol som tuna"<-1*cm)&&(point[2]<1*cm)) { + // printf ("for point= %f %f %f B= %10.10f %10.10f %10.10f \n", + // point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla); + // } + // if (evNr>evNrKriz) std::cout<<"this is the end"< > rFieldTemp(rField); + std::vector< std::vector< double > > zFieldTemp(zField); + G4bool invertR=false; + G4bool invertZ=false; + + G4cout<<"Check that the lem4Axial2DElField::Invert() function works properly!"<OpenGeometry(); + G4PhysicalVolumeStore::GetInstance()->Clean(); + G4LogicalVolumeStore::GetInstance()->Clean(); + G4SolidStore::GetInstance()->Clean(); + + ///G4cout<< "Parameter file name = "<lem4Error(FATAL,eMessage,false); + } + + + } + else ReportGeometryProblem(line); + + + pFieldMgr->SetDetectorField(lem4Field); + + G4Mag_SpinEqRhs *Mag_SpinEqRhs; + G4MagIntegratorStepper* pStepper; + Mag_SpinEqRhs = new G4Mag_SpinEqRhs(lem4Field); + + pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12); + G4ChordFinder *pChordFinder = new G4ChordFinder(lem4Field,0.01*mm,pStepper); + pFieldMgr->SetChordFinder(pChordFinder ); + pointerToField[nameOfField]=pFieldMgr; + } + + else if(strcmp(tmpString2,"electric")==0){ + + } + } + + + + else if (strcmp(tmpString1,"construct")==0){ + float x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0; + float posx,posy,posz; + char name[100]; + char mothersName[100]; + char material[100]; + // G4CSGSolid* solid=NULL; + G4VSolid* solid=NULL; + G4String solidName="sol_"; + char rotMatrix[100]="norot"; + char sensitiveDet[100]="dead"; + int volumeID=0; + char actualFieldName[100]="nofield"; + + if (strcmp(tmpString2,"tubs")==0){ + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s", + name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + solid = new G4Tubs(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg); + } + else if (strcmp(tmpString2,"box")==0){ + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %s %g %g %g %s %s", + name,&x1,&x2,&x3,material,&posx,&posy,&posz,mothersName,rotMatrix); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + solid = new G4Box(solidName,x1*mm,x2*mm,x3*mm); + } +//Addition of Cones - TS + else if (strcmp(tmpString2,"cones")==0){ + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %g %s %g %*g %*g %*s %*s", + name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,material,&posx); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %*g %g %g %s %s %s %d %s",&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + solid = new G4Cons(solidName,x1*mm,x2*mm,x3*mm,x4*mm,x5*mm,x6*deg,x7*deg); + } + else if (strcmp(tmpString2,"sphere")==0){ + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %s %g %g %g %s %s", + name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg); + } +/*! TS else if (strcmp(tmpString2,"uprofile")==0){ + // Create a U-profile geometry. x1, x2, x3 define the outer dimensions of the U-profile (as a box), + // x4 is the wall thickness of the U-profile. The centre of the U-profile + // is in the centre of the box defined by x1,x2,x3. + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %s %g %g %g %s %s", + name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + + G4double roundingErr=0.01*mm; + G4Box* box1 = new G4Box("Box1",x1*mm,x2*mm,x3*mm); + G4Box* box2 = new G4Box("Box2",(x1-x4)*mm,(x2-x4/2.+roundingErr)*mm,x3*mm); + G4RotationMatrix rot(0,0,0); + G4ThreeVector zTrans(0,x4/2.*mm,0); + G4Transform3D transform(rot,zTrans); + solid = new G4SubtractionSolid(solidName, box1, box2, transform); + } + else if (strcmp(tmpString2,"alcSupportPlate")==0){ + // Create an ALC holder geometry: x1 half-width of the holder (as a box), + // x2 half-height of the spacer + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %s %g %g %g %s %s", + name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + + G4Box* box1 = new G4Box("Box1",x1*mm,4.7*mm,130*mm); + G4Box* box2 = new G4Box("Box2",(x1-0.4)*mm,3*mm,130*mm); + G4RotationMatrix rot(0,0,0); + G4ThreeVector zTrans(0,1.3*mm,0); + G4Transform3D transform(rot,zTrans); + G4SubtractionSolid* solid_1 = new G4SubtractionSolid("Drzak", box1, box2, transform); + if (x2!=0) { + G4Box* box3 = new G4Box("Box3",12.5*mm,4.5*mm,4*mm); + G4ThreeVector zTrans2(0,(-4.7-x2)*mm,0); + G4Transform3D transform2(rot,zTrans2); + solid = new G4UnionSolid("solidName", solid_1, box3, transform2); + } + else { + solid = solid_1; + } + + } TS: These solids are not used by lem4 */ + else if (strcmp(tmpString2,"tubsbox")==0){ + // Create a tube, from which center a box is cut out. x1=box half-width; x2,x3,x4,x5 define the tube. + sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s", + name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName); + sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + G4double roundingErr=0.01*mm; // to avoid some displaying problems of the subtracted volumes + G4Box* solidInnerDetBox = new G4Box("SolidInnerDetBox",x1*mm,x1*mm,x3*mm+roundingErr); + G4Tubs* solidOuterDetTube = new G4Tubs("SolidInnerDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg); + solid = new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetBox); + } + else if (strcmp(tmpString2,"tubsboxsegm")==0){ + // Create a volume that looks like an intersection of tube and box. + char orientation[100]; + sscanf(&line[0],"%*s %*s %*s %s %s %g %g %g %g %g %s %g %g %g %s", + name,orientation,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName); + sscanf(&line[0],"%*s %*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName); + solidName+=name; + G4double roundingErr=0.01*mm; // to avoid some displaying problems of the subtracted volumes + G4Box* solidDetBox = new G4Box("SolidDetBox",x2*mm,x2*mm,x3*mm+roundingErr); + G4Tubs* solidDetTube = new G4Tubs("SolidDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg); + G4RotationMatrix* yRot = new G4RotationMatrix; + G4ThreeVector zTrans; + // if (strcmp(orientation,"D")==0) zTrans=G4ThreeVector(-x1*mm,-3*x1*mm,0); + // else if (strcmp(orientation,"U")==0) zTrans=G4ThreeVector(x1*mm,3*x1*mm,0); + // else if (strcmp(orientation,"R")==0) zTrans=G4ThreeVector(-3*x1*mm,x1*mm,0); + // else if (strcmp(orientation,"L")==0) zTrans=G4ThreeVector(3*x1*mm,-x1*mm,0); + if (strcmp(orientation,"D")==0) zTrans=G4ThreeVector((x1-x2)*mm,(-x1-x2)*mm,0); + else if (strcmp(orientation,"U")==0) zTrans=G4ThreeVector((x2-x1)*mm,(x1+x2)*mm,0); + else if (strcmp(orientation,"R")==0) zTrans=G4ThreeVector((-x1-x2)*mm,(x2-x1)*mm,0); + else if (strcmp(orientation,"L")==0) zTrans=G4ThreeVector((x1+x2)*mm,(x1-x2)*mm,0); + else { + G4cout<<"Unknown orientation of the tubsboxsegm volume!!" + <<" orientation="<SetLogicalVolumeAsSpecialSaveVolume(logicName,volumeID); + lem4RootOutput::GetRootInstance()->SetSpecialSaveVolumeDefined(); + } + } + + // Unless specified as "dead", set the volume sensitive volume: + if (strcmp(sensitiveDet,"dead")) { + if (strcmp(sensitiveDet,"lem4/ScintSD")==0) { + if(!aScintSD) { + G4SDManager* SDman = G4SDManager::GetSDMpointer(); + G4String scintSDname = sensitiveDet; + aScintSD = new lem4ScintSD( scintSDname ); + SDman->AddNewDetector( aScintSD ); + G4cout<<"lem4DetectorConstruction.cc: aScintSD added: "<GetFullPathName()<SetSensitiveDetector( aScintSD ); + } + else { + G4cout<<" lem4DetectorConstruction.cc: unknown sensitive detector \""<SetVolumeIDMapping(logicName,volumeID); + } + } + + else if (strcmp(tmpString1,"visattributes")==0){ + sscanf(&line[0],"%*s %*s %*s %s",tmpString3); + if (strncmp(tmpString2,"log_",4)==0) { + G4LogicalVolume* pLogVol = FindLogicalVolume(tmpString2); + if (pLogVol==NULL) { + G4cout << "ERROR! lem4DetectorConstruction::Construct(): Logical Volume \"" + << tmpString2 <<"\" not found!"<size(); i++) { + G4LogicalVolume* pLogVol=pLogStore->at(i); + G4String materialName=pLogVol->GetMaterial()->GetName(); + if (tmpString2==materialName) { + SetColourOfLogicalVolume(pLogVol,tmpString3); + } + } + } + } + + + // if (strcmp(tmpString3,"red" )==0) {pLogVol->SetVisAttributes(G4Colour(1,0,0));} + // else if (strcmp(tmpString3,"green" )==0) {pLogVol->SetVisAttributes(G4Colour(0,1,0));} + // else if (strcmp(tmpString3,"blue" )==0) {pLogVol->SetVisAttributes(G4Colour(0,0,1));} + // else if (strcmp(tmpString3,"white" )==0) {pLogVol->SetVisAttributes(G4Colour(1,1,1));} + // else if (strcmp(tmpString3,"yellow" )==0) {pLogVol->SetVisAttributes(G4Colour(1,1,0));} + // else if (strcmp(tmpString3,"black" )==0) {pLogVol->SetVisAttributes(G4Colour(0,0,0));} + // else if (strcmp(tmpString3,"gray" )==0) {pLogVol->SetVisAttributes(G4Colour(0.5,0.5,0.5));} + // else if (strcmp(tmpString3,"cyan" )==0) {pLogVol->SetVisAttributes(G4Colour(0,1,1));} + // else if (strcmp(tmpString3,"magenta")==0) {pLogVol->SetVisAttributes(G4Colour(1,0,1));} + // else if (strcmp(tmpString3,"invisible" )==0) {pLogVol->SetVisAttributes(G4VisAttributes::Invisible);} + // else { + // G4cout<<"ERROR: lem4DetectorConstruction.cc unknown G4VisAttributes requested:"<SetVisAttributes(tmpVisAttributes); + } + + // else if (strcmp(tmpString1,"setsensitive")==0){ + // G4SDManager* SDman = G4SDManager::GetSDMpointer(); + // if (strcmp(tmpString2,"lem4/ScintSD")==0) { + // G4String scintSDname = tmpString2; //"lem4/ScintSD"; + // if(!aScintSD) { + // aScintSD = new lem4ScintSD( scintSDname ); + // SDman->AddNewDetector( aScintSD ); + // G4cout<<"lem4DetectorConstruction.cc: aScintSD added: "<GetFullPathName()<SetSensitiveDetector( aScintSD ); + // myRootOutput->SetSensitiveDetectorMapping(tmpString3,volumeID); + // } + // else if (strcmp(tmpString2,"dummy")==0) { + // int volumeID; + // sscanf(&line[0],"%*s %*s %*s %s %d",tmpString3,&volumeID); + // myRootOutput->SetSensitiveDetectorMapping(tmpString3,volumeID); + // } + // else { + // G4cout<<" lem4DetectorConstruction.cc: unknown sensitive detector \""<lem4Error(FATAL,eMessage,false); + } + strcpy(howIsTheFieldDefined,"GLOBAL_FIELD"); + // ensure the global field is initialized + // perhaps should be placed at some separate position ? + (void)F04GlobalField::getObject(); + + char typeOfField[100]="Unset"; + float pp1=0; float pp2=0; float pp3=0; + sscanf(&line[0],"%*s %*s %*s %g %g %g %s",&pp1,&pp2,&pp3,typeOfField); + G4ThreeVector position = G4ThreeVector(pp1,pp2,pp3); + ///sscanf(&line[0],"%*s %*s %*s %s",typeOfField); //typeOfField = uniform, gaussian, or fromfile + float fieldValue = 0.000000001; + float fieldValueFinal = 0; + int fieldNrOfSteps = 0; + // if (strcmp(tmpString2,"magnetic")==0){ + + + if (strcmp(typeOfField,"fromfile")==0) { + char fieldInputFileName[100]; + char fieldTableType[100]; + char fieldType; + char logicalVolumeName[100]; + ///sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %s %s %g %g %d", + /// fieldTableType,fieldInputFileName,logicalVolumeName, + /// &fieldValue,&fieldValueFinal,&fieldNrOfSteps); + sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %s %s %g %g %d", + fieldTableType, fieldInputFileName, logicalVolumeName, + &fieldValue, &fieldValueFinal, &fieldNrOfSteps); + + // Read type of field, either electric E, or magnetic B and set + // accordingly a used-defined DEFAULT unit: E -> kV/mm (0.001), B -> T (0.001) + + fieldType = fieldTableType[2]; + + G4double EMfieldUnit = 1.; + if (fieldType == 'E') {EMfieldUnit = kilovolt/mm;}// = 0.001 -> MV is default in G4 (mm is def. length unit) + else if (fieldType == 'B') {EMfieldUnit = tesla;} // = 0.001 -> kT is default in G4! + ///G4cout<< "Assumed " << fieldType <<"-field unit: "<< EMfieldUnit << G4endl; + + + if ((fieldType != 'E') && (fieldType != 'B')) { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%c\" .",fieldType); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + + + ///sscanf(&line[0],"%*s %*s %*s %*s %s %s %s %g %g %d",fieldTableType,fieldInputFileName,logicalVolumeName, + /// &fieldValue,&fieldValueFinal,&fieldNrOfSteps); + // Find out the logical volume, to which the field will be placed: + G4LogicalVolume* logVol = FindLogicalVolume(logicalVolumeName); + if (logVol==NULL) { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): GLOBAL FIELD: Logical volume \"%s\" not found.", + logicalVolumeName); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + + // Find the coordinates of the volume in its mother + // First construct the name of the physical volume out of the logical volume name: + ///G4String stringLogName=logicalVolumeName; + ///int lastLocation = stringLogName.size () - 1; + ///G4String physicalVolumeName = "phys" + stringLogName.substr(3,lastLocation); + // G4cout<<"physicalVolumeName="<GetVolume(physicalVolumeName); + ///if (physVol==NULL) { + /// sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): GLOBAL FIELD: Physical volume \"%s\" not found.", + /// physicalVolumeName.c_str()); + /// lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + ///} + ///G4ThreeVector position = physVol->GetObjectTranslation(); + ///G4cout<<"position of the volume "<SetElementFieldName(tmpString2); + if (fieldNrOfSteps>0) { + myElementTableField->SetEventNrDependentField(fieldValue*EMfieldUnit,fieldValueFinal*EMfieldUnit,fieldNrOfSteps); + } + } + //* --- Changed by TS to check error in field steps + /* + lem4Axial2DElField* myElementTableField = /// CHECK WHETHER CORRECT!! TS + new lem4Axial2DElField(fieldInputFileName, fieldValue*kilovolt, 10*cm, 0.00001, logVol, position); + myElementTableField->SetElementFieldName(tmpString2); + if (fieldNrOfSteps>0) { + myElementTableField->SetEventNrDependentField(fieldValue*tesla,fieldValueFinal*kilovolt,fieldNrOfSteps); + ///myElementTableField->SetEventNrDependentField(true,fieldValue*tesla,fieldValueFinal*tesla,fieldNrOfSteps); + } + // FieldList* fields = F04GlobalField::getObject()->getFields(); + // if (fields) { + // G4cout<<"Kamil: Detector construction HHHHHH fields->size()="<size()<SetElementFieldName(tmpString2); + if (fieldNrOfSteps>0) { + myElementTableField->SetEventNrDependentField(fieldValue*EMfieldUnit,fieldValueFinal*EMfieldUnit,fieldNrOfSteps); + } + } + } + + else if (strncmp(fieldTableType,"3D",2)==0) { + //G4cout << "Field unit in 3D field map: " << EMfieldUnit << " - Field type: "<< fieldType << G4endl; + //else if (strcmp(fieldTableType,"3D")==0) { + lem4TabulatedElementField3D* myElementTableField = + ///new lem4TabulatedElementField3D(fieldInputFileName, fieldValue*tesla, 1*m, 1., logVol, position); + new lem4TabulatedElementField3D(fieldInputFileName, fieldType, fieldValue*EMfieldUnit, logVol, position); + myElementTableField->SetElementFieldName(tmpString2); + if (fieldNrOfSteps>0) { + myElementTableField->SetEventNrDependentField(fieldValue*EMfieldUnit,fieldValueFinal*EMfieldUnit,fieldNrOfSteps); + } + } + else { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%s\" .",fieldTableType); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + } + + + else if (strcmp(typeOfField,"uniform")==0) { + //// G4cout<<"Uniform field through global F04 call."<lem4Error(FATAL,eMessage,false); + } +*/ + float fieldValue[6]; + char logicalVolumeName[100]; + sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %g %g %g %g %g %g", ////// %g %d", + ///sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %g %g %d", + /// logicalVolumeName, &fieldValue,&fieldValueFinal,&fieldNrOfSteps); + logicalVolumeName, &fieldValue[0],&fieldValue[1],&fieldValue[2],&fieldValue[3],&fieldValue[4],&fieldValue[5]);///,&fieldValueFinal,&fieldNrOfSteps); + ///sscanf(&line[0],"%*s %*s %*s %*s %s %s %s %g %g %d",fieldTableType,fieldInputFileName,logicalVolumeName, + /// &fieldValue,&fieldValueFinal,&fieldNrOfSteps); + // Find out the logical volume, to which the field will be placed: + + //sscanf(&line[0],"%*s %*s %*s %*s %s %g %g %d", logicalVolumeName, &fieldValue, &fieldValueFinal,&fieldNrOfSteps); + // Find out the logical volume, to which the field will be placed: + G4LogicalVolume* logVol = FindLogicalVolume(logicalVolumeName); + if (logVol==NULL) { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): GLOBAL FIELD: Logical volume \"%s\" not found.", logicalVolumeName); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + + + /* ATTENTION!! THIS PIECE OF CODE WAS CAUSING PROBLEMS WITH FIELD NESTING! + // Find the coordinates of the volume in its mother + // First construct the name of the physical volume out of the logical volume name: + G4String stringLogName = logicalVolumeName; + int lastLocation = stringLogName.size () - 1; + G4String physicalVolumeName = "phys" + stringLogName.substr(3,lastLocation); + // G4cout<<"physicalVolumeName="<GetVolume(physicalVolumeName); + if (physVol==NULL) { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): GLOBAL FIELD: Physical volume \"%s\" not found.", + physicalVolumeName.c_str()); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + G4ThreeVector position = physVol->GetObjectTranslation(); + G4cout<<"The position of the volume "<GetName()<SetElementFieldName(tmpString2); ///NOT NEEDED + +/// if (fieldNrOfSteps>0) { +/// myElementField->SetEventNrDependentField(fieldValue*tesla,fieldValueFinal*tesla,fieldNrOfSteps); + ///myElementField->SetEventNrDependentField(true,fieldValue*tesla,fieldValueFinal*tesla,fieldNrOfSteps); +/// } + ////// THE FOLLOWING LINES ARE ONLY FOR CHECKS. TONI + FieldList* fields = F04GlobalField::getObject()->getFields(); + if (fields) { + ////// G4cout<<"Toni: Detector construction HHHHHH fields->size()="<size()<lem4Error(FATAL,eMessage,false); + } + +/// if (strcmp(tmpString2,"uniform")==0){ + /// nParameters=sscanf(&line[0],"%*s %*s %*s %g %g %d",&fieldValue,&fieldValueFinal,&fieldNrOfSteps); + /// SetUniformMagField(fieldValue*tesla); + /// G4cout<<"lem4DetectorConstruction.cc: Uniform Magnetic field set to "<GetFieldManager(); + G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); + if (fieldMgr==NULL) { + ReportProblemInStearingFile(line); + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): G4FieldManager not found: fieldMgr=NULL"); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + if (propagMgr==NULL) { + ReportProblemInStearingFile(line); + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): G4PropagatorInField not found: propagMgr=NULL"); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + else { + char parameterName[100]; + float parameterValue; + sscanf(&line[0],"%*s %*s %*s %s %g",parameterName,¶meterValue); + if (strcmp(parameterName,"SetDeltaIntersection")==0){ fieldMgr->SetDeltaIntersection(parameterValue*mm); } + else if (strcmp(parameterName,"SetDeltaOneStep")==0){ fieldMgr->SetDeltaOneStep(parameterValue*mm); } + else if (strcmp(parameterName,"SetMinimumEpsilonStep")==0){ fieldMgr->SetMinimumEpsilonStep(parameterValue); } + else if (strcmp(parameterName,"SetMaximumEpsilonStep")==0){ fieldMgr->SetMaximumEpsilonStep(parameterValue); } + else if (strcmp(parameterName,"SetLargestAcceptableStep")==0) { propagMgr->SetLargestAcceptableStep(parameterValue*mm); } + else if (strcmp(parameterName,"SetMaxLoopCount")==0) {propagMgr->SetMaxLoopCount(int(parameterValue)); } + else { + G4cout<<"lem4DetectorConstruction.cc: ERROR: Unknown parameterName \"" + <GetFieldManager(); + G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); + if (fieldMgr==NULL) { + G4cout<<"lem4DetectorConstruction: ERROR: Field manager not found!"<AddPointForFieldTesting(G4ThreeVector(p0,p1,p2)); + ///myGlobalField->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); + } + else { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): printFieldValueAtPoint requested, but field not found"); + lem4ErrorMessage::GetInstance()->lem4Error(SERIOUS,eMessage,false); + } + } + } + + else {ReportGeometryProblem(line);} + + } + + //! END OF GUMPLINGER IMPORTED CODE! TS + + + else if (strcmp(tmpString1,"magneticfield")==0){ //else if (strcmp(tmpString1,"magneticfield")==0) - OLD METHOD + //check that the field has not been defined in some other way + if (strcmp(howIsTheFieldDefined,"GLOBAL_FIELD")==0) { + sprintf(eMessage,"lem4DetectorConstruction.cc::Construct(): Field is already defined! howIsTheFieldDefined=\"%s\" .", + howIsTheFieldDefined); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + strcpy(howIsTheFieldDefined,"DIFFERENT_FIELDS"); + float fieldValue=0.000000001; + float fieldValueFinal=0; + int fieldNrOfSteps=0; + float fieldOption=0; + G4int nParameters=0; + if (strcmp(tmpString2,"uniform")==0){ + fieldOption=1; + nParameters=sscanf(&line[0],"%*s %*s %*s %g %g %d",&fieldValue,&fieldValueFinal,&fieldNrOfSteps); + SetUniformMagField(fieldValue*tesla); + G4cout<<"lem4DetectorConstruction.cc: Uniform Magnetic field set to "<SetPositionOffset(offset); } + } + else if (strcmp(tmpString2,"setOffsetOf3DField")==0){ + float x,y,z; + sscanf(&line[0],"%*s %*s %*s %g %g %g", &x, &y, &z); + double offset[4]; + offset[0]=x*mm; offset[1]=y*mm; offset[2]=z*mm; offset[3]=0; + if (lem4TabulatedField3D::GetInstance()==NULL) {G4cout<<"WWWWWWWWWWWWWWWWWWWWWWWWWW"<SetPositionOffset(offset); } + } + + + + else if (strcmp(tmpString2,"setparameter")==0){ // Set the accuracy parameters for the magnetic field + // First check that the magnetic field already exists: + G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); + G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); + if (fieldMgr==NULL) { + G4cout<SetDeltaIntersection(parameterValue*mm); } + else if (strcmp(parameterName,"SetDeltaOneStep")==0){ fieldMgr->SetDeltaOneStep(parameterValue*mm); } + else if (strcmp(parameterName,"SetMinimumEpsilonStep")==0){ fieldMgr->SetMinimumEpsilonStep(parameterValue); } + else if (strcmp(parameterName,"SetMaximumEpsilonStep")==0){ fieldMgr->SetMaximumEpsilonStep(parameterValue); } + else if (strcmp(parameterName,"SetLargestAcceptableStep")==0) { propagMgr->SetLargestAcceptableStep(parameterValue*mm); } + else if (strcmp(parameterName,"SetMaxLoopCount")==0) {propagMgr->SetMaxLoopCount(int(parameterValue)); } + else { + G4cout<<"lem4DetectorConstruction.cc: ERROR: Unknown parameterName \"" + <GetFieldManager(); + G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField(); + if (fieldMgr==NULL) { + G4cout<<"lem4DetectorConstruction: ERROR: Field manager not found!"<GetFieldManager(); + if(!fieldMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field + const G4Field *myfield; + myfield = fieldMgr->GetDetectorField(); + myfield -> 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); + } + } + else {ReportGeometryProblem(line);} + + if (nParameters>0) { // nParameters is >0 only if "uniform", "gaussian" or "fromfile" + // (avoids reseting the stored values for other steering options) + myRootOutput->StoreGeantParameter(0,fieldOption); // store the information about the fieldOption to Root output + myRootOutput->StoreGeantParameter(1,fieldValue); // store the information about the fieldValue to Root output + } + + // If requested, set the parameters for the time dependent magnetic field. + lem4EventAction* pointerToEventAction = lem4EventAction::GetInstance(); + if (nParameters>2) { + pointerToEventAction->SetTimeDependentField(true,fieldValue*tesla,fieldValueFinal*tesla,fieldNrOfSteps); + } + else { + if ((strcmp(tmpString2,"uniform")==0)||(strcmp(tmpString2,"gaussian")==0)||(strcmp(tmpString2,"fromfile")==0)) { + // The following line is needed just to properly fill the fieldValue into the Root output. + pointerToEventAction->SetTimeDependentField(false,fieldValue*tesla,fieldValue*tesla,1); + } + } + } + + //! Set range cut for a given volume, if requested by the macro file - Added by TS + else if (strcmp(tmpString1,"SetUserMinEkine")==0){ + G4LogicalVolume* pLogVol = FindLogicalVolume(tmpString2); + if (pLogVol==NULL) { + G4cout << "ERROR! lem4DetectorConstruction::Construct(): Logical Volume \"" + << tmpString3 <<"\" not found!"< minEnergy value ignored"<SetUserMinEkine(minEnergy*eV); + pLogVol->SetUserLimits(User_lim); + G4cout<<"\nlem4DetectorConstruction.cc: Energy limit in "< maxStep value ignored"<SetSensitiveDetector( aScintSD ); + } + else { + pLogVol->SetUserLimits(new G4UserLimits(maxStep*mm)); + G4cout<<"lem4DetectorConstruction.cc: Range cut in "< SetMyTypeOfProcesses(tmpString2); + } + + else if (strcmp(tmpString1,"storeOnlyEventsWithHits")==0){ + if (strcmp(tmpString2,"false")==0){ lem4Parameters::storeOnlyEventsWithHits = false; } + } + + else if (strcmp(tmpString1,"storeOnlyTheFirstTimeHit")==0){ + if (strcmp(tmpString2,"true")==0){ lem4Parameters::storeOnlyTheFirstTimeHit = true; } + } + + else if (strcmp(tmpString1,"signalSeparationTime")==0){ + float timeSeparation; + sscanf(&line[0],"%*s %*s %g",&timeSeparation); + lem4Parameters::signalSeparationTime = timeSeparation*nanosecond; + } + + //! TS: added to include Muonium formation and processes in the C-foil - by default true + else if (strcmp(tmpString1,"includeMuoniumProcesses")==0){ + if (strcmp(tmpString2,"false")==0){ + lem4Parameters::includeMuoniumProcesses = false; + G4cout<<"\nIgnoring Muonium formation and other Mu processes.\n"<FindOrBuildElement("Fe"); + G4Element* Ni = man->FindOrBuildElement("Ni"); + // Elements required for Brass + G4Element* Cu = man->FindOrBuildElement("Cu"); + G4Element* Zn = man->FindOrBuildElement("Zn"); + // Elements required for MCPglass + G4Element* Pb = man->FindOrBuildElement("Pb"); + G4Element* O = man->FindOrBuildElement("O" ); + G4Element* Si = man->FindOrBuildElement("Si"); + G4Element* K = man->FindOrBuildElement("K" ); + G4Element* Rb = man->FindOrBuildElement("Rb"); + G4Element* Ba = man->FindOrBuildElement("Ba"); + G4Element* As = man->FindOrBuildElement("As"); + G4Element* Cs = man->FindOrBuildElement("Cs"); + G4Element* Na = man->FindOrBuildElement("Na"); + + // Elements required for Macor + G4Element* B = man->FindOrBuildElement("B" ); + G4Element* Al = man->FindOrBuildElement("Al"); + G4Element* Mg = man->FindOrBuildElement("Mg"); + + // compounds required for MCP Macor + G4Material* MgO = new G4Material("MgO", 3.60*g/cm3, ncomponents=2); + MgO->AddElement(Mg, natoms=1); + MgO->AddElement(O, natoms=1); + + G4Material* SiO2 = new G4Material("SiO2", 2.533*g/cm3, ncomponents=2); // quartz + SiO2->AddElement(O, natoms=2); + SiO2->AddElement(Si, natoms=1); + + G4Material* Al2O3 = new G4Material("Al2O3", 3.985*g/cm3, ncomponents=2); // saphire + Al2O3->AddElement (Al, natoms=2); + Al2O3->AddElement (O, natoms=3); + + G4Material* K2O = new G4Material("K2O", 2.350*g/cm3, ncomponents=2); + K2O->AddElement(O, natoms=1); + K2O->AddElement(K, natoms=2); + + G4Material* B2O3 = new G4Material("B2O3", 2.550*g/cm3, ncomponents=2); + B2O3->AddElement (B, natoms=2); + B2O3->AddElement (O, natoms=3); + + + // Define materials from Elements: Case 2 - mixture by fractional mass + + G4Material* brass = // Common Brass + new G4Material("Brass", density=8.40*g/cm3, ncomponents=2); + brass->AddElement(Zn, fractionmass=30*perCent); + brass->AddElement(Cu, fractionmass=70*perCent); + + G4Material* steel = // Stainless Steel + new G4Material("Steel", density=7.93*g/cm3, ncomponents=3); + steel->AddElement(Ni, fractionmass=0.11); + steel->AddElement(Cr, fractionmass=0.18); + steel->AddElement(Fe, fractionmass=0.71); + + G4Material* macor= // Macor (used in the MCP detector) + new G4Material("Macor", density=2.52*g/cm3, ncomponents=5); + macor->AddMaterial(SiO2, fractionmass=0.470); // quartz + macor->AddMaterial(MgO, fractionmass=0.180); + macor->AddMaterial(Al2O3,fractionmass=0.170); // saphire + macor->AddMaterial(K2O, fractionmass=0.105); + macor->AddMaterial(B2O3, fractionmass=0.075); + + G4Material* mcpglass = // Glass of the Multi Channel Plate + new G4Material("MCPglass", density=2.0*g/cm3, ncomponents=9); + mcpglass->AddElement(Pb, fractionmass= 0.480); + mcpglass->AddElement(O, fractionmass= 0.258); + mcpglass->AddElement(Si, fractionmass= 0.182); + mcpglass->AddElement(K, fractionmass= 0.042); + mcpglass->AddElement(Rb, fractionmass= 0.018); + mcpglass->AddElement(Ba, fractionmass= 0.013); + mcpglass->AddElement(As, fractionmass= 0.004); + mcpglass->AddElement(Cs, fractionmass= 0.002); + mcpglass->AddElement(Na, fractionmass= 0.001); + +} + +//G4cout<<"Materials successfully defined."<DefineWorldVolume(Construct()); +} + +void lem4DetectorConstruction::SetUniformMagField(G4double fieldValue) +{ // If we want a uniform field along z + + G4FieldManager *pFieldMgr; + lem4MagneticField* lem4UniformField= new lem4MagneticField(); + lem4UniformField->SetFieldValue(fieldValue); + + pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); + G4Mag_SpinEqRhs *Mag_SpinEqRhs; + G4MagIntegratorStepper* pStepper; + Mag_SpinEqRhs = new G4Mag_SpinEqRhs(lem4UniformField); + pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12); + G4cout<< "DeltaStep "<GetDeltaOneStep()/mm <<"mm" <SetChordFinder( pChordFinder ); + pFieldMgr->SetDetectorField(lem4UniformField); + + // store the pointer to the lem4MagneticField in the event action, which is needed for + // the time dependent magnetic field. + lem4EventAction* pointerToEventAction = lem4EventAction::GetInstance(); + pointerToEventAction->SetPointerToMusrUniformField(lem4UniformField); + + // G4PropagatorInField* fieldPropagator + // = G4TransportationManager::GetTransportationManager() + // ->GetPropagatorInField(); + // fieldPropagator->SetMinimumEpsilonStep(1.e-5*mm); + // fieldPropagator->SetMaximumEpsilonStep(1.e-2*mm); + + //cks For better visualisation of the tracks in low density materials use SetLargestAcceptableStep() + // G4TransportationManager* tmanager = G4TransportationManager::GetTransportationManager(); + // // tmanager->GetPropagatorInField()->SetLargestAcceptableStep(1*mm); + // // tmanager->GetPropagatorInField()->SetLargestAcceptableStep(2000*mm); + // if (largestAcceptableStep>0.000000001*mm) { + // tmanager->GetPropagatorInField()->SetLargestAcceptableStep(largestAcceptableStep); + // } + // G4double LargestAcceptableStep=tmanager->GetPropagatorInField()->GetLargestAcceptableStep(); + // G4cout<< "GetLargestAcceptableStep()="<SetPointerToTabulatedField3D(myMusrField); + } + else if (strcmp(fieldTableType,"2D")==0) { + lem4TabulatedField2D* myMusrField = new lem4TabulatedField2D(inputFileName, fieldValue, 1*cm, 0.00001); + lem4Field = myMusrField; + pointerToEventAction->SetPointerToTabulatedField2D(myMusrField); + } + else { + char eMessage[200]; + sprintf(eMessage,"lem4DetectorConstruction.cc::SetMagField(): unknown type of the tabulated field \"%s\" .",fieldTableType); + lem4ErrorMessage::GetInstance()->lem4Error(FATAL,eMessage,false); + } + + // // just to test the values of the field + // for (int i=0; i<2000; i++) { + // double point[4]={0,0,i-1000,0}; + // // double point[4]={-50,50,-50,0}; + // double Bfi[3]={0,0,0}; + // lem4Field->GetFieldValue(point,Bfi); + // printf ("for point= %f %f %f B= %10.10f %10.10f %10.10f \n", + // point[0],point[1],point[2],Bfi[0]/tesla,Bfi[1]/tesla,Bfi[2]/tesla); + // } + pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); + + G4Mag_SpinEqRhs *Mag_SpinEqRhs; + G4MagIntegratorStepper* pStepper; + + Mag_SpinEqRhs = new G4Mag_SpinEqRhs(lem4Field); + pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12); + G4cout<< "DeltaStep "<GetDeltaOneStep()/mm <<"mm" <SetChordFinder( pChordFinder ); + pFieldMgr->SetDetectorField(lem4Field); +} + + + +// void lem4DetectorConstruction::SetGaussianMagField(G4double fieldValue, G4double fieldSigma, G4double dummy) { +void lem4DetectorConstruction::SetGaussianMagField(G4ThreeVector inputParam) { + //------------ Field along z, which changes intensity as a (gaussian) function of R + G4FieldManager *pFieldMgr; + // G4double fieldRMS=100*mm; + G4double fieldValue=inputParam.x()*tesla; + G4double fieldSigma=inputParam.y()*mm; + G4MagneticField* lem4Field= new lem4GaussianField(fieldValue,fieldSigma); + pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); + + G4Mag_SpinEqRhs *Mag_SpinEqRhs; + G4MagIntegratorStepper* pStepper; + + Mag_SpinEqRhs = new G4Mag_SpinEqRhs(lem4Field); + pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12); + G4cout<< "DeltaStep "<GetDeltaOneStep()/mm <<"mm" <SetChordFinder( pChordFinder ); + pFieldMgr->SetDetectorField(lem4Field); +} + + + +void lem4DetectorConstruction::ReportGeometryProblem(char myString[501]) { + G4cout<<"\nE R R O R in lem4DetectorConstruction.cc: " + <<"Unknown keyword requested in \""<< parameterFileName <<"\" :"<FindOrBuildMaterial(materialName); + //Material=G4Material::GetMaterial(materialName); +//! } + if (Material==NULL) { + G4cout<<"\nE R R O R in lem4DetectorConstruction.cc::CharToMaterial " + <<"Unknown material requested:"<entries(); i++) { + for (unsigned int i=0; isize(); i++) { + G4LogicalVolume* pLogVol=pLogStore->at(i); + G4String iLogName=pLogVol->GetName(); + if (iLogName==LogicalVolumeName) { + return pLogVol; + } + } + } + return NULL; +} + +void lem4DetectorConstruction::SetColourOfLogicalVolume(G4LogicalVolume* pLogVol,char* colour) { + if (pLogVol!=NULL) { + if (strcmp(colour,"red" )==0) {pLogVol->SetVisAttributes(G4Colour(1,0,0));} + else if (strcmp(colour,"green" )==0) {pLogVol->SetVisAttributes(G4Colour(0,1,0));} + else if (strcmp(colour,"blue" )==0) {pLogVol->SetVisAttributes(G4Colour(0,0,1));} + else if (strcmp(colour,"white" )==0) {pLogVol->SetVisAttributes(G4Colour(1,1,1));} + else if (strcmp(colour,"yellow" )==0) {pLogVol->SetVisAttributes(G4Colour(1,1,0));} + else if (strcmp(colour,"black" )==0) {pLogVol->SetVisAttributes(G4Colour(0,0,0));} + else if (strcmp(colour,"gray" )==0) {pLogVol->SetVisAttributes(G4Colour(0.5,0.5,0.5));} + else if (strcmp(colour,"cyan" )==0) {pLogVol->SetVisAttributes(G4Colour(0,1,1));} + else if (strcmp(colour,"magenta")==0) {pLogVol->SetVisAttributes(G4Colour(1,0,1));} + else if (strcmp(colour,"invisible" )==0) {pLogVol->SetVisAttributes(G4VisAttributes::Invisible);} + + else if (strcmp(colour,"blue_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.80,0.83,1));} + else if (strcmp(colour,"lightblue")==0) {pLogVol->SetVisAttributes(G4Colour(0,0.5,1));} + else if (strcmp(colour,"darkblue")==0) {pLogVol->SetVisAttributes(G4Colour(0,0.25,0.5));} + else if (strcmp(colour,"fblue_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.85,.88,0.92));} + else if (strcmp(colour,"oxsteel")==0) {pLogVol->SetVisAttributes(G4Colour(0.9,0.8,0.75));} + else if (strcmp(colour,"darkred")==0) {pLogVol->SetVisAttributes(G4Colour(0.5,0,0));} + else if (strcmp(colour,"MCP_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.5,0.2,.7));} + else if (strcmp(colour,"MACOR_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.9,0.9,.1));} + else if (strcmp(colour,"SCINT_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.5,0.5,.75));} + else if (strcmp(colour,"dSCINT_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.3,0.3,0.3));} + else if (strcmp(colour,"VTBB_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.9,0.9,.9));} + else if (strcmp(colour,"Grid_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.87,0.72,0.53));} //burlywood + else if (strcmp(colour,"RA_style")==0) {pLogVol->SetVisAttributes(G4Colour(0.8549,0.6471,0.1255));} //goldenrod + + else { + G4cout<<"ERROR: lem4DetectorConstruction::SetColourOfLogicalVolume: unknown colour requested: "< //cks -----------------------------||------------------------------- +#include "lem4EventAction.hh" // cks needed for the initialisation of the random nr. generator by event nr. +//#include +//#include +#include "globals.hh" +//#include "lem4Parameters.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + + +lem4DetectorMessenger::lem4DetectorMessenger(lem4DetectorConstruction* myDet) + :myDetector(myDet) +{ + lem4Dir = new G4UIdirectory("/lem4/"); + lem4Dir->SetGuidance("UI commands specific to this example."); + + CommandCmd = new G4UIcmdWithAString("/lem4/command",this); + CommandCmd->SetGuidance("This command will be used for the detector construction (and ignored by the default messenger."); + // CommandCmd->SetParameterName("choice",false); + CommandCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + runDir = new G4UIdirectory("/lem4/run/"); + runDir->SetGuidance("lem4 run control"); + + RunIDSetCmd = new G4UIcmdWithAnInteger("/lem4/run/runID",this); + RunIDSetCmd->SetGuidance("Set the run number"); + RunIDSetCmd->SetParameterName("something",false); + RunIDSetCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + RandomOptionCmd = new G4UIcmdWithAnInteger("/lem4/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->SetParameterName("randomOpt",false); + RandomOptionCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + HowOftenToPrintEventCmd = new G4UIcmdWithAnInteger("/lem4/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); + + detDir = new G4UIdirectory("/lem4/det/"); + detDir->SetGuidance("detector control."); + + // WhichProcessesCmd = new G4UIcmdWithAString("/lem4/det/processes",this); + // WhichProcessesCmd ->SetGuidance("Select Standard, Low Energy or Penelope processes"); + // WhichProcessesCmd ->SetParameterName("mes_processes",false); + // WhichProcessesCmd ->AvailableForStates(G4State_PreInit,G4State_Idle); + + UpdateCmd = new G4UIcmdWithoutParameter("/lem4/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); + + // FieldCmd = new G4UIcmdWithADoubleAndUnit("/lem4/det/setField",this); + // FieldCmd->SetGuidance("Define magnetic field."); + // FieldCmd->SetGuidance("Magnetic field will be in Z direction."); + // FieldCmd->SetParameterName("Bz",false); + // FieldCmd->SetUnitCategory("Magnetic flux density"); + // FieldCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + UFieldCmd = new G4UIcmdWithADoubleAndUnit("/lem4/det/setUniformField",this); + UFieldCmd->SetGuidance("Define uniform magnetic field."); + UFieldCmd->SetGuidance("Magnetic field will be in Z direction."); + UFieldCmd->SetParameterName("Bz",false); + UFieldCmd->SetUnitCategory("Magnetic flux density"); + UFieldCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + + // GFieldCmd = new G4UIcmdWithADoubleAndUnit("/lem4/det/setGaussianField",this); + // GFieldCmd = new G4UIcmdWith3VectorAndUnit("/lem4/det/setGaussianField",this); + GFieldCmd = new G4UIcmdWith3Vector("/lem4/det/setGaussianField",this); + GFieldCmd->SetGuidance("Define gaussian magnetic field: intensity in T, sigma in mm, and a dummy variable."); + GFieldCmd->SetGuidance("Magnetic field will be in Z direction."); + GFieldCmd->SetParameterName("Bz","sigmaBz","dummy",true,false); + // GFieldCmd->SetUnitCategory("Magnetic flux density"); + GFieldCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + +/// UEFieldCmd = new G4UIcmdWithADoubleAndUnit("/lem4/det/setUnifEField",this); +/// UEFieldCmd->SetGuidance("Define uniform electric field."); +/// UEFieldCmd->SetGuidance("Electric field will be in Z direction."); +/// UEFieldCmd->SetParameterName("Ez",false); +/// UEFieldCmd->SetUnitCategory("Electric field"); +/// UEFieldCmd->AvailableForStates(G4State_PreInit,G4State_Idle); + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4DetectorMessenger::~lem4DetectorMessenger() +{ + // delete WhichProcessesCmd; + delete UpdateCmd; + delete detDir; + delete lem4Dir; + delete CommandCmd; + // delete FieldCmd; + delete UFieldCmd; + delete GFieldCmd; +/// delete UEFieldCmd; + delete RunIDSetCmd; + delete RandomOptionCmd; + delete HowOftenToPrintEventCmd; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4DetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue) { + + // if( command == WhichProcessesCmd ) + // { lem4Parameters* myParameters = lem4Parameters::GetInstance(); + // lem4Parameters -> SetMyProcesses(newValue); } + // lem4Parameters* myParameters = lem4Parameters::GetInstance(); + // myParameters -> strcnpy(lem4Parameters::myProcesses,newValue,99); + // lem4Parameters::myProcesses[99]='\0'; + + if( command == UpdateCmd ) + { myDetector->UpdateGeometry(); } + + // if( command == FieldCmd ) + // { myDetector->SetMagField("fieldprofiles/input.dat",FieldCmd->GetNewDoubleValue(newValue));} + + if( command == UFieldCmd ) + { myDetector->SetUniformMagField(UFieldCmd->GetNewDoubleValue(newValue));} + + if( command == GFieldCmd ) + { myDetector->SetGaussianMagField(GFieldCmd->GetNew3VectorValue(newValue));} + +/// if( command == UEFieldCmd ) +/// { myDetector->SetUniformElField(UEFieldCmd->GetNewDoubleValue(newValue));} + + 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; + lem4EventAction::setRandomNrSeedAccordingEventNr=1; + // lem4EventAction::setMyEventNr(70); + } + else if (RandomOption == 3) { + G4cout << "*******************************************" << G4endl; + G4cout << "*** Random Seed set from external file ***" << G4endl; + G4cout << "*******************************************" << G4endl; + lem4EventAction::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); + } + // vector * seedVector = new vector; + vector * seedVector = lem4EventAction::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 + + // lem4EventAction::RandomNrInitialisers->push_back(num); + } + indata.close(); + G4cout << "End-of-file reached.." << seedVector->size()<GetNewIntValue(newValue); + lem4EventAction::nHowOftenToPrintEvent=n; + } + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/LEMuSR/src/lem4ErrorMessage.cc b/geant4/LEMuSR/src/lem4ErrorMessage.cc new file mode 100644 index 0000000..847e0fe --- /dev/null +++ b/geant4/LEMuSR/src/lem4ErrorMessage.cc @@ -0,0 +1,62 @@ +#include "lem4ErrorMessage.hh" + +lem4ErrorMessage::lem4ErrorMessage():nErrors(1) +{ + pointerToErrors=this; + severityWord[INFO]="INFO"; + severityWord[WARNING]="WARNING"; + severityWord[SERIOUS]="SERIOUS"; + severityWord[FATAL]="FATAL"; +} + +lem4ErrorMessage::~lem4ErrorMessage() {} + +lem4ErrorMessage* lem4ErrorMessage::pointerToErrors=NULL; +lem4ErrorMessage* lem4ErrorMessage::GetInstance() { + return pointerToErrors; +} + +void lem4ErrorMessage::lem4Error(SEVERITY severity, G4String message, G4bool silent) { + std::map::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<<"!!!"<1) { + G4cout<<"!!!"<::iterator it; + G4cout<<"------ ERROR SUMMARY: -------------------------------------------------------"< * lem4EventAction::RandomNrInitialisers=NULL; + +//long lem4EventAction::myEventNr=0; + +lem4EventAction::lem4EventAction() { + pointer=this; + fieldValueStart=0; + pointerToSeedVector = new vector; + 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); + // + + lem4SteppingAction::GetInstance()->DoAtTheBeginningOfEvent(); + + long thisEventNr = (long) (evt->GetEventID()); + + if (F04GlobalField::Exists()) { + F04GlobalField::getObject() -> CheckWhetherAnyNominalFieldValueNeedsToBeChanged(thisEventNr); + } + + + ///?long thisEventNr = (long) (evt->GetEventID()); + 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 "<SetFieldValue(actualFieldValue); + G4cout<<"Event "<SetFieldValue(actualFieldValue); + G4cout<<"Event "<SetFieldValue(actualFieldValue); + + // if (lem4DetectorMessenger::setRandomNrSeedAccordingEventNr) { + if (setRandomNrSeedFromFile) { + // G4cout<<"RandomNrInitialisers.size()="<size()<size()) { + G4cout <<"lem4EventAction.cc: seed will be set to="<< pointerToSeedVector->at(thisEventNr)<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="<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() <<"\t 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; iGetTrajectoryContainer()))[i]); + trj->DrawTrajectory(1000); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +vector * lem4EventAction::pointerToSeedVector=NULL; +vector * 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<<"---&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&---"<SetDetectorField(this); + GetGlobalFieldManager()->CreateChordFinder(this); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4MagneticField::lem4MagneticField(G4ThreeVector fieldVector) + : G4UniformMagField(fieldVector) +{ + GetGlobalFieldManager()->SetDetectorField(this); + GetGlobalFieldManager()->CreateChordFinder(this); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Set the value of the Global Field to fieldValue along Z +// +void lem4MagneticField::SetFieldValue(G4double fieldValue) +{ + G4UniformMagField::SetFieldValue(G4ThreeVector(0,0,fieldValue)); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Set the value of the Global Field +// +void lem4MagneticField::SetFieldValue(G4ThreeVector fieldVector) +{ + // Find the Field Manager for the global field + G4FieldManager* fieldMgr= GetGlobalFieldManager(); + + if(fieldVector!=G4ThreeVector(0.,0.,0.)) + { + G4UniformMagField::SetFieldValue(fieldVector); + fieldMgr->SetDetectorField(this); + } else { + // If the new field's value is Zero, then it is best to + // insure that it is not used for propagation. + G4MagneticField* magField = NULL; + fieldMgr->SetDetectorField(magField); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4MagneticField::~lem4MagneticField() +{ + // GetGlobalFieldManager()->SetDetectorField(0); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "G4TransportationManager.hh" + +G4FieldManager* lem4MagneticField::GetGlobalFieldManager() +{ + return G4TransportationManager::GetTransportationManager()->GetFieldManager(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/LEMuSR/src/lem4MuFormation.cc b/geant4/LEMuSR/src/lem4MuFormation.cc new file mode 100644 index 0000000..289de9e --- /dev/null +++ b/geant4/LEMuSR/src/lem4MuFormation.cc @@ -0,0 +1,106 @@ +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +// Muonium Formation according to yield.cc function (through GetYields method). +// Id : lem4MuFormation.cc, v 1.4 +// Author: Taofiq PARAISO, T. Shiroka +// Date : 2007-12 +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +#include "lem4MuFormation.hh" + +using namespace std; + +lem4MuFormation::lem4MuFormation(const G4String& name, G4ProcessType aType) + : G4VDiscreteProcess(name, aType){} + +lem4MuFormation::~lem4MuFormation(){} + +G4VParticleChange* lem4MuFormation::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 lem4MuFormation::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 lem4MuFormation::GetMeanFreePath(const G4Track&, + G4double, + G4ForceCondition* condition) +{ + *condition = Forced; + return DBL_MAX; +} + + +void lem4MuFormation::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(rndFindParticle(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 lem4MuFormation::PrepareSecondary(const G4Track& track) +{ + if(p_name=="mu+") + { + aSecondary = new G4Track(DP,track.GetGlobalTime(),track.GetPosition()); + } +} diff --git a/geant4/LEMuSR/src/lem4MuScatter.cc b/geant4/LEMuSR/src/lem4MuScatter.cc new file mode 100644 index 0000000..9f5e028 --- /dev/null +++ b/geant4/LEMuSR/src/lem4MuScatter.cc @@ -0,0 +1,81 @@ +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ +// Muonium "Scattering" +// Id : lem4MuScatter.cc, v 1.4 +// Author: Taofiq PARAISO, T. Shiroka +// Date : 2007-12 +// Notes : Simplified model for Mu scattering. Spin effects have been excluded. +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ + +#include "lem4MuScatter.hh" + +using namespace std; + +lem4MuScatter::lem4MuScatter(const G4String& name, + G4ProcessType aType) + : G4VDiscreteProcess(name, aType){} + +lem4MuScatter:: ~lem4MuScatter(){} + +/*! - 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* lem4MuScatter::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 lem4MuScatter::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 lem4MuScatter::GetMeanFreePath(const G4Track&, + G4double, + G4ForceCondition* condition) +{ + *condition = Forced; + return DBL_MAX; +} + + +void lem4MuScatter::PrepareSecondary(const G4Track& track) +{ + aSecondary = new G4Track(DP,track.GetDynamicParticle()->GetPreAssignedDecayProperTime(),track.GetPosition()); +} diff --git a/geant4/LEMuSR/src/lem4Muonium.cc b/geant4/LEMuSR/src/lem4Muonium.cc new file mode 100644 index 0000000..33b9897 --- /dev/null +++ b/geant4/LEMuSR/src/lem4Muonium.cc @@ -0,0 +1,107 @@ +// +// ******************************************************************** +// * 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: lem4Muonium.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 "lem4Muonium.hh" +#include "G4ParticleTable.hh" + +#include "MuDecayChannel.hh" +#include "G4DecayTable.hh" + +// ###################################################################### +// ### MUONIUM ### +// ###################################################################### +lem4Muonium* lem4Muonium::theInstance = 0; + +lem4Muonium* lem4Muonium::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(anInstance); + return theInstance; +} + +lem4Muonium* lem4Muonium::MuoniumDefinition() +{ + return Definition(); +} + +lem4Muonium* lem4Muonium::Muonium() +{ + return Definition(); +} + diff --git a/geant4/LEMuSR/src/lem4Parameters.cc b/geant4/LEMuSR/src/lem4Parameters.cc new file mode 100644 index 0000000..ada9b2e --- /dev/null +++ b/geant4/LEMuSR/src/lem4Parameters.cc @@ -0,0 +1,71 @@ +#include "lem4Parameters.hh" + +///lem4Parameters::lem4Parameters() +/// :typeOfProcesses("standard") +///{ +/// pointerToParameters=this; +///} + +lem4Parameters::lem4Parameters(G4String steeringFileName) + :typeOfProcesses("coulombAndMultiple") //lowenergy") +{ + pointerToParameters=this; + + // 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 lem4DetectorConstruction.cc class). + + FILE *fSteeringFile=fopen(steeringFileName.c_str(),"r"); + if (fSteeringFile==NULL) { + G4cout<<"lem4Parameters::lem4Parameters: steeringFileName=\""< Insert(new G4MuonDecayChannelWithSpin("mu+",1.00)); + G4MuonPlus::MuonPlusDefinition() -> SetDecayTable(MuonPlusDecayTable); + + //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)); + lem4Muonium::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable); + //MuoniumDecayTable ->DumpInfo(); // Info on muonium decay channels + +} + +// Mesons (only some light mesons) + +void lem4PhysicsList::ConstructMesons() +{ + G4PionPlus::PionPlusDefinition(); + G4PionMinus::PionMinusDefinition(); + G4PionZero::PionZeroDefinition(); +} + +// Baryons + +void lem4PhysicsList::ConstructBaryons() +{ + G4Proton::ProtonDefinition(); + G4AntiProton::AntiProtonDefinition(); + G4Neutron::NeutronDefinition(); + G4AntiNeutron::AntiNeutronDefinition(); +} + + + + +// Physical processes - general + +void lem4PhysicsList::ConstructProcess() +{ + AddTransportation(); + ConstructEM(); + ConstructGeneral(); +} + +// Declaration of the different processes +#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 "lem4AtRestSpinRotation.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" + +// For Muonium formation in the Carbon foil +#include "lem4MuFormation.hh" // includes the yield function Y = Y(E). + +// For a simple Muonium "scattering" when Mu hits solid materials +#include "lem4MuScatter.hh" + + +// Construction methods for the different processes + +void lem4PhysicsList::ConstructEM() +{ + G4String myTypeOfProcesses = lem4Parameters::GetInstance()->GetMyTypeOfProcesses(); + G4cout<<" uuuuuuuuu myTypeOfProcesses="<reset(); + while( (*theParticleIterator)() ){ + G4ParticleDefinition* particle = theParticleIterator->value(); + G4ProcessManager* pmanager = particle->GetProcessManager(); + G4String particleName = particle->GetParticleName(); + + if (particleName == "gamma") { // gamma + if (myTypeOfProcesses=="lowenergy") { + pmanager->AddDiscreteProcess(new G4LowEnergyPhotoElectric); + pmanager->AddDiscreteProcess(new G4LowEnergyCompton); + pmanager->AddDiscreteProcess(new G4LowEnergyGammaConversion); + pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh); + } + else if (myTypeOfProcesses=="penelope") { + pmanager->AddDiscreteProcess(new G4PenelopePhotoElectric); + pmanager->AddDiscreteProcess(new G4PenelopeCompton); + pmanager->AddDiscreteProcess(new G4PenelopeGammaConversion); + pmanager->AddDiscreteProcess(new G4PenelopeRayleigh); + } + else { // Standard G4 (default) + pmanager->AddDiscreteProcess(new G4PhotoElectricEffect); + pmanager->AddDiscreteProcess(new G4ComptonScattering); + pmanager->AddDiscreteProcess(new G4GammaConversion); + } + + } + else if (particleName == "e-") { //electron + if (myTypeOfProcesses=="lowenergy") { + pmanager->AddProcess(new G4MultipleScattering, -1, 1,1); + pmanager->AddProcess(new G4LowEnergyIonisation, -1, 2,2); + pmanager->AddProcess(new G4LowEnergyBremsstrahlung,-1,-1,3); + } + else if (myTypeOfProcesses=="penelope") { + pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1); + pmanager->AddProcess(new G4PenelopeIonisation, -1, 2, 2); + pmanager->AddProcess(new G4PenelopeBremsstrahlung, -1,-1, 3); + } + else if (myTypeOfProcesses=="coulomb") { + pmanager->AddDiscreteProcess(new G4CoulombScattering); + pmanager->AddProcess(new G4eIonisation, -1, 2,2); + pmanager->AddProcess(new G4eBremsstrahlung, -1, 3,3); + } + else { // Standard G4 (default) + pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); + pmanager->AddProcess(new G4eIonisation, -1, 2,2); + pmanager->AddProcess(new G4eBremsstrahlung, -1, 3,3); + } + } + + else if (particleName == "e+") { + + if (myTypeOfProcesses=="penelope") { + pmanager->AddProcess(new G4MultipleScattering, -1, 1, 1); + pmanager->AddProcess(new G4PenelopeIonisation, -1, 2, 2); + pmanager->AddProcess(new G4PenelopeBremsstrahlung, -1,-1, 3); + pmanager->AddProcess(new G4PenelopeAnnihilation, 0,-1, 4); + } + else if (myTypeOfProcesses=="coulomb") { + pmanager->AddDiscreteProcess(new G4CoulombScattering); + pmanager->AddProcess(new G4eIonisation, -1, 2,2); + pmanager->AddProcess(new G4eBremsstrahlung, -1, 3,3); + pmanager->AddProcess(new G4eplusAnnihilation, 0,-1,4); + } + else { // Standard G4 (default) + pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); + pmanager->AddProcess(new G4eIonisation, -1, 2,2); + pmanager->AddProcess(new G4eBremsstrahlung, -1, 3,3); + pmanager->AddProcess(new G4eplusAnnihilation, 0,-1,4); + } + } + + else if ((particleName=="mu+")||(particleName=="mu-")) { //muon + + if (myTypeOfProcesses=="coulomb") { + pmanager->AddDiscreteProcess(new G4CoulombScattering); // TS - Attention use Meyer!! + } + else if (myTypeOfProcesses=="coulombAndMultiple") { + std::cout <<"musrPhysicsList: Switching between Coulomb and Multiple scattering for mu+"<AddProcess(new G4MultipleScattering,-1, 1,1); + pmanager->AddDiscreteProcess(new G4CoulombScattering); + } + else { + pmanager->AddProcess(new G4MultipleScattering,-1, 1,1); // TS - Attention use Meyer!! + } + + if (inclMuProcesses){ // This accounts ONLY for Mu formation! + G4VProcess* aMuoniumFormation = new lem4MuFormation(); + pmanager->AddProcess(aMuoniumFormation); + pmanager->SetProcessOrdering(aMuoniumFormation,idxPostStep,2); + //pmanager->AddProcess(new lem4MuFormation, -1,-1,2); // TS Muonium formation process - shorthand + } + + pmanager->AddProcess(new G4MuIonisation, -1, 2,2); + pmanager->AddProcess(new G4MuBremsstrahlung, -1, 3,3); + pmanager->AddProcess(new G4MuPairProduction, -1, 4,4); + //cks Change to Geant 4 default muon decay with spin (instead of using the code of Taofiq). + // pmanager->AddProcess(new lem4AtRestSpinRotation, 1, -1, -1); + // lem4Decay* theDecayProcess = new lem4Decay(); + // pmanager->AddProcess(theDecayProcess); + // pmanager ->SetProcessOrderingToLast(theDecayProcess, idxAtRest); + // pmanager ->SetProcessOrdering(theDecayProcess, idxPostStep); + // + // // The spin precession (rotation) at rest is already included in + // // the G4DecayWithSpin(); + // G4AtRestSpinPrecession* theSpinPrecession = new G4AtRestSpinPrecession(); + // pmanager->AddRestProcess(theSpinPrecession); + + + + if (particle->GetParticleName()=="mu+"){ + G4DecayWithSpin* theDecayProcess = new G4DecayWithSpin(); + //lem4DecayWithSpin* theDecayProcess = new lem4DecayWithSpin(); + pmanager->AddProcess(theDecayProcess); + pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest); + pmanager->SetProcessOrdering(theDecayProcess, idxPostStep); + } + } + + + else if ((particle->GetParticleName()=="Mu")&&(inclMuProcesses)) { // other Mu processes + G4cout<<"\nAccounting for Muonium formation and other Mu processes.\n"<AddProcess(new G4MultipleScattering, -1, 1, 1); + //pmanager->AddProcess(new lem4MuScatter, -1, 1, 1); <--- ERROR ?!? + //pmanager->AddProcess(new lem4MuScatter, -1, -1, 1); //<--- ERROR ?!? + + // Muonium "scattering" + G4VProcess* aMuScatt = new lem4MuScatter(); + pmanager->AddProcess(aMuScatt); + pmanager->SetProcessOrdering(aMuScatt, idxPostStep, 1); + + // Muonium decay process + G4Decay* theDecayProcess = new G4Decay(); + //lem4DecayWithSpin* theDecayProcess = new lem4DecayWithSpin(); + pmanager->AddProcess(theDecayProcess); + pmanager->SetProcessOrderingToLast(theDecayProcess, idxAtRest); + pmanager->SetProcessOrdering(theDecayProcess, idxPostStep); + } + + + + /* + else if (particle->GetParticleName()=="Mu") + { + //! + Half of the muonium has an isotropic decay. + // + G4DecayTable* MuoniumDecayTable = new G4DecayTable(); + MuoniumDecayTable -> Insert(new G4MuonDecayChannelWithSpin("Mu",0.5)); // ?? with spin? LEMuSRMuonDecayChannel("Mu",0.5)); + MuoniumDecayTable -> Insert(new G4MuonDecayChannel("Mu",0.5)); + MuoniumDecayTable ->DumpInfo(); + MuoniumParticle::MuoniumDefinition() -> SetDecayTable(MuoniumDecayTable); + + G4ProcessManager* pmanager = particle->GetProcessManager(); + if (theLEMuSRDecayProcess->IsApplicable(*particle)) + { + pmanager ->AddProcess(theLEMuSRDecayProcess); + pmanager ->AddProcess(aStepLimiter); + + // set ordering + // pmanager ->SetProcessOrdering(aStepLimiter, idxPostStep, 1); // not needed for the muonium + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxPostStep); + pmanager ->SetProcessOrderingToLast(theLEMuSRDecayProcess, idxAtRest); + } + + pmanager->DumpInfo(); + + + } */ + + //csk + //} + + else if ((!particle->IsShortLived()) && + (particle->GetPDGCharge() != 0.0) && + (particle->GetParticleName() != "chargedgeantino")) { + // All stable, charged particles except geantino + pmanager->AddProcess(new G4MultipleScattering, -1,1,1); + if (myTypeOfProcesses=="lowenergy") { + pmanager->AddProcess(new G4hLowEnergyIonisation, -1,2,2); + } + else { + pmanager->AddProcess(new G4hIonisation, -1,2,2); + } + ///pmanager->AddProcess(new G4UserSpecialCuts, -1,-1,3); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + + + +#include "G4Decay.hh" +void lem4PhysicsList::ConstructGeneral() +{ +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +#include "G4Region.hh" +#include "G4RegionStore.hh" +#include "G4ProductionCuts.hh" + +void lem4PhysicsList::SetCuts() +{ + //G4VUserPhysicsList::SetCutsWithDefault method sets + //the default cut value for all particle types + // + SetCutsWithDefault(); + + //cks Set some cuts specific to a given region: + // G4Region* region; + // G4String regName; + // G4ProductionCuts* cuts; + // + // regName="Target_region"; + // region = G4RegionStore::GetInstance()->GetRegion(regName); + // if (region==NULL) { + // G4cout << "lem4PhysicsList.cc: ERROR! Region "<SetProductionCut(0.005*mm,G4ProductionCuts::GetIndex("mu+")); + // cuts->SetProductionCut(0.005*mm,G4ProductionCuts::GetIndex("e-")); + // region->SetProductionCuts(cuts); + //csk + + if (verboseLevel>0) DumpCutValuesTable(); + DumpCutValuesTable(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + diff --git a/geant4/LEMuSR/src/lem4PrimaryGeneratorAction.cc b/geant4/LEMuSR/src/lem4PrimaryGeneratorAction.cc new file mode 100644 index 0000000..3baa990 --- /dev/null +++ b/geant4/LEMuSR/src/lem4PrimaryGeneratorAction.cc @@ -0,0 +1,318 @@ +#include "lem4PrimaryGeneratorAction.hh" +#include "lem4DetectorConstruction.hh" +#include "lem4PrimaryGeneratorMessenger.hh" +#include "G4Event.hh" +#include "G4ParticleGun.hh" +#include "G4ParticleTable.hh" +#include "G4ParticleDefinition.hh" +#include "Randomize.hh" +#include "G4ios.hh" +#include "G4ParticleGun.hh" +#include "G4UnitsTable.hh" +#include "globals.hh" +#include "G4Gamma.hh" +#include "G4ThreeVector.hh" +#include "G4RunManager.hh" +#include "time.h" +#include +#include "lem4RootOutput.hh" //cks for storing some info in the Root output file +#include "lem4ErrorMessage.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +lem4PrimaryGeneratorAction::lem4PrimaryGeneratorAction( + lem4DetectorConstruction* lem4DC) + :lem4Detector(lem4DC), x0(0), y0(0), z0(-124*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), xPolarisIni(1.), yPolarisIni(0.), zPolarisIni(0.), + muonDecayTimeMin(-1), muonDecayTimeMax(-1), muonMeanLife(2197.03*ns), + takeMuonsFromTurtleFile(false) + //, firstCall(true) +{ + G4int n_particle = 1; + //cksdel particleGun = new lem4ParticleGun(n_particle); + particleGun = new G4ParticleGun(n_particle); + + //create a messenger for this class + gunMessenger = new lem4PrimaryGeneratorMessenger(this); + + // default particle kinematic + G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable(); + G4ParticleDefinition* particle= particleTable->FindParticle("mu+"); + particleGun->SetParticleDefinition(particle); + mu_mass = particle->GetPDGMass(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +lem4PrimaryGeneratorAction::~lem4PrimaryGeneratorAction() +{ + delete particleGun; + delete gunMessenger; + if (takeMuonsFromTurtleFile) {fclose(fTurtleFile);} +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +void lem4PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) +{ + // This function is called at the begining of event. + + // // First of all check the setting of some variables. Do it only once (for the first event). + // if (firstCall) { + // firstCall=false; + // char message[200]; + // cout<<"******************** called just once **********************"< + // 5*sqrt(xSigma*xSigma+ySigma*ySigma) ) { + // sprintf(message,"lem4PrimaryGeneratorAction::GeneratePrimaries: Too strict restriction on the generated radius."); + // lem4ErrorMessage::GetInstance()->lem4Error(FATAL,"message",true); + // } + // // Check that the restriction on the z-coordinate of the generated postion is reasonable: + // if ( ((-zMaxAllowed+z0) > 5*fabs(zSigma)) || ((zMinAllowed-z0) > 5*fabs(zSigma)) ) { + // sprintf(message,"lem4PrimaryGeneratorAction::GeneratePrimaries: Too strict restriction on the generated z-coordinate."); + // lem4ErrorMessage::GetInstance()->lem4Error(FATAL,"message",true); + // } + // } + + G4double x, y, z; + G4double p, ke_offset; + G4double xangle, yangle; + + ke_offset = 3.73*keV; // Kin. energy offset due to a 3.73 kV acceleration at the Carbon foil + + if (takeMuonsFromTurtleFile) { + char line[501]; + G4int checkNrOfCounts=0; + do { + float xTmp, yTmp, xAngleTmp, yAngleTmp, pTmp; + if (feof(fTurtleFile)) {rewind(fTurtleFile);G4cout<<"End of Turtle file."<1000) { + G4cout<<"lem4PrimaryGeneratorAction::GeneratePrimaries: Too strict requirements on the r position!"<(rMaxAllowed*rMaxAllowed) ); + z=z0; + // G4cout<<"x,y,z=("<(rMaxAllowed*rMaxAllowed))||(z>zMaxAllowed)||(z0) {p = p0;} + // if (E0>0 && p0==0) {p0 = std::sqrt(E0*E0 + 2*mu_mass*E0);} + // else if (E0>0 && p0>0) { + // G4cout<<"Define either kinetic energy or momentum, but not both!"<0) {p = G4RandGauss::shoot(p0,pSigma);} + else {p=p0;} + checkNrOfCounts++; + if (checkNrOfCounts>1000) { + G4cout<<"lem4PrimaryGeneratorAction::GeneratePrimaries: Too strict requirements on the momentum!"<pMaxAllowed)||(p0) { 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)/(1*mm);} //xSigma; } // Changed to absolute units of pitch + + 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)/(1*mm);} //ySigma; } // Changed to absolute units of pitch + + } // 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) { + G4cout<<"lem4PrimaryGeneratorAction.cc: Unpolarised muon beam not yet implemented!"< S T O P"<SetParticlePosition(G4ThreeVector(x,y,z)); + //particleGun->SetParticleMomentum(G4ThreeVector(px,py,pz)); + G4double particleEnergy = std::sqrt(p*p+mu_mass*mu_mass)-mu_mass; + particleGun->SetParticleEnergy(particleEnergy + ke_offset); + particleGun->SetParticleMomentumDirection(G4ThreeVector(px,py,pz)); + particleGun->SetParticlePolarization(G4ThreeVector(xpolaris,ypolaris,zpolaris)); + particleGun->GeneratePrimaryVertex(anEvent); + + // G4cout<<"lem4PrimaryGeneratorAction: Parameters:"<0.) { + // G4cout<<"muonDecayTimeMin="<muonDecayTimeMax)); + // + // This algorithm should be much faster then the previous one: + 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="<SetProperTime(decaytime); + } + + // Save variables into ROOT output file: + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + myRootOutput->SetInitialMuonParameters(x,y,z,px,py,pz,xpolaris,ypolaris,zpolaris); + +} + +/////////////////////////////////////////////////////////////////////// +// +// + + +void lem4PrimaryGeneratorAction::SetInitialMuonPolariz(G4ThreeVector vIniPol) +{ + G4double magnitude=vIniPol.mag(); + if(magnitude<0.00000001) { + G4cout<< "Unpolarised initial muons"<StoreGeantParameter(2,muonDecayTimeMin/microsecond); + myRootOutput->StoreGeantParameter(3,muonDecayTimeMax/microsecond); + myRootOutput->StoreGeantParameter(4,muonMeanLife/microsecond); +} + +void lem4PrimaryGeneratorAction::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;} +} diff --git a/geant4/LEMuSR/src/lem4PrimaryGeneratorMessenger.cc b/geant4/LEMuSR/src/lem4PrimaryGeneratorMessenger.cc new file mode 100644 index 0000000..674f098 --- /dev/null +++ b/geant4/LEMuSR/src/lem4PrimaryGeneratorMessenger.cc @@ -0,0 +1,136 @@ +#include "lem4PrimaryGeneratorMessenger.hh" +#include "lem4PrimaryGeneratorAction.hh" +#include "G4UIcmdWithAString.hh" +#include "G4UIcmdWithADoubleAndUnit.hh" +#include "G4UIcmdWithADouble.hh" +#include "G4UIcmdWithAnInteger.hh" +#include "G4UIcmdWith3Vector.hh" +#include "G4UIcmdWith3VectorAndUnit.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +lem4PrimaryGeneratorMessenger::lem4PrimaryGeneratorMessenger(lem4PrimaryGeneratorAction* lem4Gun) +:lem4Action(lem4Gun) +{ + 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"); + + setEnergyCmd = new G4UIcmdWithADoubleAndUnit("/gun/kenergy",this); + setEnergyCmd->SetGuidance(" Set kinetic energy of the generated muons (with unit)"); + setEnergyCmd->SetParameterName("mes_E0",true); + setEnergyCmd->SetDefaultUnit("keV"); +//setEnergyCmd->SetUnitCandidates("eV keV MeV GeV"); + + 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","mes_zangleSigma",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); + + 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); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +lem4PrimaryGeneratorMessenger::~lem4PrimaryGeneratorMessenger() +{ + delete setvertexCmd; + delete setvertexSigmaCmd; + delete setvertexBoundaryCmd; + delete setEnergyCmd; + delete setMomentumCmd; + delete setMomentumSmearingCmd; + delete setMomentumBoundaryCmd; + delete setTiltAngleCmd; + delete setSigmaTiltAngleCmd; + delete setPitchCmd; + delete setMuonPolarizCmd; + delete setMuonDecayTimeCmd; + delete setTurtleCmd; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... + +void lem4PrimaryGeneratorMessenger::SetNewValue(G4UIcommand *command, G4String newValue) +{ + if( command == setvertexCmd) + { lem4Action->Setvertex(setvertexCmd->GetNew3VectorValue(newValue));} + else if( command == setvertexSigmaCmd) + { lem4Action->SetvertexSigma(setvertexSigmaCmd->GetNew3VectorValue(newValue));} + else if( command == setvertexBoundaryCmd) + { lem4Action->SetvertexBoundary(setvertexBoundaryCmd->GetNew3VectorValue(newValue));} + else if( command == setEnergyCmd) + { lem4Action->SetEnergy(setEnergyCmd->GetNewDoubleValue(newValue));} + else if( command == setMomentumCmd) + { lem4Action->SetMomentum(setMomentumCmd->GetNewDoubleValue(newValue));} + else if( command == setMomentumSmearingCmd) + { lem4Action->SetMomentumSmearing(setMomentumSmearingCmd->GetNewDoubleValue(newValue));} + else if( command == setMomentumBoundaryCmd) + { lem4Action->SetMomentumBoundary(setMomentumBoundaryCmd->GetNew3VectorValue(newValue));} + else if( command == setTiltAngleCmd) + { lem4Action->SetTilt(setTiltAngleCmd->GetNew3VectorValue(newValue));} + else if( command == setSigmaTiltAngleCmd) + { lem4Action->SetSigmaTilt(setSigmaTiltAngleCmd->GetNew3VectorValue(newValue));} + else if( command == setPitchCmd) + { lem4Action->SetPitch(setPitchCmd->GetNewDoubleValue(newValue));} + else if( command == setMuonPolarizCmd) + { lem4Action->SetInitialMuonPolariz(setMuonPolarizCmd->GetNew3VectorValue(newValue));} + else if( command == setMuonDecayTimeCmd) + { lem4Action->SetMuonDecayTimeLimits(setMuonDecayTimeCmd->GetNew3VectorValue(newValue)); } + else if ( command == setTurtleCmd) + { lem4Action->SetTurtleInput(newValue); } +} diff --git a/geant4/LEMuSR/src/lem4RootOutput.cc b/geant4/LEMuSR/src/lem4RootOutput.cc new file mode 100644 index 0000000..45976fe --- /dev/null +++ b/geant4/LEMuSR/src/lem4RootOutput.cc @@ -0,0 +1,353 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "lem4RootOutput.hh" +#include "G4RunManager.hh" +#include "G4Run.hh" +//#include "globals.hh" +#include "lem4ErrorMessage.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4RootOutput::lem4RootOutput() { + pointerToRoot=this; + boolIsAnySpecialSaveVolumeDefined=false; + nFieldNomVal=0; + + // SensDetectorMapping["log_ScintBU"]=1; + // SensDetectorMapping["log_ScintBD"]=2; + // SensDetectorMapping["log_ScintBL"]=3; + // SensDetectorMapping["log_ScintBR"]=4; + // SensDetectorMapping["log_ScintFU"]=5; + // SensDetectorMapping["log_ScintFD"]=6; + // SensDetectorMapping["log_ScintFL"]=7; + // SensDetectorMapping["log_ScintFR"]=8; + + // SensDetectorMapping["logicScintBU"]=1; + // SensDetectorMapping["logicScintBD"]=2; + // SensDetectorMapping["logicScintBL"]=3; + // SensDetectorMapping["logicScintBR"]=4; + // SensDetectorMapping["logicScintFU"]=5; + // SensDetectorMapping["logicScintFD"]=6; + // SensDetectorMapping["logicScintFL"]=7; + // SensDetectorMapping["logicScintFR"]=8; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4RootOutput::~lem4RootOutput() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4RootOutput* lem4RootOutput::pointerToRoot=0; +lem4RootOutput* lem4RootOutput::GetRootInstance() { + return pointerToRoot; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +G4bool lem4RootOutput::store_runID = true; +G4bool lem4RootOutput::store_eventID = true; +G4bool lem4RootOutput::store_BFieldAtDecay = true; +G4bool lem4RootOutput::store_muIniPosX = true; +G4bool lem4RootOutput::store_muIniPosY = true; +G4bool lem4RootOutput::store_muIniPosZ = true; +G4bool lem4RootOutput::store_muIniMomX = true; +G4bool lem4RootOutput::store_muIniMomY = true; +G4bool lem4RootOutput::store_muIniMomZ = true; +G4bool lem4RootOutput::store_muIniPolX = true; +G4bool lem4RootOutput::store_muIniPolY = true; +G4bool lem4RootOutput::store_muIniPolZ = true; +G4bool lem4RootOutput::store_muDecayDetID= true; +G4bool lem4RootOutput::store_muDecayPosX = true; +G4bool lem4RootOutput::store_muDecayPosY = true; +G4bool lem4RootOutput::store_muDecayPosZ = true; +G4bool lem4RootOutput::store_muDecayTime = true; +G4bool lem4RootOutput::store_muDecayPolX = true; +G4bool lem4RootOutput::store_muDecayPolY = true; +G4bool lem4RootOutput::store_muDecayPolZ = true; +G4bool lem4RootOutput::store_muTargetTime = true; +G4bool lem4RootOutput::store_muTargetPolX = true; +G4bool lem4RootOutput::store_muTargetPolY = true; +G4bool lem4RootOutput::store_muTargetPolZ = true; +G4bool lem4RootOutput::store_posIniMomX = true; +G4bool lem4RootOutput::store_posIniMomY = true; +G4bool lem4RootOutput::store_posIniMomZ = true; +G4bool lem4RootOutput::store_globalTime = true; +G4bool lem4RootOutput::store_fieldValue = true; +G4bool lem4RootOutput::store_det_ID = true; +G4bool lem4RootOutput::store_det_edep = true; +G4bool lem4RootOutput::store_det_edep_el = true; +G4bool lem4RootOutput::store_det_edep_pos = true; +G4bool lem4RootOutput::store_det_edep_gam = true; +G4bool lem4RootOutput::store_det_edep_mup = true; +G4bool lem4RootOutput::store_det_nsteps = true; +G4bool lem4RootOutput::store_det_length = true; +G4bool lem4RootOutput::store_det_start = true; +G4bool lem4RootOutput::store_det_end = true; +G4bool lem4RootOutput::store_det_x = true; +G4bool lem4RootOutput::store_det_y = true; +G4bool lem4RootOutput::store_det_z = true; +G4bool lem4RootOutput::store_fieldNomVal = true; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4RootOutput::BeginOfRunAction() { + G4cout << "lem4RootOutput::BeginOfRunAction() Defining Root tree and branches ..."<GetCurrentRun()->GetRunID(); + char RootOutputFileName[200]; + sprintf(RootOutputFileName, "data/lem4_%i.root", tmpRunNr); + // G4cout<<"RootOutputFileName="<Branch("runID",&runID_t,"runID/I");} + if (store_eventID) {rootTree->Branch("eventID",&eventID_t,"eventID/I");} + 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_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");} + + // rootTree->Branch("nSubTracks",&nSubTracks,"nSubTracks/I"); + // rootTree->Branch("particleID",&particleID_t,"particleID[nSubTracks]/I"); + // rootTree->Branch("trackID",&trackID_t,"trackID[nSubTracks]/I"); + // rootTree->Branch("logicalVolumeID",&logicalVolumeID_t,"logicalVolumeID[nSubTracks]/I"); + // rootTree->Branch("edep",&edep_t,"edep[nSubTracks]/D"); + // rootTree->Branch("x",&x_t,"x[nSubTracks]/D"); + // rootTree->Branch("y",&y_t,"y[nSubTracks]/D"); + // rootTree->Branch("z",&z_t,"z[nSubTracks]/D"); + // rootTree->Branch("post_x",&post_x_t,"post_x[nSubTracks]/D"); + // rootTree->Branch("post_y",&post_y_t,"post_y[nSubTracks]/D"); + // rootTree->Branch("post_z",&post_z_t,"post_z[nSubTracks]/D"); + // rootTree->Branch("kineticEnergy",&kineticEnergy_t,"kineticEnergy[nSubTracks]/D"); + // rootTree->Branch("stepLength",&stepLength_t,"stepLength[nSubTracks]/D"); + // rootTree->Branch("nrOfSubsteps",&nrOfSubsteps_t,"nrOfSubsteps[nSubTracks]/I"); + + 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 (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 << "lem4RootOutput::BeginOfRunAction() Root tree and branches were defined.\n"<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; iClose(); + G4cout<<"lem4RootOutput::EndOfRunAction() - Root tree written out.\n"<Fill(atan2(ttt.y(),ttt.x())); + // htest6->Fill(atan2(sqrt(ttt.x()*ttt.x()+ttt.y()*ttt.y()),ttt.z())); + htest5->Fill(atan2(posIniMomy_t,posIniMomx_t)); + htest6->Fill(atan2(sqrt(posIniMomx_t*posIniMomx_t+posIniMomy_t*posIniMomy_t),posIniMomz_t)); + // G4double PolModulus=sqrt(muDecayPolX_t*muDecayPolX_t+muDecayPolY_t*muDecayPolY_t+muDecayPolZ_t*muDecayPolZ_t); + // G4double PolModulus0=sqrt(muTargetPolX_t*muTargetPolX_t+muTargetPolY_t*muTargetPolY_t+muTargetPolZ_t*muTargetPolZ_t); + // if ((muPolx_t>1.0001)||((muPolx_t<-1.0001)&&(muPolx_t>-999.))||(muPolx0_t>1.0001)||((muPolx0_t<-1.0001)&&(muPolx0_t>-999.))) { + // if ((PolModulus>1.001)||((PolModulus0>1.001)&&(muTargetTime_t>-999))) { + // G4cout <<"lem4RootOutput.cc: Strange polarisation in run="<lem4Error(FATAL,message,false); + } + else{ + SensDetectorMapping[logivol]=volumeID; + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4int lem4RootOutput::ConvertVolumeToID(std::string logivol) { + G4int volumeID = SensDetectorMapping[logivol]; + if (volumeID==0) { + char message[100]; + sprintf(message,"lem4RootOutput::ConvertVolumeToID: No ID number assigned to sensitive volume %s .",logivol.c_str()); + lem4ErrorMessage::GetInstance()->lem4Error(SERIOUS,"message",true); + } + return volumeID; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4RootOutput::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,"lem4RootOutput.cc::SetSaveDetectorInfo(): more \"save\" hits then allowed: save_nMax=%i",save_nMax); + lem4ErrorMessage::GetInstance()->lem4Error(SERIOUS,message,true); + } + else { + save_detID[save_n]=ID; + save_particleID[save_n]=particleID; + save_ke[save_n]=ke/keV; + 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 lem4RootOutput::SetFieldNomVal(G4int i, G4double value) { + if (ilem4Error(SERIOUS,message,true); + } +} + diff --git a/geant4/LEMuSR/src/lem4RunAction.cc b/geant4/LEMuSR/src/lem4RunAction.cc new file mode 100644 index 0000000..5da8072 --- /dev/null +++ b/geant4/LEMuSR/src/lem4RunAction.cc @@ -0,0 +1,85 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +// Make G4Timer appear first! +#include "G4Timer.hh" +#include "lem4RunAction.hh" +#include "lem4EventAction.hh" +#include "G4Run.hh" +#include "lem4ErrorMessage.hh" +#include "F04GlobalField.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4RunAction::lem4RunAction() { + timer = new G4Timer; +} + + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4RunAction::~lem4RunAction() { + delete timer; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4RunAction::BeginOfRunAction(const G4Run* aRun) { + timer->Start(); + G4int run_id= aRun->GetRunID(); + // if (run_id%100 == 0) { + G4cout << "### Run " << run_id << G4endl; + // } + // Set nr. of events that will be processed in this run to lem4EventAction class, as it will + // be needed for the time dependent magnetic field; + lem4EventAction* pointerToEventAction = lem4EventAction::GetInstance(); + pointerToEventAction->SetNumberOfEventsInThisRun(aRun->GetNumberOfEventToBeProcessed()); + // + lem4RootOutput::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: ----------------\n"<begin(); i!=fields->end(); ++i) { + (*i)->construct(); + // Get the nominal field value for the given field and store it in the Root output + G4double nomFieldValue = (*i)->GetNominalFieldValue(); + lem4RootOutput::GetRootInstance()->SetFieldNomVal(jjj,nomFieldValue); + jjj++; + } + G4cout<<"-----------------------------------------------------------------"<PrintFieldAtRequestedPoints(); + } + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4RunAction::EndOfRunAction(const G4Run* aRun) { + lem4RootOutput::GetRootInstance()->EndOfRunAction(); + lem4ErrorMessage::GetInstance()->PrintErrorSummary(); + timer->Stop(); + G4cout << "lem4RunAction::EndOfRunAction:"< +#include +#include + +G4Allocator lem4ScintHitAllocator; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ScintHit::lem4ScintHit() { +} + +G4int lem4ScintHit::ScintMultihit=0; +G4int lem4ScintHit::runIDoldScint=-1; +G4int lem4ScintHit::eventIDoldScint=-1; +G4int lem4ScintHit::NIS=0; +G4int lem4ScintHit::ScintChamberNbold=-1; +G4int lem4ScintHit::verboseLevel=0; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ScintHit::~lem4ScintHit() { + //save the Tree header. The file will be automatically closed + //when going out of the function scope + // rootTree->Write(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ScintHit::lem4ScintHit(const lem4ScintHit& right) + : G4VHit() +{ + particleName = right.particleName; + trackID = right.trackID; + edep = right.edep; + pre_pos = right.pre_pos; + pol = right.pol; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +const lem4ScintHit& lem4ScintHit::operator=(const lem4ScintHit& 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 lem4ScintHit::operator==(const lem4ScintHit& right) const +{ + return (this==&right) ? 1 : 0; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4ScintHit::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 lem4ScintHit::Print() +{ + if (verboseLevel>2) G4cout<<"VERBOSE 3: Kamil: lem4ScintHit::Print()"<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<<" "< // needed for the sort() function +#include "G4VProcess.hh" // needed for the degugging message of the process name +#include "G4RunManager.hh" +#include "G4Run.hh" +#include "lem4Parameters.hh" +#include "lem4ErrorMessage.hh" +#include + +//bool myREMOVEfunction (int i,int j) { return (i::iterator i1, std::map::iterator m2) { +// return ( (*i1).first()<(*i2).second() ); +//} +// +//bool timeOrdering2 (std::pair p1, std::pair p2) { +// return ( p1.first()1) G4cout<<"VERBOSE 2: lem4ScintSD::Initialize\n"; + // Positron_momentum_already_stored=0; + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + myRootOutput->ClearAllRootVariables(); + + scintCollection = new lem4ScintHitsCollection + (SensitiveDetectorName,collectionName[0]); + static G4int HCID = -1; + if(HCID<0) { + HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); + if (verboseLevel>1) G4cout<<"VERBOSE 2: lem4ScintSD::HCID was <0\n, now HCID="<AddHitsCollection( HCID, scintCollection ); + myStoreOnlyEventsWithHits = lem4Parameters::storeOnlyEventsWithHits; + mySignalSeparationTime = lem4Parameters::signalSeparationTime; + myStoreOnlyTheFirstTimeHit= lem4Parameters::storeOnlyTheFirstTimeHit; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4bool lem4ScintSD::ProcessHits(G4Step* aStep,G4TouchableHistory*) +{ + if (verboseLevel>1) G4cout<<"VERBOSE 2: lem4ScintSD::ProcessHits\n"; + G4double edep = aStep->GetTotalEnergyDeposit(); + if(edep==0.) { + return false; + } + + G4Track* aTrack = aStep->GetTrack(); + G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName(); + // G4cout <<"actualVolume=="<SetTrackStatus(fSuspend); + // } + + // 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; + } + } + + + //cks Ignore hits from the particles that had been created at the time when the muon + // was stopping in the target - not a correct thing to do, however no better way + // to suppres wrong timing information of the hit found yet. + // lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + // G4double muonArrivedInTarget = myRootOutput->GetTimeInTarget(); + // if (muonArrivedInTarget>0.) { + // G4double timeOfThisTrack = aTrack->GetGlobalTime(); + // if (fabs(muonArrivedInTarget-timeOfThisTrack)<(2*ns)) { + // return false; + // } + // } + // else { + // // G4cout<<"Kamil: in the hit before muon has arrived into the target"<GetGlobalTime()) < (0.03*microsecond) ) return false; + lem4ScintHit* newHit = new lem4ScintHit(); + // G4RunManager* fRunManager = G4RunManager::GetRunManager(); + // newHit->SetEventID(fRunManager->GetCurrentEvent()->GetEventID()); + // newHit->SetRunID(fRunManager->GetCurrentRun()->GetRunID()); + newHit->SetParticleName (aTrack->GetDefinition()->GetParticleName()); + newHit->SetParticleID (aTrack->GetDefinition()->GetPDGEncoding()); + // if (aTrack->GetDefinition()->GetPDGEncoding()==-13) { + // if ((aTrack->GetGlobalTime())<(0.015*microsecond)) { + // G4cout<<"PROCESS: detected particle="<GetDefinition()->GetParticleName() + // <<" TIME="<GetGlobalTime()<GetPostStepPoint()->GetProcessDefinedStep(); + // // G4cout<<"GetProcessName()="<GetProcessName()<GetCreatorProcess(); + // G4cout<<"Process that limitted this step: "; interestingProcess->DumpInfo(); + // if (creatorProcess!=NULL) {G4cout<<"Process that created this particle: "; creatorProcess->DumpInfo();} + // } + newHit->SetTrackID (aTrack->GetTrackID()); + newHit->SetEdep (edep); + newHit->SetPrePos (aStep->GetPreStepPoint()->GetPosition()); + newHit->SetPostPos (aStep->GetPostStepPoint()->GetPosition()); + newHit->SetPol (aTrack->GetPolarization()); + newHit->SetLogVolName(aTrack->GetVolume()->GetLogicalVolume()->GetName()); + newHit->SetGlobTime (aTrack->GetGlobalTime()); + + /// Printout to check proper assignment of Global Time! + //G4cout<<"StepTime="<GetDeltaTime()<<"\t"; + //G4cout<<"GlobalTime="<GetGlobalTime()<<" LocalTime="<GetLocalTime()<<" ProperTime="<GetProperTime(); + //G4cout<<"\t trackID="<GetTrackID() <<"\t actualVolume="<IsFirstStepInVolume() only available in Geant version >= 4.8.2 ! + // newHit->SetFirstStepInVolumeFlag(aStep->IsFirstStepInVolume()); + // newHit->SetLastStepInVolumeFlag(aStep->IsLastStepInVolume()); + newHit->SetKineticEnergy(aTrack->GetVertexKineticEnergy()); + // newHit->SetKineticEnergy(aTrack->GetKineticEnergy()+edep); + newHit->SetStepLength (aStep->GetStepLength()); + // G4double BFieldAtOrigin[6] = {0.,0.,0.,0.,0.,0.}; + // G4double Origin[4] ={0.,0.,0.,0.}; + // G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); + // if (fMgr!=NULL) { + // if(!fMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field + // fMgr->GetDetectorField()->GetFieldValue(Origin,BFieldAtOrigin); + // } + // else{ + // } + // // G4cout<<"Kamil: pole="<SetBField(BFieldAtOrigin); + scintCollection->insert( newHit ); + // newHit->Print(); + newHit->Draw(); + return true; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4ScintSD::EndOfEvent(G4HCofThisEvent*) { + if (verboseLevel>1) { + G4cout<<"VERBOSE 2: lem4ScintSD::EndOfEvent"<entries(); + G4cout << "\n-------->Hits Collection: in this event they are " << NbHits + << " hits in the scint chambers: " << G4endl; + // for (G4int i=0;iPrint(); + } + + // Positron_momentum_already_stored=0; + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + + G4RunManager* fRunManager = G4RunManager::GetRunManager(); + myRootOutput->SetRunID(fRunManager->GetCurrentRun()->GetRunID()); + myRootOutput->SetEventID(fRunManager->GetCurrentEvent()->GetEventID()); + + G4int NbHits = scintCollection->entries(); + if ((NbHits<=0)&&(myStoreOnlyEventsWithHits)) { + myRootOutput->ClearAllRootVariables(); + return; + } + + // Sort out hits and fill them into root + if (NbHits>0) { + lem4ScintHit* firstHit = (*scintCollection)[0]; + myRootOutput->SetGlobTime(firstHit->GetGlobalTime()); + + const G4int det_IDmax = lem4RootOutput::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]; + // for (G4int ii=0; ii myvector (myints, myints+8); // 32 71 12 45 26 80 53 33 + // std::vector::iterator it; + // sort (myvector.begin(), myvector.end(), myREMOVEfunction); + // // print out content: + // G4cout << "myvector contains:"; + // for (it=myvector.begin(); it!=myvector.end(); ++it) + // G4cout << " " << *it; + // G4cout << G4endl; + + + // 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::map myHitTimeMapping; + std::map::iterator it; + for (G4int i=0; iGetGlobalTime(); + // G4cout<<"Hit nr "<first; + // G4cout << " Value:" << it->second << "\n"; + G4int ii = it->second; // ii is the index of the hits, which is sorted according to time + lem4ScintHit* aHit = (*scintCollection)[ii]; + G4String aHitVolumeName = aHit->GetLogVolName(); + G4int aHitVolumeID = myRootOutput->ConvertVolumeToID(aHitVolumeName); + G4double aHitTime = aHit->GetGlobalTime(); + + // Loop over all already defined signals and check whether the hit falls into any of them + G4bool signalAssigned=false; + for (G4int j=0; jGetEdep(); + 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,"lem4ScintSD.cc::EndOfEvent(): untreated particle \"%s\" deposited energy.",aParticleName.c_str()); + lem4ErrorMessage::GetInstance()->lem4Error(WARNING,message,true); + } + 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,"lem4ScintSD.cc::EndOfEvent(): number of signals exceeds maximal allowed value."); + lem4ErrorMessage::GetInstance()->lem4Error(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,"lem4ScintSD.cc::EndOfEvent(): UNTREATED PARTICLE \"%s\" deposited energy.",aParticleName.c_str()); + lem4ErrorMessage::GetInstance()->lem4Error(WARNING,message,true); + } + G4ThreeVector prePos = aHit->GetPrePos(); + det_x[nSignals]=prePos.x(); + det_y[nSignals]=prePos.y(); + det_z[nSignals]=prePos.z(); + 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 mySignalMapping; + std::map::iterator itt; + for (G4int i=0; i(-det_edep[i],i) ); + } + + // Write out the signals (sorted according to energy) to the lem4RootOutput class: + G4int jj=-1; + for (itt=mySignalMapping.begin(); itt!=mySignalMapping.end(); itt++) { + jj++; + G4int ii = itt->second; + myRootOutput->SetDetectorInfo(jj,det_ID[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]); + } + + } //end "if (NbHits>0)" + + myRootOutput->FillEvent(); + myRootOutput->ClearAllRootVariables(); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/LEMuSR/src/lem4ShieldHit.cc b/geant4/LEMuSR/src/lem4ShieldHit.cc new file mode 100644 index 0000000..cd8f649 --- /dev/null +++ b/geant4/LEMuSR/src/lem4ShieldHit.cc @@ -0,0 +1,99 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "lem4ShieldHit.hh" +#include "G4UnitsTable.hh" +#include "G4VVisManager.hh" +#include "G4Circle.hh" +#include "G4Colour.hh" +#include "G4VisAttributes.hh" +#include +#include +#include "G4ios.hh" +#include "G4MagneticField.hh" +#include "G4FieldManager.hh" +#include "G4TransportationManager.hh" +#include +#include "globals.hh" +#include "G4Transform3D.hh" +#include "G4ProcessManager.hh" +#include "G4Track.hh" +#include "G4ThreeVector.hh" +#include "G4Event.hh" +#include "G4EventManager.hh" +#include "G4RunManager.hh" +#include "G4Run.hh" + +G4Allocator lem4ShieldHitAllocator; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ShieldHit::lem4ShieldHit() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ShieldHit::~lem4ShieldHit() {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ShieldHit::lem4ShieldHit(const lem4ShieldHit& right) + : G4VHit() +{ + + trackID = right.trackID; + particle_name = right.particle_name; + edep = right.edep; + pos = right.pos; + pol = right.pol; + logicalvolume = right.logicalvolume; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +const lem4ShieldHit& lem4ShieldHit::operator=(const lem4ShieldHit& right) +{ + + trackID = right.trackID; + particle_name = right.particle_name; + edep = right.edep; + pos = right.pos; + pol = right.pol; + logicalvolume=right.logicalvolume; + return *this; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4int lem4ShieldHit::operator==(const lem4ShieldHit& right) const +{ + return (this==&right) ? 1 : 0; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4ShieldHit::Draw() +{ + G4VVisManager* pVVisManager = G4VVisManager::GetConcreteInstance(); + if(pVVisManager) + { + G4Circle circle(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 lem4ShieldHit::Print() +{ + + G4EventManager* fEventManager = G4EventManager::GetEventManager(); + fEventManager->AbortCurrentEvent(); + +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + diff --git a/geant4/LEMuSR/src/lem4ShieldSD.cc b/geant4/LEMuSR/src/lem4ShieldSD.cc new file mode 100644 index 0000000..905520c --- /dev/null +++ b/geant4/LEMuSR/src/lem4ShieldSD.cc @@ -0,0 +1,82 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "lem4ShieldSD.hh" +#include "G4HCofThisEvent.hh" +#include "G4Step.hh" +#include "G4ThreeVector.hh" +#include "G4SDManager.hh" +#include "G4ios.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + /// + /// The Shield sensitive detector. This object collects the + /// information of a particle's interaction with this part + /// of the geometry + /// + +lem4ShieldSD::lem4ShieldSD(G4String name) +:G4VSensitiveDetector(name) +{ + G4String HCname; + collectionName.insert(HCname="shieldCollection"); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4ShieldSD::~lem4ShieldSD(){ } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4ShieldSD::Initialize(G4HCofThisEvent* HCE) +{ + shieldCollection = new lem4ShieldHitsCollection + (SensitiveDetectorName,collectionName[0]); + static G4int HCID = -1; + if(HCID<0) + { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); } + HCE->AddHitsCollection( HCID, shieldCollection ); +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +G4bool lem4ShieldSD::ProcessHits(G4Step* aStep,G4TouchableHistory*) +{ + /// + /// Collect the information regarding the interaction of a particle + /// and the sensitive detector. + /// + G4double edep = aStep->GetTotalEnergyDeposit(); + + if(edep==0.) return false; + + lem4ShieldHit* newHit = new lem4ShieldHit(); + newHit->SetParticleName (aStep->GetTrack()->GetDefinition()->GetParticleName()); + newHit->SetTrackID (aStep->GetTrack()->GetTrackID()); + newHit->SetEdep (edep); + newHit->SetPos (aStep->GetPostStepPoint()->GetPosition()); + newHit->SetPol (aStep->GetTrack()->GetPolarization()); + newHit->SetLogVolName (aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetName()); + shieldCollection->insert( newHit ); + + newHit->Print(); + newHit->Draw(); + + + + return true; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4ShieldSD::EndOfEvent(G4HCofThisEvent*) +{ + if (verboseLevel>0) { + G4int NbHits = shieldCollection->entries(); + G4cout << "\n-------->Hits Collection: in this event they are " << NbHits + << " hits: " << G4endl; + for (G4int i=0;iPrint(); + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + diff --git a/geant4/LEMuSR/src/lem4SteppingAction.cc b/geant4/LEMuSR/src/lem4SteppingAction.cc new file mode 100644 index 0000000..840dea1 --- /dev/null +++ b/geant4/LEMuSR/src/lem4SteppingAction.cc @@ -0,0 +1,337 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +#include "lem4SteppingAction.hh" +//cksdel#include "lem4RunAction.hh" +#include "G4SteppingManager.hh" +#include "G4UnitsTable.hh" +#include "lem4RootOutput.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 "lem4ErrorMessage.hh" +#include "lem4Parameters.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4SteppingAction::lem4SteppingAction() { + pointer=this; + // oldEventID=-1000; + muPlusProcessManager=NULL; + multipleToCoulombScatteringIsPossible=false; + coulombScatteringIsActive=false; + multipleScatteringIndex=1000; + coulombScatteringIndex=1000; + boolIsAnySpecialSaveVolumeDefined = false; + lastActualVolume="Unset"; +} + +lem4SteppingAction::~lem4SteppingAction() { +} + + + lem4SteppingAction* lem4SteppingAction::pointer=0; + lem4SteppingAction* lem4SteppingAction::GetInstance() +{ + return pointer; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4SteppingAction::DoAtTheBeginningOfEvent() { + // G4cout<<"lem4SteppingAction::DoAtTheBeginningOfEvent: KAMIL"<GetCurrentEvent()->GetEventID(); + // nrOfSteps++; + // if (debugOldEventID != debugEventID) { + // G4cout<<"lem4SteppingAction: Event nr. "<GetTrack(); + // if (nrOfSteps > 100000) { + // G4cout<<"lem4SteppingAction.cc : event nr. "<100000 ==> KILL THE CURRENT TRACK with" + // <GetDynamicParticle()->GetDefinition()->GetParticleName() <<" at" + // <GetPostStepPoint()->GetPosition() < 100000 ==> TRACK KILLED",true); + // + // G4cout<<"lem4SteppingAction.cc : event nr. "< 100000 ==> KILL THIS TRACK: " <GetDynamicParticle()->GetDefinition()->GetParticleName() <<" at" + // <GetPostStepPoint()->GetPosition() + // <<" E_kin_vertex="<GetVertexKineticEnergy()/MeV<<" MeV."<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); + } + + + if (aTrack->GetDefinition()) { + G4ParticleDefinition* p_definition = aTrack->GetDynamicParticle()->GetDefinition(); + G4String p_name = p_definition->GetParticleName(); + // G4ProcessManager* p_manager = p_definition->GetProcessManager(); + // G4String p_name=aTrack->GetDynamicParticle()->GetDefinition()->GetParticleName(); + G4String actualVolume=aTrack->GetVolume()->GetLogicalVolume()->GetName(); + + + +//! First Commented by TS. Currently not in use! or not strictly necessary + + + // This are the data just for the radioactive decay (when using the radioactive source): + if ((lem4Parameters::boolG4GeneralParticleSource)) { + // &&(!radioactiveElectronAlreadySavedInThisEvent)) { + if (aTrack->GetTrackID() != 1 ){ + if (aTrack->GetCreatorProcess()->GetProcessName() == "RadioactiveDecay") { + // if (aStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() == "RadioactiveDecay") { did not work out + // emitted particles + // if (fTrack->GetDefinition()->GetParticleType() != "nucleus") { + // G4String particleName = fTrack->GetDefinition()->GetParticleName(); + if (aTrack->GetDefinition()->GetParticleName()=="e-") { + // G4double particleName = G4double (fTrack->GetDefinition()->GetPDGEncoding()); + // G4double time = fStep->GetPreStepPoint()->GetGlobalTime() ; + // //- fStep->GetPreStepPoint()->GetLocalTime(); + // G4double weight = fStep->GetPreStepPoint()->GetWeight() ; + // G4double energy = fStep->GetPreStepPoint()->GetKineticEnergy(); + // + // exrdmAnalysisManager::getInstance()->AddParticle(particleName, energy, weight, time); + //hovno + // G4cout<<"aTrack->GetCurrentStepNumber()="<GetCurrentStepNumber()<<" "; + if (aTrack->GetCurrentStepNumber()==1) { + G4double electron_kinetic_energy=aStep->GetPreStepPoint()->GetKineticEnergy(); + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + myRootOutput->htest4->Fill(electron_kinetic_energy); + // G4cout<<"First kinetic energy="<GetPreStepPoint()->GetKineticEnergy()<<" "; + // G4cout<<"GetKineticEnergy()="<GetKineticEnergy()<GetPreStepPoint()->GetMomentum(); + // G4ThreeVector radioactivePositronMomentum = aTrack->GetMomentum(); + // G4double Apx=radioactivePositronMomentum.x(); + // G4double Apy=radioactivePositronMomentum.y(); + // G4double Apz=radioactivePositronMomentum.z(); + // lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + // myRootOutput->SetInitialMuonParameters(0,0,0,Apx,Apy,Apz,0,0,0); + // radioactiveElectronAlreadySavedInThisEvent=true; + // // G4cout<<"Apx,Apy,Apz="<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()/keV; + G4double x_save=aStep->GetPreStepPoint()->GetPosition().x()/mm; + G4double y_save=aStep->GetPreStepPoint()->GetPosition().y()/mm; + G4double z_save=aStep->GetPreStepPoint()->GetPosition().z()/mm; + G4double px_save=aStep->GetPreStepPoint()->GetMomentum().x()/MeV; + G4double py_save=aStep->GetPreStepPoint()->GetMomentum().y()/MeV; + G4double pz_save=aStep->GetPreStepPoint()->GetMomentum().z()/MeV; + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,px_save,py_save,pz_save); + } + } + } + + + // Store information about when mu+ or Mu enter the target for the fist time + // in a given event (i.e. the code has to be called just once during the event). + if ((p_name == "mu+") || (p_name == "Mu")) { + if (actualVolume=="log_target") { + // G4RunManager* fRunManager = G4RunManager::GetRunManager(); + // G4int eventID = fRunManager->GetCurrentEvent()->GetEventID(); + // if (oldEventID != eventID) { + // oldEventID=eventID; + if (!muAlreadyWasInTargetInThisEvent) { + muAlreadyWasInTargetInThisEvent=true; + lem4RootOutput* myRootOutput = lem4RootOutput::GetRootInstance(); + myRootOutput->SetPolInTarget(aTrack->GetPolarization()); + myRootOutput->SetTimeInTarget(aTrack->GetGlobalTime()); + //G4cout<<"particle = "<GetProcessManager(); + G4ProcessVector* v = muPlusProcessManager->GetProcessList(); + size_t np = v->size(); + for(size_t i=0; iGetProcessName(); + if( (*v)[i]->GetProcessName() == "msc" ) { + multipleScatteringIndex = i; + } + if( (*v)[i]->GetProcessName() == "eCoulombScat" ) { + coulombScatteringIndex = i; + } + } + // std::cout<SetProcessActivation(multipleScatteringIndex, true); + muPlusProcessManager->SetProcessActivation(coulombScatteringIndex, false); + } + else { // At least one of the "multiple scattering" and "coulomb scattering" was not found. + // Perhaps they were not defined in the physics list ? + if (lem4Parameters::GetInstance()->GetMyTypeOfProcesses()=="coulombAndMultiple") { + if (coulombScatteringIndex==1000) {std::cout<<"lem4SteppingAction: Coulomb scattering process for mu+ not found!"<lem4Error(SERIOUS, + "lem4SteppingAction: It is not possible to switch between Coulomb and multiple scattering.",true); + + muPlusProcessManager->DumpInfo(); + } + } + } + + // switch between G4MultipleScattering and G4CoulombScattering + if (multipleToCoulombScatteringIsPossible) { + if (strncmp(actualVolume.c_str(),"log_coulomb",11)==0) { + if (!coulombScatteringIsActive) { + // std::cout<<"Activating coulomb; volume= "<SetProcessActivation(multipleScatteringIndex, false); + muPlusProcessManager->SetProcessActivation(coulombScatteringIndex, true); + coulombScatteringIsActive=true; + } + } + else if (coulombScatteringIsActive) { + // std::cout<<"Deactivating coulomb; volume= "<SetProcessActivation(multipleScatteringIndex, true); + muPlusProcessManager->SetProcessActivation(coulombScatteringIndex, false); + coulombScatteringIsActive=false; + } + } + + // Store the information about the muon decay into the Root Class. + // Pick up process "DecayWithSpin": + const G4VProcess* process = aStep->GetPostStepPoint()->GetProcessDefinedStep(); + if (process!=NULL) { + G4String processName = process->GetProcessName(); + if (processName=="DecayWithSpin") { + lem4Parameters::field_DecayWithSpin=true; + // std::cout<<"lem4SteppingAction: DecayWithSpin"<GetGlobalTime(); + myRootOutput->SetDecayTime(timeOfDecay_tmp); + myRootOutput->SetDecayPolarisation(aTrack->GetPolarization()); + G4ThreeVector positionOfDecay_tmp = aStep->GetPostStepPoint()->GetPosition(); + myRootOutput->SetDecayPosition(positionOfDecay_tmp); + // store the detector ID in which the muon decayed + myRootOutput->SetDecayDetectorID(actualVolume); + // store the information about the magnetic field at the place where the muon decays + G4double BFieldAtOrigin[6] = {0.,0.,0.,0.,0.,0.}; + G4double PointOfDecay[4] ={positionOfDecay_tmp.x(),positionOfDecay_tmp.y(),positionOfDecay_tmp.z(),timeOfDecay_tmp}; + G4FieldManager *fMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager(); + if (fMgr!=NULL) { + if(!fMgr->DoesFieldChangeEnergy()) { //then we have a magnetic field + fMgr->GetDetectorField()->GetFieldValue(PointOfDecay,BFieldAtOrigin); + myRootOutput->SetBField(BFieldAtOrigin); + } + else{ + } + } + + // store the information about the emerging positron + G4TrackVector* secondary = fpSteppingManager->GetSecondary(); + G4int n_secondaries= (*secondary).size(); + for (G4int i=0; iGetDefinition()->GetParticleName()) == "e+" ) { + myRootOutput->SetInitialPositronMomentum((*secondary)[i]->GetMomentum()); + + // These two times seem to be equivalent: (*secondary)[i]->GetGlobalTime() + // aTrack->GetGlobalTime() + // G4cout <<"Positron initial momentum="<<(*secondary)[i]->GetMomentum()<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 lem4SteppingAction::SetLogicalVolumeAsSpecialSaveVolume(G4String logicName, G4int volumeID) { + boolIsAnySpecialSaveVolumeDefined = true; + saveVolumeMapping[logicName]=volumeID; +} diff --git a/geant4/LEMuSR/src/lem4SteppingVerbose.cc b/geant4/LEMuSR/src/lem4SteppingVerbose.cc new file mode 100644 index 0000000..68d3fc2 --- /dev/null +++ b/geant4/LEMuSR/src/lem4SteppingVerbose.cc @@ -0,0 +1,151 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#include "lem4SteppingVerbose.hh" + +#include "G4SteppingManager.hh" +#include "G4UnitsTable.hh" + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4SteppingVerbose::lem4SteppingVerbose() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4SteppingVerbose::~lem4SteppingVerbose() +{} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4SteppingVerbose::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 lem4SteppingVerbose::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...... diff --git a/geant4/LEMuSR/src/lem4TabulatedElementField2D.cc b/geant4/LEMuSR/src/lem4TabulatedElementField2D.cc new file mode 100644 index 0000000..e0778c3 --- /dev/null +++ b/geant4/LEMuSR/src/lem4TabulatedElementField2D.cc @@ -0,0 +1,402 @@ +#include "lem4TabulatedElementField2D.hh" +//#include "lem4EventAction.hh" // cks: just to get the event nr. if needed. +#include "lem4Parameters.hh" /// Comment by TS - is this really needed? +///#include "lem4ErrorMessage.hh" -> PUT SOME LIMIT DURING WRONG UNIT CONVERSIONS +#include "G4UnitsTable.hh" +#include + +//lem4TabulatedElementField2D::lem4TabulatedElementField2D(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), +// ffieldValue(fieldValue) + + +///G4double lenUnit, G4double fieldNormalisation, + +lem4TabulatedElementField2D::lem4TabulatedElementField2D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), +fldType(fieldType), ffieldValue(fieldValue) +{ + G4cout << "\n-----------------------------------------------------------"; + + // The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla) + // NOTE: Should be the same as EMfieldUnit defined in DetectorConstruction.cc! + if (fldType == 'E') {G4cout << ("\n Electric field 2D"); fUnit = "kV/mm"; fieUnit= kilovolt/mm;} + else if (fldType == 'B') {G4cout << ("\n Magnetic field 2D"); fUnit = "T"; fieUnit = tesla;} + + G4cout << "\n-----------------------------------------------------------" << G4endl; + G4cout << "\n ---> Reading the "<< fldType <<"-field map from " << filename << " ... " << G4endl; + + // Open the file for reading + std::ifstream file( filename ); + + ///char buffer[256]; //not needed with file peek + // Ignore first blank line + // file.getline(buffer,256); + + // Read table dimensions + //G4String lUnit, fUnit; + //G4double cFact; + //fieldNormalisation = cFact; + //lenUnit = 1*cm, + //fieldNormalisation = 0.00001; + // 101 1 1001 cm tesla 0.00001 + //G4double fieldNorm, lenNorm; + //G4cout<< lUnit << " XXXX " << fUnit << G4endl; + //G4cout<< "Tesla in G4 is "<< tesla << ", while kV/mm is " << kilovolt/mm << G4endl; + + ///file >> nr >> nDummy >> nz >> lUnit >> fUnit >> fieldNormalisation; // Note dodgy order + // Read the file header + file >> nr >> nz >> lUnit >> fieldNormalisation; + // Add manually the length unit and norm. factor in the field-map file if missing! + + G4cout << " The grid consists of [" << nr << " x " << nz << "] r and z values" << G4endl; + + // Get the field map LENGTH unit and its value in G4 native units. + lenUnit = G4UnitDefinition::GetValueOf(lUnit); + G4cout << " Field length unit = " << lUnit << ", in G4 native units = " << lenUnit << G4endl; + + /* + // Get the EM field AMPLITUDE unit and its value in G4 native units. + // Separate the E case from B, since E presents complex units (e.g. kV/mm) + if (fldType == 'E') { + G4double volNorm, lenNorm; + std::string::size_type pos = fUnit.find("/"); + std::string v_unit = fUnit.substr(0, pos); // Voltage unit + std::string l_unit = fUnit.substr(pos+1); // Length unit + + volNorm = G4UnitDefinition::GetValueOf(v_unit); + lenNorm = G4UnitDefinition::GetValueOf(l_unit); + ///G4cout<< v_unit << " " << l_unit << G4endl; + ///G4cout<< volNorm << " " << lenNorm << G4endl; + + fNorm = volNorm/lenNorm; // Electric field unit = [Voltage]/[Length] + G4cout << " Electric field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl; + } + + else if (fldType == 'B') { + fNorm = G4UnitDefinition::GetValueOf(fUnit); + G4cout << " Magnetic field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl; + } + */ + G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl; + + + ///fieldValue = fieldValue*fNorm/tesla; // Convert + + ///fieldValue = fieldValue/fNorm; // Convert + //double lenUnit = meter; + //double fieldUnit= tesla; + + //lenNorm = G4UnitDefinition::GetValueOf(leng_unit); + //lenNorm = G4UnitDefinition::GetValueOf(lUnit); + //fieldNorm = G4UnitDefinition::GetValueOf(fNorm); + + + //G4double GetNewDoubleValue(const char*) + //fieldNorm = G4UIcmdWithADouble::SetNewDoubleValue(fNorm); + //fieldNorm = G4UIcmdWithADoubleAndUnit::GetNewUnitValue(fNorm); //GetNewDoubleRawValue(fNorm); + + +/* + + G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit + G4cout << " Field Unit = " << leng_unit << G4endl; //lUnit + + G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit + G4cout << " Conv. Unit = " << lenNorm << G4endl; + G4cout << " Field Unit = " << fNorm << G4endl; + G4cout << " Conv. Unit = " << fieldNorm << G4endl; + G4cout << " Norm. Fact = " << cFact << G4endl; + + */ + /// FIELD NORMALISATION FACTOR INSIDE FIELD MAPS IS SUCH THAT (MAX. FIELD VALUE)*FNORMFACT = 1! + /// CHANGE SIGN TO REVERT FIELD DIRECTION! THE FIELD MAP HAS NO UNITS. + + ///G4cout << " Field set to "<< fieldValue/fNorm << " " << fUnit << G4endl; ///???? + ///G4cout << " Field set to "<< fieldValue*fNorm << " " << fUnit << G4endl; + ///G4cout << " Field set to "<< fieldValue << " " << fUnit << G4endl; + + + // Set up storage space for the 2D table + rField.resize(nr); + zField.resize(nr); + int ir, iz; + for (ir=0; ir> rval >> zval >> fr >> fz; + //G4cout << rval << " v " << zval << " v " << fr << " v " << fz << G4endl; + if ( ir==0 && iz==0 ) { + minimumr = rval * lenUnit; + minimumz = zval * lenUnit; + } + rField[ir][iz] = fr*fieldNormalisation; + zField[ir][iz] = fz*fieldNormalisation; + } + } + file.close(); + + maximumr = rval * lenUnit; + maximumz = zval * lenUnit; + + G4String volumeName = logVolume->GetName().substr(4); + + G4cout << " ... done reading " << G4endl; + G4cout << "\n ---> " << fldType << "-field in volume "<< volumeName + << " set to: " << fieldValue/fieUnit << " " << fUnit << G4endl; + G4cout << "\n ---> Assumed order: r, z, "<< fldType <<"r, "<< fldType <<"z " + << "\n Min values r, z: " + << minimumr/cm << " " << minimumz/cm << " cm " + << "\n Max values r, z: " + << maximumr/cm << " " << maximumz/cm << " cm" << G4endl; + + + // Should really check that the limits are not the wrong way around. + if (maximumr < minimumr) {Invert("r");} + if (maximumz < minimumz) {Invert("z");} + + if ((maximumr < minimumr) || (maximumz < minimumz)) { + G4cout << "\n ---> After reordering:" + << "\n Min values r, z: " + << minimumr/cm << " " << minimumz/cm << " cm " + << "\n Max values r, z: " + << maximumr/cm << " " << maximumz/cm << " cm " << G4endl; + } + + dr = maximumr - minimumr; + dz = maximumz - minimumz; + + G4cout << " Range of values: " + << dr/cm << " cm (in r) and " << dz/cm << " cm (in z)" + << "\n-----------------------------------------------------------\n" << G4endl; + +} + + +void lem4TabulatedElementField2D::addFieldValue(const G4double point[4], + G4double *field ) const +{ + + ///G4cout << "addFieldValue is being called" << G4endl; + G4double EMf[3]; // EM field value obtained from the field map + ///TS + ///G4cout<<"TTfileTT Global coord. ("<GetName()<0) ? 1.:-1.; + +// G4cout<<"Global points= "<evNrKriz) std::cout<<"bol som tu"<(rdindex); + int zindex = static_cast(zdindex); + + //cks The following check is necessary - even though rindex 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 ((rindex<0)||(rindex>(nr-2))) { + std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="< > rFieldTemp(rField); + std::vector< std::vector< double > > zFieldTemp(zField); + G4bool invertR=false; + G4bool invertZ=false; + + G4cout<<"Check that the lem4TabulatedElementField2D::Invert() function works properly!"< +//using namespace std; + + ///G4double G4UIcmdWithADouble::GetNewDoubleValue(const char* paramString) + + + ///**********************/// + //char str[] ="- This, a sample string."; + ///char str[] ="kilovolt/mm"; + ///char * pch; + ///printf ("Splitting string \"%s\" into tokens:\n",str); + //pch = strtok (str," ,.-"); + ///pch = strtok (str,"*/"); + ///while (pch != NULL) + ///{ + ///printf ("%s\n",pch); + //pch = strtok (NULL, " ,.-"); + ///pch = strtok (NULL,"*/"); + ///} + + ///**********************/// + + //char str[] ="kilovolt/mm"; + //char * pch; + //printf ("Splitting string \"%s\" into tokens:\n",str); + //pch = strtok (str," ,.-"); + //pch = strtok (str,"/"); + //for (int i=0; i<2; x++) { + //while (pch != NULL) + //{ +// printf ("%s\n",pch); + // pch[] + //} + //pch = strtok (NULL,"*/"); + //} + + //int x; + //for(x=1; x <= 100; x++) { + //printf("%d ", x); + //} + + ///**********************/// + // Only if field is E! + ///char str[] = "This is a sample string kilovolt/mm"; + ///char * pch; + ///printf ("Looking for the '/' character in \"%s\"...\n",str); + ///pch = strchr(str,'/'); + ///while (pch!=NULL) + ///{ + ///printf ("found at %d\n",pch-str+1); + ///pch=strchr(pch+1,'s'); + ///} + + ///**********************/// + diff --git a/geant4/LEMuSR/src/lem4TabulatedElementField2Df.cc b/geant4/LEMuSR/src/lem4TabulatedElementField2Df.cc new file mode 100644 index 0000000..41d90c9 --- /dev/null +++ b/geant4/LEMuSR/src/lem4TabulatedElementField2Df.cc @@ -0,0 +1,420 @@ +#include "lem4TabulatedElementField2Df.hh" +//#include "lem4EventAction.hh" // cks: just to get the event nr. if needed. +#include "lem4Parameters.hh" /// Comment by TS - is this really needed? +///#include "lem4ErrorMessage.hh" -> PUT SOME LIMIT DURING WRONG UNIT CONVERSIONS +#include "G4UnitsTable.hh" +#include + +/// TS ATTENTION: Since B and E are respectively Axial and Polar vectors, their +/// transformation upon inversion of coordinates are DIFFERENT!! +/// A Polar vector reverses sign when the coordinate axes are reversed. +/// An Axial vector does not reverse sign when the coordinate axes are reversed. + + +/// This version of the 2d axial field map contains an additional two-fold symmetry at z=0. +/// Therefore the maps read by this file (ending in 1Dq) are 1/4 of the whole plane map. + +//lem4TabulatedElementField2D::lem4TabulatedElementField2D(const char* filename, G4double fieldValue, G4double lenUnit, G4double fieldNormalisation, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), +// ffieldValue(fieldValue) + + +///G4double lenUnit, G4double fieldNormalisation, + +lem4TabulatedElementField2Df::lem4TabulatedElementField2Df(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), +fldType(fieldType), ffieldValue(fieldValue) +{ + G4cout << "\n-----------------------------------------------------------"; + + // The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla) + // NOTE: Should be the same as EMfieldUnit defined in DetectorConstruction.cc! + if (fldType == 'E') {G4cout << ("\n Electric field 2D - Folded (r, z > 0)"); fUnit = "kV/mm"; fieUnit= kilovolt/mm;} + else if (fldType == 'B') {G4cout << ("\n Magnetic field 2D - Folded (r, z > 0)"); fUnit = "T"; fieUnit = tesla;} + + G4cout << "\n-----------------------------------------------------------" << G4endl; + G4cout << "\n ---> Reading the "<< fldType <<"-field map from " << filename << " ... " << G4endl; + + // Open the file for reading + std::ifstream file( filename ); + + ///char buffer[256]; //not needed with file peek + // Ignore first blank line + // file.getline(buffer,256); + + // Read table dimensions + //G4String lUnit, fUnit; + //G4double cFact; + //fieldNormalisation = cFact; + //lenUnit = 1*cm, + //fieldNormalisation = 0.00001; + // 101 1 1001 cm tesla 0.00001 + //G4double fieldNorm, lenNorm; + //G4cout<< lUnit << " XXXX " << fUnit << G4endl; + //G4cout<< "Tesla in G4 is "<< tesla << ", while kV/mm is " << kilovolt/mm << G4endl; + + ///file >> nr >> nDummy >> nz >> lUnit >> fUnit >> fieldNormalisation; // Note dodgy order + // Read the file header + file >> nr >> nz >> lUnit >> fieldNormalisation; + // Add manually the length unit and norm. factor in the field-map file if missing! + + G4cout << " The grid consists of [" << nr << " x " << nz << "] r and z values" << G4endl; + + // Get the field map LENGTH unit and its value in G4 native units. + lenUnit = G4UnitDefinition::GetValueOf(lUnit); + G4cout << " Field length unit = " << lUnit << ", in G4 native units = " << lenUnit << G4endl; + + /* + // Get the EM field AMPLITUDE unit and its value in G4 native units. + // Separate the E case from B, since E presents complex units (e.g. kV/mm) + if (fldType == 'E') { + G4double volNorm, lenNorm; + std::string::size_type pos = fUnit.find("/"); + std::string v_unit = fUnit.substr(0, pos); // Voltage unit + std::string l_unit = fUnit.substr(pos+1); // Length unit + + volNorm = G4UnitDefinition::GetValueOf(v_unit); + lenNorm = G4UnitDefinition::GetValueOf(l_unit); + ///G4cout<< v_unit << " " << l_unit << G4endl; + ///G4cout<< volNorm << " " << lenNorm << G4endl; + + fNorm = volNorm/lenNorm; // Electric field unit = [Voltage]/[Length] + G4cout << " Electric field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl; + } + + else if (fldType == 'B') { + fNorm = G4UnitDefinition::GetValueOf(fUnit); + G4cout << " Magnetic field unit = " << fUnit << ", in G4 native units = " << fNorm << G4endl; + } + */ + G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl; + + + ///fieldValue = fieldValue*fNorm/tesla; // Convert + + ///fieldValue = fieldValue/fNorm; // Convert + //double lenUnit = meter; + //double fieldUnit= tesla; + + //lenNorm = G4UnitDefinition::GetValueOf(leng_unit); + //lenNorm = G4UnitDefinition::GetValueOf(lUnit); + //fieldNorm = G4UnitDefinition::GetValueOf(fNorm); + + + //G4double GetNewDoubleValue(const char*) + //fieldNorm = G4UIcmdWithADouble::SetNewDoubleValue(fNorm); + //fieldNorm = G4UIcmdWithADoubleAndUnit::GetNewUnitValue(fNorm); //GetNewDoubleRawValue(fNorm); + + +/* + + G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit + G4cout << " Field Unit = " << leng_unit << G4endl; //lUnit + + G4cout << " Length Unit = " << leng_unit << G4endl; //lUnit + G4cout << " Conv. Unit = " << lenNorm << G4endl; + G4cout << " Field Unit = " << fNorm << G4endl; + G4cout << " Conv. Unit = " << fieldNorm << G4endl; + G4cout << " Norm. Fact = " << cFact << G4endl; + + */ + /// FIELD NORMALISATION FACTOR INSIDE FIELD MAPS IS SUCH THAT (MAX. FIELD VALUE)*FNORMFACT = 1! + /// CHANGE SIGN TO REVERT FIELD DIRECTION! THE FIELD MAP HAS NO UNITS. + + ///G4cout << " Field set to "<< fieldValue/fNorm << " " << fUnit << G4endl; ///???? + ///G4cout << " Field set to "<< fieldValue*fNorm << " " << fUnit << G4endl; + ///G4cout << " Field set to "<< fieldValue << " " << fUnit << G4endl; + + // Set up storage space for the 2D table + rField.resize(nr); + zField.resize(nr); + int ir, iz; + for (ir=0; ir> rval >> zval >> fr >> fz; + if ( ir==0 && iz==0 ) { + minimumr = rval * lenUnit; + minimumz = zval * lenUnit; + } + rField[ir][iz] = fr*fieldNormalisation; + zField[ir][iz] = fz*fieldNormalisation; + } + } + file.close(); + + maximumr = rval * lenUnit; + maximumz = zval * lenUnit; + + G4String volumeName = logVolume->GetName().substr(4); + + G4cout << " ... done reading " << G4endl; + G4cout << "\n ---> " << fldType << "-field in volume "<< volumeName + << " set to: " << fieldValue/fieUnit << " " << fUnit << G4endl; + G4cout << "\n ---> Assumed order: r, z, "<< fldType <<"r, "<< fldType <<"z " + << "\n Min values r, z: " + << minimumr/cm << " " << minimumz/cm << " cm " + << "\n Max values r, z: " + << maximumr/cm << " " << maximumz/cm << " cm" << G4endl; + + + // Should really check that the limits are not the wrong way around. + if (maximumr < minimumr) {Invert("r");} + if (maximumz < minimumz) {Invert("z");} + + if ((maximumr < minimumr) || (maximumz < minimumz)) { + G4cout << "\n ---> After reordering:" + << "\n Min values r, z: " + << minimumr/cm << " " << minimumz/cm << " cm " + << "\n Max values r, z: " + << maximumr/cm << " " << maximumz/cm << " cm " << G4endl; + } + + dr = maximumr - minimumr; + dz = maximumz - minimumz; + + G4cout << " Range of values: " + << dr/cm << " cm (in r) and " << dz/cm << " cm (in z)" + << "\n-----------------------------------------------------------\n" << G4endl; + +} + + +void lem4TabulatedElementField2Df::addFieldValue(const G4double point[4], + G4double *field ) const +{ + G4double EMf[3]; // EM field value obtained from the field map + G4ThreeVector global(point[0],point[1],point[2]); + G4ThreeVector local; + + local = global2local.TransformPoint(global); + + /// The folded 2D field map contains only positive r and z! + double r = sqrt(local.x()*local.x()+local.y()*local.y()); + double z = fabs(local.z()); + /// Store the z_sign and use it when extending map to negative z-values + double z_sign =(local.z()>0) ? 1.:-1.; + + ///double z = 0; /// An alternative approach + /// E is a Polar vectorn, while B is an axial one (pseudovector)! TS + ///if (fldType == 'E') {z = local.z();} + ///else if (fldType == 'B') {z = fabs(local.z());} + ///G4cout<<"Global points= "<< z << G4endl; + +// G4cout<<"Global points= "<evNrKriz) std::cout<<"bol som tu"<(rdindex); + int zindex = static_cast(zdindex); + + //cks The following check is necessary - even though rindex 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 ((rindex<0)||(rindex>(nr-2))) { + std::cout<<"SERIOUS PROBLEM: rindex out of range! rindex="< > rFieldTemp(rField); + std::vector< std::vector< double > > zFieldTemp(zField); + G4bool invertR=false; + G4bool invertZ=false; + + G4cout<<"Check that the lem4TabulatedElementField2Df::Invert() function works properly!"< +//using namespace std; + + ///G4double G4UIcmdWithADouble::GetNewDoubleValue(const char* paramString) + + + ///**********************/// + //char str[] ="- This, a sample string."; + ///char str[] ="kilovolt/mm"; + ///char * pch; + ///printf ("Splitting string \"%s\" into tokens:\n",str); + //pch = strtok (str," ,.-"); + ///pch = strtok (str,"*/"); + ///while (pch != NULL) + ///{ + ///printf ("%s\n",pch); + //pch = strtok (NULL, " ,.-"); + ///pch = strtok (NULL,"*/"); + ///} + + ///**********************/// + + //char str[] ="kilovolt/mm"; + //char * pch; + //printf ("Splitting string \"%s\" into tokens:\n",str); + //pch = strtok (str," ,.-"); + //pch = strtok (str,"/"); + //for (int i=0; i<2; x++) { + //while (pch != NULL) + //{ +// printf ("%s\n",pch); + // pch[] + //} + //pch = strtok (NULL,"*/"); + //} + + //int x; + //for(x=1; x <= 100; x++) { + //printf("%d ", x); + //} + + ///**********************/// + // Only if field is E! + ///char str[] = "This is a sample string kilovolt/mm"; + ///char * pch; + ///printf ("Looking for the '/' character in \"%s\"...\n",str); + ///pch = strchr(str,'/'); + ///while (pch!=NULL) + ///{ + ///printf ("found at %d\n",pch-str+1); + ///pch=strchr(pch+1,'s'); + ///} + + ///**********************/// + diff --git a/geant4/LEMuSR/src/lem4TabulatedElementField3D.cc b/geant4/LEMuSR/src/lem4TabulatedElementField3D.cc new file mode 100644 index 0000000..75db7b3 --- /dev/null +++ b/geant4/LEMuSR/src/lem4TabulatedElementField3D.cc @@ -0,0 +1,327 @@ +#include "lem4TabulatedElementField3D.hh" +#include "lem4Parameters.hh" + +#include "G4UnitsTable.hh" +#include + +lem4TabulatedElementField3D::lem4TabulatedElementField3D(const char* filename, const char fieldType, G4double fieldValue, G4LogicalVolume* logVolume, G4ThreeVector positionOfTheCenter) : F04ElementField(positionOfTheCenter, logVolume), +fldType(fieldType), ffieldValue(fieldValue) +{ + G4cout << "\n-----------------------------------------------------------"; + + // The DEFAULT user-defined field units for E and B (kilovolt/mm and tesla) + // NOTE: Should be the same as EMfieldUnit defined in DetectorConstruction.cc! + if (fldType == 'E') {G4cout << ("\n Electric field 3D"); fUnit = "kV/mm"; fieUnit= kilovolt/mm;} + else if (fldType == 'B') {G4cout << ("\n Magnetic field 3D"); fUnit = "T"; fieUnit = tesla;} + + G4cout << "\n-----------------------------------------------------------" << G4endl; + G4cout << "\n ---> Reading the "<< fldType <<"-field map from " << filename << " ... " << G4endl; + + // Open the file for reading + std::ifstream file( filename ); + + // Read the first line of the file header (contains the metadata) + char buffer[256]; + file.getline(buffer,256); + + //G4cout << "The first line of file "<< filename << " is:\n"<< buffer << G4endl; + // Read the file header + ///file >> nx >> ny >> nz >> lUnit >> fieldNormalisation; // OLD VERSION + + // Read the number of arguments and decide filetype - 3 or 6 columns + int n_arg = sscanf (buffer,"%d %d %d %s %lf %lf %lf %lf %lf %lf %lf", + &nx, &ny, &nz, lUnit, &fieldNormalisation, + &minimumx, &maximumx, &minimumy, &maximumy, &minimumz, &maximumz); + // Add manually the length unit and norm. factor in the field-map file if missing! + + //printf ("The following data were read in: \n"); + //printf ("nx = %d, ny = %d, nz = %d, Unit = %s, Norm = %f,\n minx = %f, maxx = %f, miny = %f, maxy = %f, minz = %f, maxz = %f\n",nx,ny,nz,lUnit,fieldNormalisation,minimumx, maximumx, minimumy, maximumy, minimumz, maximumz); + ///printf ("nx = %d, ny = %d, nz = %d, Unit = %s, Norm = %f,\n minx = %f\n",nx,ny,nz,lUnit,fieldNormalisation,minimumx); + //printf ("Number of args read: %d\n", n_arg); + + G4cout << " The grid consists of [" << nx << " x " << ny << " x " << nz << "] x, y, z values" << G4endl; + + // Get the field map LENGTH unit and its value in G4 native units. + lenUnit = G4UnitDefinition::GetValueOf(lUnit); + G4cout << " Field length unit = " << lUnit << ", in G4 native units = " << lenUnit << G4endl; + G4cout << " Field map normalisation factor = " << fieldNormalisation << G4endl; + + // Set up storage space for the 3D table: notice the order! + xField.resize(nx); + yField.resize(nx); + zField.resize(nx); + int ix, iy, iz; + for (ix=0; ix> xval >> yval >> zval >> fx >> fy >> fz; + //fileout << fx <<" \t"<< fy << " \t"<< fz << std::endl; // for rewriting files + if ( ix==0 && iy==0 && iz==0 ) { + minimumx = xval * lenUnit; + minimumy = yval * lenUnit; + minimumz = zval * lenUnit; + } + } + else if (n_arg == 11) { + file >> fx >> fy >> fz; // Read ONLY field values + } + + xField[ix][iy][iz] = fx*fieldNormalisation; + yField[ix][iy][iz] = fy*fieldNormalisation; + zField[ix][iy][iz] = fz*fieldNormalisation; + } + } + } + //fileout.close(); // for rewriting files + file.close(); + + if (n_arg == 5) { + maximumx = xval * lenUnit; + maximumy = yval * lenUnit; + maximumz = zval * lenUnit; + } + else if (n_arg == 11) { + minimumx = minimumx * lenUnit; + minimumy = minimumy * lenUnit; + minimumz = minimumz * lenUnit; + maximumx = maximumx * lenUnit; + maximumy = maximumy * lenUnit; + maximumz = maximumz * lenUnit; + } + + + G4String volumeName = logVolume->GetName().substr(4); + + G4cout << " ... done reading " << G4endl; + G4cout << "\n ---> " << fldType << "-field in volume "<< volumeName + << " set to: " << fieldValue/fieUnit << " " << fUnit << G4endl; + if (n_arg == 5) { G4cout << "\n ---> Assumed order (6 col.): x, y, z, ";} + else if (n_arg == 11) { G4cout << "\n ---> Assumed order (3 col.): ";} + G4cout << fldType <<"x, "<< fldType <<"y, "<< fldType <<"z " + << "\n Min values x, y, z: " + << minimumx/cm << " " << minimumy/cm << " " << minimumz/cm << " cm " + << "\n Max values x, y, z: " + << maximumx/cm << " " << maximumy/cm << " " << maximumz/cm << " cm" << G4endl; + + // Should really check that the limits are not the wrong way around. + if (maximumx < minimumx) {Invert("x");} + if (maximumy < minimumy) {Invert("y");} + if (maximumz < minimumz) {Invert("z");} + + if ((maximumx < minimumx) || (maximumy < minimumy)|| (maximumz < minimumz)) { + G4cout << "\n ---> After reordering:" + << "\n Min values x, y, z: " + << minimumx/cm << " " << minimumy/cm << " " << minimumz/cm << " cm " + << "\n Max values x, y, z: " + << maximumx/cm << " " << maximumy/cm << " " << maximumz/cm << " cm " << G4endl; + } + + dx = maximumx - minimumx; + dy = maximumy - minimumy; + dz = maximumz - minimumz; + + G4cout << " Range of values: " + << dx/cm << " cm (in x), " << dy/cm << " cm (in y) and " << dz/cm << " cm (in z)." + << "\n-----------------------------------------------------------\n" << G4endl; + +} + +void lem4TabulatedElementField3D::addFieldValue(const G4double point[4], + G4double *field ) const +{ + G4double EMf[3]; // EM field value obtained from the field map + + 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= "<minimumx && xminimumy && yminimumz && z(xdindex); + int yindex = static_cast(ydindex); + int zindex = static_cast(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="< " "Reading the field grid from " << filename << " ... " << G4endl; + positionOffset[0]=0; + positionOffset[1]=0; + positionOffset[2]=0; + positionOffset[3]=0; + std::ifstream file( filename ); // Open the file for reading. + + // Ignore first blank line + char buffer[256]; + // file.getline(buffer,256); + + // Read table dimensions + int nDummy; + file >> nx >> nDummy >> nz; // Note dodgy order + + G4cout << " [ Number of values x,z: " + << nx << " " << nz << " ] " + << G4endl; + + // Set up storage space for table + xField.resize( nx ); + zField.resize( nx ); + int ix, iz; + for (ix=0; ix> xval >> yval >> zval >> bx >> bz >> permeability; + if ( ix==0 && iz==0 ) { + minx = xval * lenUnit; + minz = zval * lenUnit; + } + //xField[ix][iy][iz] = bx * fieldUnit * fieldValue; + //zField[ix][iy][iz] = bz * fieldUnit * fieldValue; + // xField[ix][iy][iz] = bx * fieldValue; + // zField[ix][iy][iz] = bz * fieldValue; + xField[ix][iz] = bx*fieldNormalisation; + zField[ix][iz] = bz*fieldNormalisation; + } + } + file.close(); + + maxx = xval * lenUnit; + maxz = zval * lenUnit; + + G4cout << "\n ---> ... done reading " << G4endl; + + // G4cout << " Read values of field from file " << filename << G4endl; + G4cout << " ---> assumed the order: x, y, z, Bx, Bz " + << "\n ---> Min values x,y,z: " + << minx/cm << " " << minz/cm << " cm " + << "\n ---> Max values x,y,z: " + << maxx/cm << " " << maxz/cm << " cm " << G4endl; + + + // Should really check that the limits are not the wrong way around. + if (maxx < minx) {std::swap(maxx,minx); invertX = true;} + if (maxz < minz) {std::swap(maxz,minz); invertZ = true;} + G4cout << "\nAfter reordering if neccesary" + << "\n ---> Min values x,y,z: " + << minx/cm << " " << minz/cm << " cm " + << " \n ---> Max values x,y,z: " + << maxx/cm << " " << maxz/cm << " cm "; + + dx = maxx - minx; + dz = maxz - minz; + G4cout << "\n ---> Dif values x,z (range): " + << dx/cm << " " << dz/cm << " cm in z " + << "\n-----------------------------------------------------------" << G4endl; +} + +void lem4TabulatedField2D::GetFieldValue(const double point[4], + double *Bfield ) const +{ + // G4cout<<"Tabulated: Field requested at point ("<GetLatestEventNr(); + // G4int evNrKriz=6795; //457; //14250 + // if (evNr==evNrKriz) { + // std::cout<<"evNr="<0) ? 1.:-1.; + //if (evNr==evNrKriz) std::cout<<"x ="<(xdindex); + int zindex = static_cast(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="< + // // G4String processName = G4EventManager::GetEventManager()-> + // GetTrackingManager()->GetTrack()->GetStep()-> + // GetPostStepPoint()->GetProcessDefinedStep()-> + // GetProcessName(); + // The following seems to be OK: + // if (lem4Parameters::field_DecayWithSpin) { + // lem4Parameters::field_DecayWithSpin=false; + // + // } + // else { // Muon is not decaying, use a fast method of interpolation. + // Interpolate between the neighbouring points + double Bfield_R = + xField[xindex ][zindex ] * (1-xlocal) * (1-zlocal) + + xField[xindex ][zindex+1] * (1-xlocal) * zlocal + + xField[xindex+1][zindex ] * xlocal * (1-zlocal) + + xField[xindex+1][zindex+1] * xlocal * zlocal ; + Bfield[0] = (x>0) ? Bfield_R * (pppoint[0]/x) : 0.; + Bfield[1] = (x>0) ? Bfield_R * (pppoint[1]/x) : 0.; + Bfield[2] = + zField[xindex ][zindex ] * (1-xlocal) * (1-zlocal) + + zField[xindex ][zindex+1] * (1-xlocal) * zlocal + + zField[xindex+1][zindex ] * xlocal * (1-zlocal) + + zField[xindex+1][zindex+1] * xlocal * zlocal ; + + Bfield[0] = Bfield[0] * ffieldValue * z_sign; + Bfield[1] = Bfield[1] * ffieldValue * z_sign; + Bfield[2] = Bfield[2] * ffieldValue; + // } + } + + // Set some small field if field is almost zero (to avoid internal problems of Geant). + if (sqrt(Bfield[0]*Bfield[0]+Bfield[1]*Bfield[1]+Bfield[2]*Bfield[2])<0.00001*tesla) { + // if (evNr>evNrKriz) std::cout<<"bol som tuna"<-0.1*mm)&&(point[2]<0.1*mm)) { + // printf ("for point= %f %f %f B= %10.10f %10.10f %10.10f \n", + // point[0],point[1],point[2],Bfield[0]/tesla,Bfield[1]/tesla,Bfield[2]/tesla); + // } + // if (evNr>evNrKriz) std::cout<<"this is the end"< " "Reading the field grid from " << filename << " ... " << G4endl; + ifstream file( filename ); // Open the file for reading. + + // Ignore first blank line + char buffer[256]; + file.getline(buffer,256); + + // Read table dimensions + file >> nx >> ny >> nz; // Note dodgy order + + G4cout << " [ Number of values x,y,z: " + << nx << " " << ny << " " << nz << " ] " + << endl; + + // Set up storage space for table + xField.resize( nx ); + yField.resize( nx ); + zField.resize( nx ); + int ix, iy, iz; + for (ix=0; ix> xval >> yval >> zval >> bx >> by >> bz >> permeability; + if ( ix==0 && iy==0 && iz==0 ) { + minx = xval * lenUnit; + miny = yval * lenUnit; + minz = zval * lenUnit; + } + //xField[ix][iy][iz] = bx * fieldUnit * fieldValue; + //yField[ix][iy][iz] = by * fieldUnit * fieldValue; + //zField[ix][iy][iz] = bz * fieldUnit * fieldValue; + // xField[ix][iy][iz] = bx * fieldValue; + // yField[ix][iy][iz] = by * fieldValue; + // zField[ix][iy][iz] = bz * fieldValue; + xField[ix][iy][iz] = bx; + yField[ix][iy][iz] = by; + zField[ix][iy][iz] = bz; + } + } + } + file.close(); + + maxx = xval * lenUnit; + maxy = yval * lenUnit; + maxz = zval * lenUnit; + + G4cout << "\n ---> ... done reading " << endl; + + // G4cout << " Read values of field from file " << filename << endl; + G4cout << " ---> assumed the order: x, y, z, Bx, By, Bz " + << "\n ---> Min values x,y,z: " + << minx/cm << " " << miny/cm << " " << minz/cm << " cm " + << "\n ---> Max values x,y,z: " + << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm " << endl; + + + // Should really check that the limits are not the wrong way around. + if (maxx < minx) {swap(maxx,minx); invertX = true;} + if (maxy < miny) {swap(maxy,miny); invertY = true;} + if (maxz < minz) {swap(maxz,minz); invertZ = true;} + G4cout << "\nAfter reordering if neccesary" + << "\n ---> Min values x,y,z: " + << minx/cm << " " << miny/cm << " " << minz/cm << " cm " + << " \n ---> Max values x,y,z: " + << maxx/cm << " " << maxy/cm << " " << maxz/cm << " cm "; + + dx = maxx - minx; + dy = maxy - miny; + dz = maxz - minz; + G4cout << "\n ---> Dif values x,y,z (range): " + << dx/cm << " " << dy/cm << " " << dz/cm << " cm in z " + << "\n-----------------------------------------------------------" << endl; +} + +void lem4TabulatedField3D::GetFieldValue(const double point[4], + double *Bfield ) const +{ + Bfield[0]=0.; Bfield[1]=0.; Bfield[2]=0.; + + double x = point[0]+positionOffset[0]; + double y = point[1]+positionOffset[1]; + double z = point[2]+positionOffset[2]; + + + // Check that the point is within the defined region + if ( x>minx && xminy && yminz && z(xdindex); + int yindex = static_cast(ydindex); + int zindex = static_cast(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="<-1*cm)&&(point[2]<1*cm)) { + // G4cout<<"Field at point ("<GetSolid())->GetZHalfLength(); + fieldWidth = 2.*((G4Box*)lvolume->GetSolid())->GetXHalfLength(); + fieldHeight = 2.*((G4Box*)lvolume->GetSolid())->GetYHalfLength(); + + G4cout << "\n-----------------------------------------------------------" + << "\n Uniform electromagnetic field" + << "\n-----------------------------------------------------------" + << 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)); +} + + + +void lem4UniformField::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 lem4UniformField::GetNominalFieldValue() { + ///G4double EMfield[6] = {0,0,0,0,0,0}; + ///for (int i = 0; i < 6; i++){ + //return EMfield[i]; + return 0; + ///G4double nomFieldValue = (*i)->GetNominalFieldValue(); + ///EMfield[i] = Bfield[i]; + //} + ///return Bfield; ///ffieldValue; +} + +void lem4UniformField::SetNominalFieldValue(G4double newFieldValue){ ///[6]) { + ///Bfield = newFieldValue; ///ffieldValue=newFieldValue; + ///for (int i = 0; i < 6; i++){ + ///EMfield[i] = newFieldValue[i]; + ///} + G4cout<<"lem4UniformField.cc: SetNominalFieldValue method is NOT defined in this case!\n Dummy field value "<< newFieldValue << G4endl; + //G4cout<<"lem4UnifromField.cc: ffieldValue changed to="<< Bfield/tesla<<"T "< fieldLength/2.0 || std::fabs(local.x()) > fieldWidth/2.0 || std::fabs(local.y()) > fieldHeight/2.0); +} + +G4bool lem4UniformField::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); +} diff --git a/geant4/LEMuSR/src/lem4VisManager.cc b/geant4/LEMuSR/src/lem4VisManager.cc new file mode 100644 index 0000000..6f439cd --- /dev/null +++ b/geant4/LEMuSR/src/lem4VisManager.cc @@ -0,0 +1,115 @@ +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#ifdef G4VIS_USE + +#include "lem4VisManager.hh" + +// Supported drivers... + +// Not needing external packages or libraries... +#include "G4ASCIITree.hh" +#include "G4DAWNFILE.hh" +//#include "G4GAGTree.hh" // Class removed from ver. 4.9 +#include "G4HepRepFile.hh" +#include "G4HepRep.hh" +#include "G4RayTracer.hh" +#include "G4VRML1File.hh" +#include "G4VRML2File.hh" + +// Needing external packages or libraries... + +#ifdef G4VIS_USE_DAWN +#include "G4FukuiRenderer.hh" +#endif + +#ifdef G4VIS_USE_OPENGLX +#include "G4OpenGLImmediateX.hh" +#include "G4OpenGLStoredX.hh" +#endif + +#ifdef G4VIS_USE_OPENGLWIN32 +#include "G4OpenGLImmediateWin32.hh" +#include "G4OpenGLStoredWin32.hh" +#endif + +#ifdef G4VIS_USE_OPENGLXM +#include "G4OpenGLImmediateXm.hh" +#include "G4OpenGLStoredXm.hh" +#endif + +#ifdef G4VIS_USE_OIX +#include "G4OpenInventorX.hh" +#endif + +#ifdef G4VIS_USE_OIWIN32 +#include "G4OpenInventorWin32.hh" +#endif + +#ifdef G4VIS_USE_VRML +#include "G4VRML1.hh" +#include "G4VRML2.hh" +#endif + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +lem4VisManager::lem4VisManager () {} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +void lem4VisManager::RegisterGraphicsSystems () { + + // Graphics Systems not needing external packages or libraries... + RegisterGraphicsSystem (new G4ASCIITree); + RegisterGraphicsSystem (new G4DAWNFILE); +//RegisterGraphicsSystem (new G4GAGTree); // Class removed from ver. 4.9 + RegisterGraphicsSystem (new G4HepRepFile); + RegisterGraphicsSystem (new G4HepRep); + RegisterGraphicsSystem (new G4RayTracer); + RegisterGraphicsSystem (new G4VRML1File); + RegisterGraphicsSystem (new G4VRML2File); + + // Graphics systems needing external packages or libraries... + +#ifdef G4VIS_USE_DAWN + RegisterGraphicsSystem (new G4FukuiRenderer); +#endif + +#ifdef G4VIS_USE_OPENGLX + RegisterGraphicsSystem (new G4OpenGLImmediateX); + RegisterGraphicsSystem (new G4OpenGLStoredX); +#endif + +#ifdef G4VIS_USE_OPENGLWIN32 + RegisterGraphicsSystem (new G4OpenGLImmediateWin32); + RegisterGraphicsSystem (new G4OpenGLStoredWin32); +#endif + +#ifdef G4VIS_USE_OPENGLXM + RegisterGraphicsSystem (new G4OpenGLImmediateXm); + RegisterGraphicsSystem (new G4OpenGLStoredXm); +#endif + +#ifdef G4VIS_USE_OIX + RegisterGraphicsSystem (new G4OpenInventorX); +#endif + +#ifdef G4VIS_USE_OIWIN32 + RegisterGraphicsSystem (new G4OpenInventorWin32); +#endif + +#ifdef G4VIS_USE_VRML + RegisterGraphicsSystem (new G4VRML1); + RegisterGraphicsSystem (new G4VRML2); +#endif + + if (fVerbose > 0) { + G4cout << + "\nYou have successfully chosen to use the following graphics systems." + << G4endl; + PrintAvailableGraphicsSystems (); + } +} + +#endif + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... diff --git a/geant4/LEMuSR/src/meyer.cc b/geant4/LEMuSR/src/meyer.cc new file mode 100644 index 0000000..f37b639 --- /dev/null +++ b/geant4/LEMuSR/src/meyer.cc @@ -0,0 +1,1032 @@ +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$//* +// LOW ENERGY MUON SPIN RELAXATION, ROTATION, RADIATION Geant4 SIMULATION +// ID : MEYER.cc , v 1.0 +// AUTHOR: Taofiq PARAISO +// DATE : 2005-04 +// +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$// +// +// & &&&&&&&&&& &&&&&&& &&&&&&&& +// & & && && & && +// & & & & & & && +// & &&&&&&& & & &&&&&& &&&&&&&& +// & & & && & & && +// & & && & & && && & & +// &&&&&&&&&& &&&&&&&&&& & &&&&& && &&&&&&& & && +// & +// & +// & +// & +// MEYER +/* + fIRST IMPLEMENTATION BY ANLSEM,H. IN FORTRAN + C++ CONVERSION T.K.PARAISO 04-2005 + + !!! IMPORTANT !!! + + Notice: + Tables definition changes between FORTRAN and C++: + 1/ Fortran indices start at 1 and C++ indices start at 0 + 2/ Tables are defined as table[column][row] in Fortran + table[row][column] in c++ + + usefull reference + http://gershwin.ens.fr/vdaniel/Doc-Locale/Langages-Program-Scientific/Fortran/Tutorial/arrays.htm + +*/ +//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$// + + +#include "meyer.h" +#include +#include +#include +#include +#include +using namespace std; + +meyer::meyer() +{;} + +meyer::~meyer() +{;} + +void meyer::GFunctions(double* g1,double* g2, double tau) +{ + + //Diese Routine gibt in Abhaengigkeit von der reduzierten Dicke 'tau' + //Funktionswerte fuer g1 und g2 zurueck. g1 und g2 sind dabei die von + //Meyer angegebenen tabellierten Funktionen fuer die Berechnung von Halbwerts- + //breiten von Streuwinkelverteilungen. (L.Meyer, phys.stat.sol. (b) 44, 253 + //(1971)) + + + double help; + + int i; + + + double tau_[] = {0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, + 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, + 10.0, 12.0, 14.0, 16.0, 18.0, 20.0 }; + + double g1_[] = {0.050,0.115,0.183,0.245,0.305,0.363,0.419,0.473,0.525,0.575, + 0.689,0.799,0.905,1.010,1.100,1.190,1.370,1.540,1.700,1.850, + 1.990,2.270,2.540,2.800,3.050,3.290 }; + + double g2_[] = {0.00,1.25,0.91,0.79,0.73,0.69,0.65,0.63,0.61,0.59, + 0.56,0.53,0.50,0.47,0.45,0.43,0.40,0.37,0.34,0.32, + 0.30,0.26,0.22,0.18,0.15,0.13 }; + + + if (tau kann ich nicht ... => STOP"< Tabelle A + thetaSchlangeMax = 4.0; + } + else if (tau<=8.) + { + //! => Tabelle B + thetaSchlangeMax = 7.0; + } + else if (tau<=20.) + { + //! => Tabelle C + thetaSchlangeMax = 20.0; + } + else + { + std::cout<< "Subroutine ''Get_F_Function_Meyer'':"< kann ich nicht ... => STOP"<50.) + { + thetaStep = .5; + } + + else if (thetaMax>25) + { + thetaStep = .25; + } + else if (thetaMax>12.5) + { + thetaStep = .125; + } + else + { + thetaStep = .0625; + } + + + //Tabelle der F-Werte erstellen: + + nBin = 0; + std::cout<<"thetamax = "<nBinMax) + { + std::cout<< "nBin > nBinMax => EXIT"; + break; + } + + value[nBin] = sin(theta)*F; + + fValues[nBin+1] = F; // ! fuer Testzwecke + fValuesFolded[nBin+1] = sin(theta/180*M_PI)*F;// ! fuer Testzwecke + + + }// end of do loop + + + //Berechnen der Flaecheninhalte der einzelnen Kanaele sowie der Integrale: + + bigtheta:for( i = 1;i<= nBin; i++) + { + area[i] = (value[i]+value[i-1])/2.* thetaStep; + integ[i] = integ[i-1] + area[i]; + } + + + //Normiere totale Flaeche auf 1: + + rHelp = integ[nBin]; + for( i = 1; i<=nBin; i++) + { + value[i] = value[i] / rHelp; + area[i] = area[i] / rHelp; + integ[i] = integ[i] / rHelp; + } + + + //vorerst noch: gib Tabelle in Datei und Histogrammfile aus: + + //! Berechne die Werte fuer theta=0: + + F_Functions_Meyer(tau,0.,&f1,&f2); + F = Meyer_faktor4*Meyer_faktor4 * Ekin*Ekin /2 /M_PI * (f1 - Meyer_faktor3*f2);// TAO, Anselm was: Meyer_faktor5 * Ekin*Ekin * (f1 - Meyer_faktor3*f2); + fValues[1] = F; + fValuesFolded[1] = 0.; + + //! Gib die Werte in das Tabellenfile aus: + + /* ofstream Mprint("testmeyer.out"); + theta = thetaStep; + if (!Mprint.is_open()) exit(8); + for( i = 1; i<=nBin+1;i++) + { + Mprint << theta<< " "<< fValues[i]/fValues[1]<<" " << fValuesFolded[i]< Reihe mit hoeherem Index + //iColumn = 2 => Reihe mit kleinerem Index + + + iColumn = 1; + + // 5 continue; + do{ + + if (column_<=8) + { + //! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + //! Werte aus 1. Tabelle: 0.2 <= tau <= 1.8 + + column = column_; + // std::cout<<"thetaSchlange = "< * lem4EventAction::RandomNrInitialisers=NULL; + +//long lem4EventAction::myEventNr=0; + +lem4EventAction::lem4EventAction() { + pointer=this; + fieldValueStart=0; + pointerToSeedVector = new vector; + 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"<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:"< 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 "<SetFieldValue(actualFieldValue); + G4cout<<"Event "<SetFieldValue(actualFieldValue); + G4cout<<"Event "<SetFieldValue(actualFieldValue); + + // if (lem4DetectorMessenger::setRandomNrSeedAccordingEventNr) { + if (setRandomNrSeedFromFile) { + // G4cout<<"RandomNrInitialisers.size()="<size()<size()) { + G4cout <<"lem4EventAction.cc: seed will be set to="<< pointerToSeedVector->at(thisEventNr)<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="<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; iGetTrajectoryContainer()))[i]); + trj->DrawTrajectory(1000); + } + } +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +vector * lem4EventAction::pointerToSeedVector=NULL; +vector * 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<<"---&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&---"< +#include +#include +#include + +/* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + 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 < ABORTED!" < 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; +}