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

11.20.2009. - Kamil Sedlak

Browse Source 
This svn update includes several smaller corrections and updates accumulated
over last few months:
1) For the Geant4.9.2 it was necessary to remove the privately modified file
   src/G4EqEMFieldWithSpin.cc.  Our corrections in this file (and also in
   the file src/G4DecayWithSpin.cc) were adopted by the Geant developers into
   the official Geant code, and therefore these two files were deleted.
   However, if one uses older version of Geant (i.e. Geant4.9.1 or older),
   one should rename the G4EqEMFieldWithSpin.cc_for_Geant4.9.1_and_older in
   his/her src directory to G4EqEMFieldWithSpin.cc to correct for a Geant
   bug.
2) Implementation of save_polx, save_poly and save_polz variables
3) Implementation of the field map normalisation within the field map
   itself (was already possible for the field-map formats generated by Toni,
   now it is extended also for the field maps generated by OPERA).
4) Possibility to swap and invert x and y axis read out from the TURTLE
   file.
5) Perhaps some other tiny changes
This commit is contained in:
sedlak
2009-11-20 10:29:02 +00:00
parent d05e77dbe5
commit b5594b4513
15 changed files with 183 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

0
src/G4DecayWithSpin.cc → src/G4DecayWithSpin.cc_for_Geant4.9.1_and_older
View File

0
src/G4EqEMFieldWithSpin.cc → src/G4EqEMFieldWithSpin.cc_for_Geant4.9.1_and_older
View File

9
src/musrDetectorConstruction.cc
View File

@@ -462,6 +462,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if ((strcmp(name,"Target")==0)||(strcmp(name,"target")==0)) {
musrRootOutput::store_muTargetTime = true; musrRootOutput::store_muTargetPolX = true;
musrRootOutput::store_muTargetPolY = true; musrRootOutput::store_muTargetPolZ = true;
musrRootOutput::store_muTargetMomX = true;
musrRootOutput::store_muTargetMomY = true; musrRootOutput::store_muTargetMomZ = true;
}
if ((strcmp(name,"M0")==0)||(strcmp(name,"m0")==0)) {
musrRootOutput::store_muM0Time = true; musrRootOutput::store_muM0PolX = true;
@@ -577,7 +579,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(typeOfField,"quadrupole")==0) {
float halfLength, fieldRadius, gradientValue, gradientValueFinal, fringeFactor;
int gradientNrOfSteps;
int gradientNrOfSteps = 0;
char logicalVolumeName[100];
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %g %g %g %s %g %g %d",&halfLength,&fieldRadius,&fringeFactor,logicalVolumeName,
&gradientValue,&gradientValueFinal,&gradientNrOfSteps);
@@ -676,7 +678,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
G4LogicalVolume* pLogVol = FindLogicalVolume(tmpString2);
if (pLogVol==NULL) {
G4cout << "ERROR! musrDetectorConstruction::Construct(): SetUserLimits: Logical Volume \""
<< tmpString3 <<"\" not found!"<<G4endl<<"S T O P F O R C E D"<<G4endl;
<< tmpString2 <<"\" not found!"<<G4endl<<"S T O P F O R C E D"<<G4endl;
ReportGeometryProblem(line);
}
// float maxStep, ;
@@ -871,6 +873,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"muTargetPolX")==0) {musrRootOutput::store_muTargetPolX = false;}
if (strcmp(tmpString2,"muTargetPolY")==0) {musrRootOutput::store_muTargetPolY = false;}
if (strcmp(tmpString2,"muTargetPolZ")==0) {musrRootOutput::store_muTargetPolZ = false;}
if (strcmp(tmpString2,"muTargetMomX")==0) {musrRootOutput::store_muTargetMomX = false;}
if (strcmp(tmpString2,"muTargetMomY")==0) {musrRootOutput::store_muTargetMomY = false;}
if (strcmp(tmpString2,"muTargetMomZ")==0) {musrRootOutput::store_muTargetMomZ = false;}
if (strcmp(tmpString2,"muM0Time")==0) {musrRootOutput::store_muM0Time = false;}
if (strcmp(tmpString2,"muM0PolX")==0) {musrRootOutput::store_muM0PolX = false;}
if (strcmp(tmpString2,"muM0PolY")==0) {musrRootOutput::store_muM0PolY = false;}

