From 29e49736ea5d43a6b29458b096c7a55a0a036051 Mon Sep 17 00:00:00 2001
From: Kamil Sedlak <kamil.sedlak@psi.ch>
Date: Mon, 15 Mar 2010 13:19:31 +0000
Subject: [PATCH] 15.3.2010  Kamil Sedlak   1) Implemented a posibility to add
 "Models" to physics processes      (this should be a first step towards the
 implementation of the      Meyer scattering.  The documentation has not been
 updated yet.   2) Minor change - implementation of a special volume for GPD

---
 src/musrDetectorConstruction.cc | 42 ++++++++++++++++++++++++++-
 src/musrPhysicsList.cc          | 50 +++++++++++++++++++++++++++++++--
 2 files changed, 89 insertions(+), 3 deletions(-)

diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc
index d4bc3f2..722f1ef 100644
--- a/src/musrDetectorConstruction.cc
+++ b/src/musrDetectorConstruction.cc
@@ -354,9 +354,49 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct()  {
 		  <<"  orientation="<<orientation<<G4endl;
 	    ReportGeometryProblem(line);
 	  }
-
 	  solid =  new G4IntersectionSolid(solidName, solidDetTube, solidDetBox, yRot, zTrans);
 	}
+	else if (strcmp(tmpString2,"GPDsampleHolderA")==0){
+	  // First part of the GPD sample holder, where posx, posy, posz = centre of the whole (long) tube 
+	  //                                                      (=111.25mm below the centre of the holes)
+	  // 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 %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;
+	  //	  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("SolidOuterDetTube",0.*mm,x2*mm,x3*mm,x4*deg,x5*deg);
+	  //	  solid =  new G4SubtractionSolid(solidName, solidOuterDetTube, solidInnerDetBox);
+
+
+	  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",sensitiveDet,&volumeID);
+	  solidName+=name;
+	  G4double upperPartHalfHeight = 65*mm;
+	  G4double lowerPartHalfHeight = (287.5-65)/2*mm;
+	  G4double halfHeight = upperPartHalfHeight+lowerPartHalfHeight;
+	  G4double innerR = 37.*mm; G4double outerR = 37.5*mm;
+	  
+	  //	  G4Box* solidDetBox = new G4Box("SolidDetBox",x2*mm,x2*mm,x3*mm+roundingErr);
+	  G4Tubs* solidGPDTubeA1 = new G4Tubs("SolidGPDTubeA1",innerR,outerR,halfHeight,0,360);
+	  G4Box*  solidGPDBoxA2  = new G4Box ("SolidGPDBoxA2",10*mm,38*mm, 70*mm);
+	  G4Box*  solidGPDBoxA5  = new G4Box ("solidGPDBoxA5",5*mm,40*mm,30*mm); 
+
+	  G4RotationMatrix* yRotA12 = new G4RotationMatrix();
+	  G4ThreeVector zTransA12( 30*mm,0, (111.25+20)*mm);
+	  G4ThreeVector zTransA13(-30*mm,0, (111.25+20)*mm);
+	  G4SubtractionSolid* solidA12 =  new G4SubtractionSolid("solidA12", solidGPDTubeA1, solidGPDBoxA2, yRotA12, zTransA12);
+	  G4SubtractionSolid* solidA123 =  new G4SubtractionSolid("solidA123", solidA12 , solidGPDBoxA2, yRotA12, zTransA13);
+	  G4ThreeVector zTransA5(0,0,111.25*mm);
+	  solid =  new G4SubtractionSolid(solidName, solidA123, solidGPDBoxA5, yRotA12, zTransA5);
+	  //G4cout<<"solidName="<<solidName<<"  solidA123="<<solidA123<<"   solidGPDBoxA5="<<solidGPDBoxA5<<"   yRotA12="<<yRotA12<<"  zTransA5="<<zTransA5<<G4endl;
+	  //solid =  new G4SubtractionSolid(solidName, solidGPDTubeA1, solidGPDBoxA2, yRotA12, zTransA12);
+	  
+	  G4cout<<"Debug 80"<<G4endl;
+	}
+
 	else ReportGeometryProblem(line);
 
 	G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm);
diff --git a/src/musrPhysicsList.cc b/src/musrPhysicsList.cc
index d634e4b..bb31700 100644
--- a/src/musrPhysicsList.cc
+++ b/src/musrPhysicsList.cc
@@ -234,6 +234,12 @@ void musrPhysicsList::ConstructProcess()
 // For a simple Muonium "scattering" when Mu hits solid materials
 #include "musrMuScatter.hh"
 