36
src/musrPhysicsList.cc
View File

@@ -51,6 +51,9 @@
musrPhysicsList::musrPhysicsList(): G4VUserPhysicsList()
{
defaultCutValue = 0.1*mm;
cutForGamma = 0.1*mm;
cutForElectron = 0.1*mm;
cutForMuon = 0.01*mm;
SetVerboseLevel(0);
}
@@ -248,6 +251,7 @@ void musrPhysicsList::ConstructEM()
if ( (strcmp(tmpString0,"/musr/ignore")!=0)&&(strcmp(tmpString0,"/musr/command")!=0) ) continue;
if (strcmp(tmpString1,"process")!=0) continue;
float tmpCutValue;
if ((strcmp(tmpString2,"addProcess")==0)||(strcmp(tmpString2,"addDiscreteProcess")==0)) {
char charParticleName[100], charProcessName[100];
sscanf(&line[0],"%*s %*s %s %s %s",tmpString2,charParticleName,charProcessName);
@@ -346,6 +350,23 @@ void musrPhysicsList::ConstructEM()
}
}
else if (strcmp(tmpString2,"cutForGamma")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
cutForGamma = tmpCutValue;
}
else if (strcmp(tmpString2,"cutForElectron")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
cutForElectron = tmpCutValue;
}
//cks UNFORTUNATELY cutForMuon does not work!!!
// else if (strcmp(tmpString2,"cutForMuon")==0) {
// sscanf(&line[0],"%*s %*s %*s %g",&tmpCutValue);
// cutForMuon = tmpCutValue;
// }
else ReportProblemWithProcessDefinition(line);
}
}
@@ -440,7 +461,20 @@ void musrPhysicsList::SetCuts()
//G4VUserPhysicsList::SetCutsWithDefault method sets
//the default cut value for all particle types
//
SetCutsWithDefault();
//cks 6.10.2009 SetCutsWithDefault();
// set cut values for gamma at first and for e- second and next for e+,
// because some processes for e+/e- need cut values for gamma
SetCutValue(cutForGamma, "gamma");
SetCutValue(cutForElectron, "e-");
SetCutValue(cutForElectron, "e+");
//cks - This command does not work for muons: SetCutValue(cutForMuon, "mu-");
//cks - This command does not work for muons: SetCutValue(cutForMuon, "mu+");
// G4cout<<"Kamil: cutForGamma = "<<cutForGamma/mm<<" mm"<<G4endl;
// G4cout<<"Kamil: cutForElectron = "<<cutForElectron/mm<<" mm"<<G4endl;
// G4cout<<"Kamil: cutForMuons = "<<cutForMuon/mm<<" mm"<<G4endl;
if (verboseLevel>0) DumpCutValuesTable();
DumpCutValuesTable();

37
src/musrPrimaryGeneratorAction.cc
View File

@@ -82,6 +82,7 @@ musrPrimaryGeneratorAction::musrPrimaryGeneratorAction(
// cks Implement also alpha and proton particles for the simulation of Juan Pablo Urrego
alphaParticle= particleTable->FindParticle("alpha");
protonParticle= particleTable->FindParticle("proton");
turtleInterpretAxes="undefined";
// csk
G4int n_particle = 1;
@@ -149,6 +150,7 @@ void musrPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
}
numberOfGeneratedEvents++;
sscanf(&line[0],"%g %g %g %g %g %g %g %d %d",&xTmp,&xAngleTmp,&yTmp,&yAngleTmp,&pTmp,&dummy1,&dummy2,&Ztmp,&Atmp);
if (turtleInterpretAxes!="undefined") swapTheAxisInTurtle(xTmp,xAngleTmp,yTmp,yAngleTmp);
if (boolPrintInfoAboutGeneratedParticles) {
G4cout<<"musrPrimaryGeneratorAction::GeneratePrimaries: Turtle input for this event: "
<<xTmp<<", "<<xAngleTmp<<" "<<yTmp<<" "<<yAngleTmp<<" "<< pTmp<<G4endl;
@@ -475,3 +477,38 @@ void musrPrimaryGeneratorAction::SetPrimaryParticule(G4String particleName) {
exit(1);
}
}
//===============================================================================
void musrPrimaryGeneratorAction::swapTheAxisInTurtle(float& x_x, float& x_xprime ,float& y_y, float& y_yprime) {
// G4cout<<"turtleInterpretAxes = "<<turtleInterpretAxes<<G4endl;
if (turtleInterpretAxes=="-xy") {
x_x = -x_x; x_xprime = -x_xprime;
}
else if (turtleInterpretAxes=="x-y") {
y_y = -y_y; y_yprime = -y_yprime;
}
else if (turtleInterpretAxes=="-x-y") {
x_x = -x_x; x_xprime = -x_xprime; y_y = -y_y; y_yprime = -y_yprime;
}
else if (turtleInterpretAxes=="yx") {
float tmpX = x_x; float tmpXprime = x_xprime;
x_x = y_y; x_xprime = y_yprime; y_y = tmpX; y_yprime = tmpXprime;
}
else if (turtleInterpretAxes=="-yx") {
float tmpX = x_x; float tmpXprime = x_xprime;
x_x = y_y; x_xprime = y_yprime; y_y = -tmpX; y_yprime = -tmpXprime;
}
else if (turtleInterpretAxes=="y-x") {
float tmpX = x_x; float tmpXprime = x_xprime;
x_x = -y_y; x_xprime = -y_yprime; y_y = tmpX; y_yprime = tmpXprime;
}
else if (turtleInterpretAxes=="-y-x") {
float tmpX = x_x; float tmpXprime = x_xprime;
x_x = -y_y; x_xprime = -y_yprime; y_y = -tmpX; y_yprime = -tmpXprime;
}
else {
G4cout<<"musrPrimaryGeneratorAction::swapTheAxisInTurtle: Not known how to inpterpret turtleInterpretAxes="<<turtleInterpretAxes<<G4endl;
G4cout<<"S T O P F O R C E D" << G4endl;
exit(1);
}
}

10
src/musrPrimaryGeneratorMessenger.cc
View File