+// For models
+#include "G4ProcessTable.hh"
+#include "G4WentzelVIModel.hh"
+#include "G4UrbanMscModel90.hh"
+#include "G4UrbanMscModel92.hh"
+#include "G4UrbanMscModel93.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
@@ -258,7 +264,7 @@ void musrPhysicsList::ConstructEM()
       if (strcmp(tmpString1,"process")!=0)  continue; 
       
       float tmpCutValue;
-      if ((strcmp(tmpString2,"addProcess")==0)||(strcmp(tmpString2,"addDiscreteProcess")==0)) {
+      if ((strcmp(tmpString2,"addProcess")==0)||(strcmp(tmpString2,"addDiscreteProcess")==0)||(strcmp(tmpString2,"addModel")==0)) {
 	char charParticleName[100], charProcessName[100];
 	sscanf(&line[0],"%*s %*s %s %s %s",tmpString2,charParticleName,charProcessName);
 	G4cout<<"musrPhysicsList: Defining process "<<charProcessName<<" for "<<charParticleName<<G4endl;
@@ -267,7 +273,7 @@ void musrPhysicsList::ConstructEM()
 	G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable() -> FindParticle(stringParticleName);
 	//      G4cout<<"particleDefinition of "<<stringParticleName<<" = "<<particleDefinition<<G4endl;
 	if (particleDefinition==NULL) {
-	  sprintf(eMessage,"musrPhysicsList:  Partile \"%s\" not found in G4ParticleTable when trying to assign process  \"%s\".",
+	  sprintf(eMessage,"musrPhysicsList:  Partile \"%s\" not found in G4ParticleTable when trying to find or assign process  \"%s\".",
 		  charParticleName,charProcessName);
 	  musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
 	}
@@ -357,6 +363,46 @@ void musrPhysicsList::ConstructEM()
 	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
 	  }
 	}
+	else if (strcmp(tmpString2,"addModel")==0) {
+	  G4ProcessTable* processTable = G4ProcessTable::GetProcessTable();
+	  //	  G4VProcess* myProc =  processTable->FindProcess(charProcessName,particleDefinition);
+	  // pManager
+	  
+	  G4int modelPriority;
+	  char charModelName[100]="";
+	  sscanf(&line[0],"%*s %*s %*s %*s %*s %s %d",charModelName,&modelPriority);
+	  G4String stringModelName = charModelName;
+	  G4cout<<" Adding model "<<charModelName<<" with priority "<<modelPriority<<" to process "<<charProcessName<<"."<<G4endl;
+	  //	  if (myProc==NULL) {
+	  //	    sprintf(eMessage,"musrPhysicsList:  Trying to add model \"%s\" , but the process \"%s\" not found in G4ProcessTable.",
+	  //		    charModelName,charProcessName);
+	  //	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
+	  //	  }
+	  if      ((stringModelName=="G4WentzelVIModel")&&(stringProcessName=="G4MuMultipleScattering")) {
+	    G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
+	    mmm->AddEmModel(modelPriority, new G4WentzelVIModel());
+	  }
+	  else if ((stringModelName=="G4UrbanMscModel90")&&(stringProcessName=="G4MuMultipleScattering")) {
+	    G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
+	    mmm->AddEmModel(modelPriority, new G4UrbanMscModel90());
+	  }
+	  else if ((stringModelName=="G4UrbanMscModel92")&&(stringProcessName=="G4MuMultipleScattering")) {
+	    G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
+	    mmm->AddEmModel(modelPriority, new G4UrbanMscModel92());
+	  }
+	  else if ((stringModelName=="G4UrbanMscModel93")&&(stringProcessName=="G4MuMultipleScattering")) {
+	    G4MuMultipleScattering* mmm = (G4MuMultipleScattering*) processTable->FindProcess("muMsc",particleDefinition);
+	    mmm->AddEmModel(modelPriority, new G4UrbanMscModel93());
+	  }
+
+	  else {
+	    sprintf(eMessage,"musrPhysicsList:  Model \"%s\" is not implemented for \"%s\" in musrPhysicsList.cc for addModel.  It can be easily added.",
+		    charModelName,charProcessName);
+	    musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
+	  }
+
+	}
+
       }
 
       else if (strcmp(tmpString2,"cutForGamma")==0) {