@@ -130,6 +130,13 @@ musrPrimaryGeneratorMessenger::musrPrimaryGeneratorMessenger(musrPrimaryGenerato
setTurtleZ0Cmd->SetParameterName("mes_z0Turtle",true);
setTurtleZ0Cmd->SetDefaultUnit("mm");
setTurtleInterpretAxesCmd = new G4UIcmdWithAString("/gun/turtleInterpretAxes",this);
setTurtleInterpretAxesCmd->SetGuidance("Specify whether the x and y axes of the position in TURTLE should be interpretted differently.");
setTurtleInterpretAxesCmd->SetGuidance("The following options are supported: x-y, -xy, -x-y, yx, y-x, -yx, -y-x .");
setTurtleInterpretAxesCmd->SetGuidance("Example: the option y-x means that first four coordinates in the TURTLE input file");
setTurtleInterpretAxesCmd->SetGuidance(" are interpreded as y, yprime, -x, -xprime.");
setTurtleInterpretAxesCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
setTurtleMomentumBite = new G4UIcmdWith3Vector("/gun/turtleMomentumBite",this);
setTurtleMomentumBite->SetGuidance(" Modify smearing of the turtle momentum bite. The first value is the mean momentum in MeV/c,");
setTurtleMomentumBite->SetGuidance(" the second value is the smearing factor in per cent, by which the momentum bite,");
@@ -168,6 +175,7 @@ musrPrimaryGeneratorMessenger::~musrPrimaryGeneratorMessenger()
delete setMuonDecayTimeCmd;
delete setTurtleCmd;
delete setTurtleZ0Cmd;
delete setTurtleInterpretAxesCmd;
delete setTurtleMomentumBite;
delete setTurtleEventNrCmd;
}
@@ -212,6 +220,8 @@ void musrPrimaryGeneratorMessenger::SetNewValue(G4UIcommand * command,G4String n
{ musrAction->SetTurtleInput(newValue); }
if( command == setTurtleZ0Cmd)
{ musrAction->SetTurtleZ0(setTurtleZ0Cmd->GetNewDoubleValue(newValue)); }
if( command == setTurtleInterpretAxesCmd)
{ musrAction->SetTurtleInterpretAxes(newValue);}
if( command == setTurtleMomentumBite)
{ musrAction-> SetTurtleMomentumBite(setTurtleMomentumBite->GetNew3VectorValue(newValue)); }
if( command == setTurtleEventNrCmd)

16
src/musrRootOutput.cc
View File

@@ -87,6 +87,9 @@ G4bool musrRootOutput::store_muTargetTime = false;
G4bool musrRootOutput::store_muTargetPolX = false;
G4bool musrRootOutput::store_muTargetPolY = false;
G4bool musrRootOutput::store_muTargetPolZ = false;
G4bool musrRootOutput::store_muTargetMomX = false;
G4bool musrRootOutput::store_muTargetMomY = false;
G4bool musrRootOutput::store_muTargetMomZ = false;
G4bool musrRootOutput::store_muM0Time = false;
G4bool musrRootOutput::store_muM0PolX = false;
G4bool musrRootOutput::store_muM0PolY = false;
@@ -177,6 +180,9 @@ void musrRootOutput::BeginOfRunAction() {
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_muTargetMomX) {rootTree->Branch("muTargetMomX",&muTargetMomX_t,"muTargetMomX/D");}
if (store_muTargetMomY) {rootTree->Branch("muTargetMomY",&muTargetMomY_t,"muTargetMomY/D");}
if (store_muTargetMomZ) {rootTree->Branch("muTargetMomZ",&muTargetMomZ_t,"muTargetMomZ/D");}
if (store_muM0Time) {rootTree->Branch("muM0Time",&muM0Time_t,"muM0Time/D");}
if (store_muM0PolX) {rootTree->Branch("muM0PolX",&muM0PolX_t,"muM0PolX/D");}
if (store_muM0PolY) {rootTree->Branch("muM0PolY",&muM0PolY_t,"muM0PolY/D");}
@@ -248,6 +254,9 @@ void musrRootOutput::BeginOfRunAction() {
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");
rootTree->Branch("save_polx",&save_polx,"save_polx[save_n]/D");
rootTree->Branch("save_poly",&save_poly,"save_poly[save_n]/D");
rootTree->Branch("save_polz",&save_polz,"save_polz[save_n]/D");
}
// htest1 = new TH1F("htest1","The debugging histogram 1",50,-4.,4.);
@@ -313,6 +322,7 @@ void musrRootOutput::ClearAllRootVariables() {
muDecayDetID_t=-1000;
muDecayPolX_t=-1000; muDecayPolY_t=-1000; muDecayPolZ_t=-1000;
muTargetTime_t=-1000; muTargetPolX_t=-1000; muTargetPolY_t=-1000; muTargetPolZ_t=-1000;
muTargetMomX_t=-1000; muTargetMomY_t=-1000; muTargetMomZ_t=-1000;
muM0Time_t=-1000; muM0PolX_t=-1000; muM0PolY_t=-1000; muM0PolZ_t=-1000;
muM1Time_t=-1000; muM1PolX_t=-1000; muM1PolY_t=-1000; muM1PolZ_t=-1000;
muM2Time_t=-1000; muM2PolX_t=-1000; muM2PolY_t=-1000; muM2PolZ_t=-1000;
@@ -371,7 +381,7 @@ G4int musrRootOutput::ConvertProcessToID(std::string processName) {
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
void musrRootOutput::SetSaveDetectorInfo (G4int ID, G4int particleID, G4double ke,
G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz) {
G4double x, G4double y, G4double z, G4double px, G4double py, G4double pz, G4double polx, G4double poly, G4double polz) {
if (save_n>=save_nMax) {
char message[200];
sprintf(message,"musrRootOutput.cc::SetSaveDetectorInfo(): more \"save\" hits then allowed: save_nMax=%i",save_nMax);
@@ -387,7 +397,9 @@ void musrRootOutput::SetSaveDetectorInfo (G4int ID, G4int particleID, G4double k
save_px[save_n]=px/MeV;
save_py[save_n]=py/MeV;
save_pz[save_n]=pz/MeV;
save_polx[save_n]=polx;
save_poly[save_n]=poly;
save_polz[save_n]=polz;
save_n++;
}
}

9
src/musrSteppingAction.cc
View File

@@ -208,7 +208,10 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
G4double px_save=aStep->GetPreStepPoint()->GetMomentum().x();
G4double py_save=aStep->GetPreStepPoint()->GetMomentum().y();
G4double pz_save=aStep->GetPreStepPoint()->GetMomentum().z();
myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,px_save,py_save,pz_save);
G4double polx_save=aStep->GetPreStepPoint()->GetPolarization().x();
G4double poly_save=aStep->GetPreStepPoint()->GetPolarization().y();
G4double polz_save=aStep->GetPreStepPoint()->GetPolarization().z();
myRootOutput->SetSaveDetectorInfo(tmpVolumeID,particle_id_save,ke_save,x_save,y_save,z_save,px_save,py_save,pz_save,polx_save,poly_save,polz_save);
}
}
}
@@ -222,6 +225,10 @@ void musrSteppingAction::UserSteppingAction(const G4Step* aStep) {
muAlreadyWasInTargetInThisEvent=true;
myRootOutput->SetPolInTarget(aTrack->GetPolarization());
myRootOutput->SetTimeInTarget(aTrack->GetGlobalTime());
G4ThreeVector p_beforeEnteringTarget = aTrack->GetMomentum() - aStep->GetDeltaMomentum(); // must be - sign because DeltaMomentum is negative
// std::cout<<"(aStep->GetDeltaMomentum()).z() = "<<(aStep->GetDeltaMomentum()).z()<<std::endl;
myRootOutput->SetMomentumInTarget(p_beforeEnteringTarget);
// myRootOutput->SetMomentumInTarget(aTrack->GetMomentum());
}
}
else if ((actualVolume=="log_M0")||(actualVolume=="log_m0")) {

17
src/musrTabulatedElementField.cc
View File

@@ -66,6 +66,8 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
}
char buffer[256];
char tmpString1[100]="Unset";
float fvalue;
G4bool boolMinimaAndMaximaDefinedInTheFile = false;
if (fldType=='E') G4cout<<" Electric field ";
if (fldType=='B') G4cout<<" Magnetic field ";
@@ -118,6 +120,11 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
// be the last line of the header.
do {
file.getline(buffer,256);
sscanf(&buffer[0],"%s %g",tmpString1,&fvalue);
if (strcmp(tmpString1,"fieldNormalisation")==0) {
fieldNormalisation = fvalue;
G4cout << "DEBUG: musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
}
} while ( buffer[1]!='0');
}
else if ((strcmp(fieldTableType,"2DE")==0)||(strcmp(fieldTableType,"2DB")==0)||(strcmp(fieldTableType,"2DEf")==0)) {
@@ -142,6 +149,11 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
file >> nx >> nz >> nDummy;
do {
file.getline(buffer,256);
sscanf(&buffer[0],"%s %g",tmpString1,&fvalue);
if (strcmp(tmpString1,"fieldNormalisation")==0) {
fieldNormalisation = fvalue;
G4cout << "DEBUG: musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
}
} while ( buffer[1]!='0');
}
else if ((strcmp(fieldTableType,"2D")==0)||(strcmp(fieldTableType,"2DBOpera")==0)) {
@@ -154,6 +166,11 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
// G4cout << nx <<" "<< nDummy <<" "<< nz<<G4endl;
do {
file.getline(buffer,256);
sscanf(&buffer[0],"%s %g",tmpString1,&fvalue);
if (strcmp(tmpString1,"fieldNormalisation")==0) {
fieldNormalisation = fvalue;
G4cout << "DEBUG: musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
}
} while ( buffer[1]!='0');
}
else {