#include "sr1DetectorConstruction.hh"
#include "sr1DetectorMessenger.hh"
#include "sr1MagneticField.hh"
#include "sr1TabulatedField3D.hh"
#include "sr1TabulatedField2D.hh"
#include "sr1TabulatedElementField2Df.hh" // Read a 2D axial field map folded along r-axis (only positive z)
#include "sr1TabulatedElementField2D.hh"  // Read a 2D axial field map
#include "sr1TabulatedElementField3D.hh"  // Read a 3D field map
///#include "sr1Axial2DElField.hh"        // ADDED to account for Gumplinger field method
#include "sr1UniformField.hh"             // Uniform EM field (B,E)


#include "sr1GaussianField.hh"

//#include "G4EqMagElectricField.hh" //Added by TS to account for electrical fields
//#include "G4UniformElectricField.hh"
//#include "G4ElectricField.hh"



#include "sr1ScintSD.hh"
#include "sr1EventAction.hh"

#include "G4GeometryManager.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4LogicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4SDManager.hh"
#include "G4PVParameterised.hh"
#include "G4UnionSolid.hh"
#include "G4IntersectionSolid.hh"
#include "G4SubtractionSolid.hh"
//#include "G4BooleanSolid.hh"

#include "G4Material.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Cons.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4Tubs.hh"
#include "G4Sphere.hh"
#include "G4Trd.hh"
#include "G4PVDivision.hh"

#include "G4ChordFinder.hh"
#include "G4Mag_SpinEqRhs.hh"
#include "G4ClassicalRK4.hh"
#include "G4MagIntegratorStepper.hh"
#include "G4TransportationManager.hh"
//#include "G4FieldManager.hh"
//#include "G4FieldTrack.hh"   //cks for visualisation
#include "G4PropagatorInField.hh" //cks for visualisation
#include "G4UserLimits.hh"

// In order to implement overlaping field as was done by 
// Gumplinger in examples/extended/field/field04/
#include "F04GlobalField.hh"

#include "G4VisAttributes.hh"
#include "G4Colour.hh"
#include "G4ios.hh"

#include "sr1RootOutput.hh"   //cks for storing some info in the Root output file
#include "sr1Parameters.hh"
#include "sr1ErrorMessage.hh"
#include "sr1SteppingAction.hh" // CHANGED. TS: What is its purpose?

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

sr1DetectorConstruction::sr1DetectorConstruction()
:parameterFileName("Unset"), checkOverlap(true), aScintSD(0)
{  
  DefineMaterials();
  detectorMessenger = new sr1DetectorMessenger(this);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
sr1DetectorConstruction::~sr1DetectorConstruction()
{
  delete detectorMessenger;             
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

#include "G4LogicalVolumeStore.hh"
G4VPhysicalVolume* sr1DetectorConstruction::Construct()  {
  // Clean old geometry, if any
  //
  G4GeometryManager::GetInstance()->OpenGeometry();
  G4PhysicalVolumeStore::GetInstance()->Clean();
  G4LogicalVolumeStore::GetInstance()->Clean();
  G4SolidStore::GetInstance()->Clean();

  ///G4cout<< "Parameter file name = "<<parameterFileName<<G4endl;
  char howIsTheFieldDefined[100] = "NOT_DEFINED_YET";     // can be set to "GLOBAL_FIELD" or "DIFFERENT_FIELDS"

//  G4double roundingErr=0.001*mm;     
//  G4double roundingErr=0.01*mm;        // 0.01mm precision is OK for displaying subtracted volumes, while 0.001mm is not
//----------------------

  sr1RootOutput* myRootOutput = sr1RootOutput::GetRootInstance();

  G4VPhysicalVolume* pointerToWorldVolume=NULL;
  //  Read detector configuration parameters from the steering file:
  FILE *fSteeringFile=fopen(parameterFileName.c_str(),"r");
  if (fSteeringFile==NULL) {
    if (parameterFileName=="Unset") {
      G4cout<<"sr1DetectorConstruction:  No input macro file specified"<<G4endl;
    }
    else {
      G4cout << "E R R O R :    Failed to open detector configuration file " << parameterFileName << G4endl;
    }
    G4cout << "S T O P    F O R C E D" << G4endl;
    exit(1);
  }
  G4cout << "Detector configuration file \"" << parameterFileName << "\" was opened."<<G4endl;
  char  line[501];

    
  // Write out the macro configuration file to the standard output (file/screen):
  G4bool PrintMacro = true;
  if (PrintMacro) {
    G4cout << "Content of the configuration file used for this run:"<<G4endl<<G4endl; 
      while (!feof(fSteeringFile)) {
        fgets(line,500,fSteeringFile);
	  if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) { // Do not print comments
              G4cout << line;
	      }
      }
    G4cout << G4endl<<"End of the configuration file.\n"<<G4endl;
    }
 
  
  // Main Loop
  rewind (fSteeringFile);
  char eMessage[200];
  while (!feof(fSteeringFile)) {
    fgets(line,500,fSteeringFile);
    //    if (line[0]!='*') {
    if ((line[0]!='#')&&(line[0]!='\n')&&(line[0]!='\r')) {
      char tmpString0[100]="Unset";
      sscanf(&line[0],"%s",tmpString0);
      if (strcmp(tmpString0,"/sr1/command")!=0) { break; }
      

      char tmpString1[100]="Unset",tmpString2[100]="Unset",tmpString3[100]="Unset";
      sscanf(&line[0],"%*s %s %s",tmpString1,tmpString2);

      // Define rotation matrix, which might be later on used for some volumes
      if (strcmp(tmpString1,"rotation")==0){
	char matrixName[100]="Unset";
	float par1,par2,par3,par4=0	;
	float phi,theta,psi;
	float vx,vy,vz,delta;
	sscanf(&line[0],"%*s %*s %s %g %g %g %g",matrixName,&par1,&par2,&par3,&par4);
//	sscanf(&line[0],"%*s %*s %s %g %g %g",matrixName,&phi,&theta,&psi);
	if (par4==0) {
	  phi=par1; theta=par2; psi=par3;
	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(phi*deg,theta*deg,psi*deg);
	  pointerToRotationMatrix[matrixName]=pRotMatrix;
	}
	else {
	  vx=par1; vy=par2; vz=par3; delta=par4;
	  G4RotationMatrix* pRotMatrix = new G4RotationMatrix(G4ThreeVector(vx,vy,vz),delta*deg);
	  pointerToRotationMatrix[matrixName]=pRotMatrix;
	}
//	pointerToRotationMatrix[matrixName]=pRotMatrix;
	//	G4cout<<"matrixName="<<matrixName<<"  Rot. matrix. nr.="<<pointerToRotationMatrix.size()
	//	      <<"   pointer="<<pRotMatrix<<"      pointerToRotationMatrix="<<pointerToRotationMatrix[matrixName]<<G4endl;
      // Added by TS for rotation not according to Euler, but to local axes!
      }     
      
      
      // Define magnetic or electric field, which will be later on used for a volume
      else if (strcmp(tmpString1,"field")==0){
        char typeOfField[100]="Unset";
        char nameOfField[100]="Unset";
        float fieldValue=0.000000001;

        if (strcmp(tmpString2,"magnetic")==0){
          G4FieldManager* pFieldMgr = new G4FieldManager();
	  //          sr1MagneticField*  sr1Field = NULL;
	  G4MagneticField* sr1Field = NULL;

	  sscanf(&line[0],"%*s %*s %*s %s",typeOfField);
          if (strcmp(typeOfField,"uniform")==0) {
	    sscanf(&line[0],"%*s %*s %*s %*s %s %g",nameOfField,&fieldValue);
             sr1Field = new sr1MagneticField();
          }
          else if (strcmp(typeOfField,"gaussian")==0) {
//           sscanf(&line[0],"%*s %*s %*s %*s %s %g %g",nameOfField,&fieldValue,&sigma);
//           sr1Field = new sr1MagneticField();
          }
	  else if (strcmp(typeOfField,"fromfile")==0) {
	    char fieldInputFileName[100];
	    char fieldTableType[100];
	    sscanf(&line[0],"%*s %*s %*s %*s %s %s %s %g",fieldTableType,fieldInputFileName,nameOfField,&fieldValue);

	    if (strcmp(fieldTableType,"2D")==0) {
	      sr1Field = new sr1TabulatedField2D(fieldInputFileName, fieldValue*tesla, 1*cm, 0.00001);
	    }
	    else {
	      sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%s\" .",fieldTableType);
	      sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	    }


	  }
          else  ReportGeometryProblem(line);


          pFieldMgr->SetDetectorField(sr1Field);

          G4Mag_SpinEqRhs *Mag_SpinEqRhs;
          G4MagIntegratorStepper* pStepper;
          Mag_SpinEqRhs = new G4Mag_SpinEqRhs(sr1Field);

          pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12);
          G4ChordFinder *pChordFinder = new G4ChordFinder(sr1Field,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 sr1 */
	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="<<orientation<<G4endl;
	    ReportGeometryProblem(line);
	  }

	  solid =  new G4IntersectionSolid(solidName, solidDetTube, solidDetBox, yRot, zTrans);
	}
	else ReportGeometryProblem(line);

	G4ThreeVector position = G4ThreeVector (posx*mm,posy*mm,posz*mm);
	//	G4cout << "New volume: "<<tmpString2<<" "<<name<<", solid geometry parameters="<<x1<<" "<<x2<<" "<<x3<<" "<<x4<<" "<<x5
	//	       << ", position="<<position<<", mother="<<mothersName<<G4endl;


	G4LogicalVolume* mothersVolume;
	if (strcmp(name,"World")==0) { mothersVolume=NULL; }
	else {
	  mothersVolume	= FindLogicalVolume(mothersName);
	  if (mothersVolume==NULL) {
	    G4cout << "ERROR! sr1DetectorConstruction::Construct():"<<G4endl;
	    G4cout << "       Logical Volume \"" << mothersName <<"\" not found!"<<G4endl;
	    G4cout << "       (requested for the definition of the volume \""<<name<<"\")"<<G4endl;
	    G4cout << "S T O P    F O R C E D"<<G4endl;
	    exit(1);
	  }
	}
	
	G4cout<<"Vol_ID = "<<volumeID<<"\t ";	
	G4RotationMatrix* pRot = pointerToRotationMatrix[rotMatrix];
	//	G4cout << "rotMatrix="<<rotMatrix<<"  pRot="<<pRot<<"        ";
	if (pRot!=NULL) {G4cout<<"    rotated, ";}          // OK, rotation matrix found.
	else if ((strcmp(rotMatrix,"norot"))) {             // Problem, matrix not found.
	  G4cout <<"Rotation Matrix \""<<rotMatrix<<"\" not found (not defined at this point)!"
		 <<"    Check the geometry!"<<G4endl;
	  G4cout <<"S T O P    F O R C E D"<<G4endl;
	  ReportGeometryProblem(line);
	}
	else {G4cout<<"non-rotated, ";}
	
	if (strcmp(sensitiveDet,"dead")) {G4cout<<"    sensitive volume:  ";}
	else {G4cout<<"non-sensitive volume:  ";}
	
	G4FieldManager* pFieldMan = pointerToField[actualFieldName];
	if (pFieldMan!=NULL) {G4cout<< actualFieldName << " ON, ";}
	else if ((strcmp(actualFieldName,"nofield"))) {
	  G4cout <<"Field Manager \""<<actualFieldName<<"\" not found (not defined at this point)!"
		 <<"    Check the geometry!"<<G4endl;
	  G4cout <<"S T O P    F O R C E D"<<G4endl;
	  ReportGeometryProblem(line);
        }

	G4Material* MATERIAL=CharToMaterial(material);
	G4String logicName="log_";  logicName+=name;
	G4LogicalVolume* logicVolume = new G4LogicalVolume(solid,MATERIAL,logicName,pFieldMan,0,0);
	//	G4cout << "sr1DetectorConstruction.cc: logicVolume created ="<<logicVolume<<G4endl;
	G4String physName="phys_";  physName+=name;
	G4VPhysicalVolume* physiVolume = new G4PVPlacement(pRot,              // rotation 
					position,
					logicVolume,    // its logical volume
					physName,       // its name
					mothersVolume,  // its mother  volume
					false,           // no boolean operations
					0,               // no field specific to volume
					checkOverlap);  // Check for overlap
	//	G4cout << "sr1DetectorConstruction.cc physiVolume created ="<<physiVolume<<G4endl;
	if (strcmp(name,"World")==0)  { pointerToWorldVolume=physiVolume; G4cout<<"World Volume"<<endl;}
	
	// If the name of the volume contains the string "save", it will be treated specifically 
        // (a position of where a particle went through the volume will be saved)
	char tmpMagicKeyword[5] = "save"; 
	if (CheckIfStringContainsSubstring(name,tmpMagicKeyword)) {
	  if (volumeID!=0) {    // do not store hits in special "save" volume if ID of this volume is 0 
	                        // (due to difficulties to distinguish between ID=0 and no save volume when using std::map)
	    sr1SteppingAction* mySteppingAction = sr1SteppingAction::GetInstance();
	    mySteppingAction->SetLogicalVolumeAsSpecialSaveVolume(logicName,volumeID);
	    sr1RootOutput::GetRootInstance()->SetSpecialSaveVolumeDefined();
	  }
	}

	// Unless specified as "dead", set the volume sensitive volume:
	if (strcmp(sensitiveDet,"dead")) {
	  if (strcmp(sensitiveDet,"sr1/ScintSD")==0) {
	    if(!aScintSD) {
	      G4SDManager* SDman = G4SDManager::GetSDMpointer();
	      G4String scintSDname = sensitiveDet;
	      aScintSD = new sr1ScintSD( scintSDname );
	      SDman->AddNewDetector( aScintSD );
	      G4cout<<"sr1DetectorConstruction.cc:  aScintSD added: "<<aScintSD->GetFullPathName()<<G4endl;
	    }
	    logicVolume->SetSensitiveDetector( aScintSD );
	  }
	  else {
	    G4cout<<" sr1DetectorConstruction.cc: unknown sensitive detector \""<<sensitiveDet<<"\" requested!"<<G4endl;
	    G4cout<<"\""<<line<<"\""<<G4endl;
	    G4cout<<"        S T O P    F O R C E D !!!"<<G4endl;
	    exit(1);
	  }
	}

	//  Set the volume ID (used only for the Root output).
	if (volumeID!=0) {
	  myRootOutput->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! sr1DetectorConstruction::Construct():   Logical Volume \""
		   << tmpString2 <<"\" not found!"<<G4endl<<"S T O P    F O R C E D"<<G4endl;
	    exit(1);
	  }	  
	  SetColourOfLogicalVolume(pLogVol,tmpString3);
	}
	else {
	  G4LogicalVolumeStore* pLogStore = G4LogicalVolumeStore::GetInstance();
	  if (pLogStore==NULL) {
	    G4cout<<"ERROR:  sr1DetectorConstruction.cc:    G4LogicalVolumeStore::GetInstance()  not found!"<<G4endl;
	  }
	  else {
	    for (unsigned int i=0; i<pLogStore->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: sr1DetectorConstruction.cc   unknown G4VisAttributes requested:"<<G4endl;
	//	  G4cout<<"\""<<line<<"\""<<G4endl;
	//	  G4cout<<"    Ignoring and continuing"<<G4endl;
	//	}
	//	      G4VisAttributes* tmpVisAttributes = new G4VisAttributes(visibility,colour);
	//	      pLogVol->SetVisAttributes(tmpVisAttributes);
      }
      
      //      else if (strcmp(tmpString1,"setsensitive")==0){
      //	G4SDManager* SDman = G4SDManager::GetSDMpointer();
      //	if (strcmp(tmpString2,"sr1/ScintSD")==0) {
      //	  G4String scintSDname = tmpString2;   //"sr1/ScintSD";
      //	  if(!aScintSD) {
      //	    aScintSD = new sr1ScintSD( scintSDname );
      //	    SDman->AddNewDetector( aScintSD );
      //	    G4cout<<"sr1DetectorConstruction.cc:  aScintSD added: "<<aScintSD->GetFullPathName()<<G4endl;
      //	  }
      //	  int volumeID;	  
      //	  sscanf(&line[0],"%*s %*s %*s %s %d",tmpString3,&volumeID);
      //	  G4LogicalVolume* pLogVol = FindLogicalVolume(tmpString3);
      //	  if (pLogVol==NULL) {
      //	    G4cout << "ERROR! sr1DetectorConstruction::Construct():   Logical Volume \""
      //		   << tmpString3 <<"\" not found!"<<G4endl<<"S T O P    F O R C E D"<<G4endl;
      //	    exit(1);
      //	  }
      //	  pLogVol->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<<" sr1DetectorConstruction.cc: unknown sensitive detector \""<<tmpString2<<"\" requested!"<<G4endl;
      //	  G4cout<<"\""<<line<<"\""<<G4endl;
      //	  G4cout<<"        S T O P    F O R C E D !!!"<<G4endl;
      //	  exit(1);
      //	}
      //      }

      
      //! VERY IMPORTANT: THIS PART HAS BEEN CHANGED TO ALLOW FOR THE GUMPLINGER CODE! TS
      
      else if (strcmp(tmpString1,"globalfield")==0){
	//check that the field has not been defined in some other way
	if (strcmp(howIsTheFieldDefined,"DIFFERENT_FIELDS")==0) {
	  sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): Field already defined!  howIsTheFieldDefined=\"%s\" .",
		  howIsTheFieldDefined);
	  sr1ErrorMessage::GetInstance()->sr1Error(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,"sr1DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%c\" .",fieldType);
	    sr1ErrorMessage::GetInstance()->sr1Error(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,"sr1DetectorConstruction.cc::Construct(): GLOBAL FIELD: Logical volume \"%s\" not found.",
		    logicalVolumeName);
	    sr1ErrorMessage::GetInstance()->sr1Error(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="<<physicalVolumeName<<G4endl;
	  //     Now find the physical volume and its position within its mother:
	  ///G4VPhysicalVolume* physVol = G4PhysicalVolumeStore::GetInstance()->GetVolume(physicalVolumeName);
	  ///if (physVol==NULL) {
	  ///  sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): GLOBAL FIELD: Physical volume \"%s\" not found.",
	  ///	    physicalVolumeName.c_str());
	  ///  sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	  ///}
	  ///G4ThreeVector position = physVol->GetObjectTranslation();
	  ///G4cout<<"position of the volume "<<physicalVolumeName<<" in its mother is: "
	  ///	<<position.x()<<", "<<position.y()<<", "<<position.z()<<G4endl;
	  
	  // Now call the relevant electric field //magnetic field
	  if (strncmp(fieldTableType,"2D",2)==0) {
	  ///if (strcmp(fieldTableType,"2D")==0) {
	    //sr1TabulatedElementField2D*  myElementTableField = 
	      //new sr1TabulatedElementField2D(fieldInputFileName, fieldValue*tesla, 1*cm, 0.00001, logVol, position);
			  
             // decimeters, volts/m factor = 1E-5, and offset = -567 mm //Default unit in G4 is MegaVolt
	     
	     
	     ///if (strncmp(fieldTableType,"f",4)==0) 
	     ///
	     ///found=str.find('.');
	     //std::string::size_type pos = fieldTableType.find("f");
	     
	     char fieldFolded = fieldTableType[3]; // In case of symmetric (i.e. folded) 2D axial field
	     if (fieldFolded == 'f') {
	     sr1TabulatedElementField2Df*  myElementTableField = 
	         new sr1TabulatedElementField2Df(fieldInputFileName, fieldType, fieldValue*EMfieldUnit, logVol, position);
	     myElementTableField->SetElementFieldName(tmpString2);
	      if (fieldNrOfSteps>0) {
	        myElementTableField->SetEventNrDependentField(fieldValue*EMfieldUnit,fieldValueFinal*EMfieldUnit,fieldNrOfSteps);
	     }
	     }
	     //* --- Changed by TS to check error in field steps
	     /* 
	    sr1Axial2DElField*  myElementTableField =	    /// CHECK WHETHER CORRECT!! TS
	      new sr1Axial2DElField(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()="<<fields->size()<<G4endl;
	    //	    }
             */	  
	    else if (fieldFolded != 'f') { // 2D axial field
	     sr1TabulatedElementField2D*  myElementTableField = 
	        new sr1TabulatedElementField2D(fieldInputFileName, fieldType, fieldValue*EMfieldUnit, logVol, position);
	    myElementTableField->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) {
            sr1TabulatedElementField3D*  myElementTableField =
	   ///new sr1TabulatedElementField3D(fieldInputFileName, fieldValue*tesla, 1*m, 1., logVol, position);
              new sr1TabulatedElementField3D(fieldInputFileName, fieldType, fieldValue*EMfieldUnit, logVol, position);
            myElementTableField->SetElementFieldName(tmpString2);
            if (fieldNrOfSteps>0) {
              myElementTableField->SetEventNrDependentField(fieldValue*EMfieldUnit,fieldValueFinal*EMfieldUnit,fieldNrOfSteps);
            }
	  }
	  else {
	    sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%s\" .",fieldTableType);
	    sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	  }
	}

	
	else if (strcmp(typeOfField,"uniform")==0) {
	   //// G4cout<<"Uniform field through global F04 call."<<G4endl;
          
/*!	  if (strcmp(typeOfField,"uniform")==0) {
	    sscanf(&line[0],"%*s %*s %*s %*s %s %g",nameOfField,&fieldValue);
             sr1Field = new sr1MagneticField();
          }
	  else if (strcmp(typeOfField,"fromfile")==0) {
	    char fieldInputFileName[100];
	    char fieldTableType[100];
	    sscanf(&line[0],"%*s %*s %*s %*s %s %s %s %g",      fieldTableType,fieldInputFileName,nameOfField,&fieldValue);
///	    /sr1/command 1 globalfield 2 MCPSfield 3 fromfile 4        2D_E 5        OI_10T_grid_xz.table 6 log_MCPV 7 0.004 8
	    
	    
	    if (strcmp(fieldTableType,"2D")==0) {
	      sr1Field = new sr1TabulatedField2D(fieldInputFileName, fieldValue*tesla, 1*cm, 0.00001);
	    }
	    else {
	      sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): unknown type of the tabulated field \"%s\" .",fieldTableType);
	      sr1ErrorMessage::GetInstance()->sr1Error(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,"sr1DetectorConstruction.cc::Construct(): GLOBAL FIELD: Logical volume \"%s\" not found.", logicalVolumeName);
	    sr1ErrorMessage::GetInstance()->sr1Error(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="<<physicalVolumeName<<G4endl;
	  // Now find the physical volume and its position within its mother:
	  G4VPhysicalVolume* physVol = G4PhysicalVolumeStore::GetInstance()->GetVolume(physicalVolumeName);
	  if (physVol==NULL) {
	    sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): GLOBAL FIELD: Physical volume \"%s\" not found.",
		    physicalVolumeName.c_str());
	    sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	  }
	  G4ThreeVector position = physVol->GetObjectTranslation();
	  G4cout<<"The position of the volume "<<physicalVolumeName<<" in its mother is: "
		<<position.x()<<", "<<position.y()<<", "<<position.z()<<G4endl;
	  
	*/
		
		
		
		
	  // Now call the relevant magnetic field
	  G4double fieldValue_tmp[6] = {
	  fieldValue[0]*tesla,        fieldValue[1]*tesla,        fieldValue[2]*tesla,
	  fieldValue[3]*(kilovolt/mm),fieldValue[4]*(kilovolt/mm),fieldValue[5]*(kilovolt/mm)};
	  
	  sr1UniformField* myElementField = 
	      new sr1UniformField(fieldValue_tmp, logVol, position);
	      ///new sr1UniformField(fieldValue*tesla, logVol, position);
	  ///G4cout<<"The LOGICAL volume is named "<<logVol->GetName()<<G4endl;
	  myElementField->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()="<<fields->size()<<G4endl;
	    	          }
	  else {
	    sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): unknown type of error?.");
	    sr1ErrorMessage::GetInstance()->sr1Error(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<<"sr1DetectorConstruction.cc:   Uniform Magnetic field set to "<<fieldValue<<" T."<<G4endl;
	  /// G4cout<<"nr. of param="<<nParameters<<G4endl;
	  
	  ///nParameters=sscanf(&line[0],"%*s %*s %*s %s %s %g %g %d",
          ///                            fieldTableType,fieldInputFileName,&fieldValue,&fieldValueFinal,&fieldNrOfSteps);
	  ///SetMagField(fieldTableType,fieldInputFileName,fieldValue*tesla);
	 /// }
	  
	  
	}


	else if (strcmp(tmpString2,"setparameter")==0){
	  // First check that the magnetic field already exists:
	  G4FieldManager*      fieldMgr  = G4TransportationManager::GetTransportationManager()->GetFieldManager();
	  G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField();
	  if (fieldMgr==NULL) {
	    ReportProblemInStearingFile(line);
	    sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): G4FieldManager not found: fieldMgr=NULL");
	    sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	  }
	  if (propagMgr==NULL) {
	    ReportProblemInStearingFile(line);
	    sprintf(eMessage,"sr1DetectorConstruction.cc::Construct(): G4PropagatorInField not found: propagMgr=NULL");
	    sr1ErrorMessage::GetInstance()->sr1Error(FATAL,eMessage,false);
	  }
	  else {
	    char parameterName[100];
	    float parameterValue;
	    sscanf(&line[0],"%*s %*s %*s %s %g",parameterName,&parameterValue);
	    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<<"sr1DetectorConstruction.cc: ERROR: Unknown parameterName \""
		    <<parameterName<<"\" ."<<G4endl;
	      ReportGeometryProblem(line);
	    }
	  }
	}

	else if (strcmp(tmpString2,"printparameters")==0){    // Print the accuracy parameters for the magnetic field
	  G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
	  G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField();
	  if (fieldMgr==NULL) {
	    G4cout<<"sr1DetectorConstruction: ERROR: Field manager not found!"<<G4endl;
	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
	  }
	  else { 
	    G4cout<<"GetDeltaIntersection()    = "<<fieldMgr->GetDeltaIntersection()/mm<<" mm"<<G4endl;
	    G4cout<<"GetDeltaOneStep()         = "<<fieldMgr->GetDeltaOneStep()/mm<<" mm"<<G4endl;
	    G4cout<<"GetMinimumEpsilonStep()   = "<<fieldMgr->GetMinimumEpsilonStep()<<G4endl;
	    G4cout<<"GetMaximumEpsilonStep()   = "<<fieldMgr->GetMaximumEpsilonStep()<<G4endl;
	  }

	  if (propagMgr==NULL) {
	    G4cout<<"sr1DetectorConstruction: ERROR: G4PropagatorInField not found!"<<G4endl;
	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
	  }
	  else {
	    G4cout<<"GetLargestAcceptableStep()= "<<propagMgr->GetLargestAcceptableStep()/mm<<" mm"<<G4endl;
	    G4cout<<"GetMaxLoopCount()         = "<<propagMgr->GetMaxLoopCount()<<G4endl;
	  }	    
	}

	else if (strcmp(tmpString2,"printFieldValueAtPoint")==0){    // Print the fieldvalue at the given point
	  float p0, p1, p2;
	  sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2);
	  ///double point[4]={p0,p1,p2,0};  double Bfi[6]={0,0,0,0,0,0}; //See comment below! TS
	  if (F04GlobalField::Exists()) {
	    F04GlobalField* myGlobalField = F04GlobalField::getObject();
	    if (myGlobalField) {
	      // The global field is not properly initialised at this point.  Therefore the points have to
              // be stored in  "F04GlobalField" object  and printed later on in musrRunAction.cc .
	      myGlobalField->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,"sr1DetectorConstruction.cc::Construct(): printFieldValueAtPoint requested, but field not found");
	      sr1ErrorMessage::GetInstance()->sr1Error(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,"sr1DetectorConstruction.cc::Construct(): Field is already defined!  howIsTheFieldDefined=\"%s\" .",
		  howIsTheFieldDefined);
	  sr1ErrorMessage::GetInstance()->sr1Error(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<<"sr1DetectorConstruction.cc:   Uniform Magnetic field set to "<<fieldValue<<" T."<<G4endl;
	  G4cout<<"nr. of param="<<nParameters<<G4endl;
	}
	else if (strcmp(tmpString2,"gaussian")==0){
	  fieldOption=2;
	  float sigma=1000;
	  nParameters=sscanf(&line[0],"%*s %*s %*s %g %g %g %d",&fieldValue,&sigma,&fieldValueFinal,&fieldNrOfSteps);
	  G4ThreeVector inputParam=G4ThreeVector(fieldValue,sigma,0.);
	  SetGaussianMagField(inputParam);
	  G4cout<<"sr1DetectorConstruction.cc:   Gaussian Magnetic field set to "<<fieldValue
		<<" T. with sigma= "<<sigma<<" mm."<<G4endl;
	}
	else if (strcmp(tmpString2,"fromfile")==0){
	  fieldOption=3;
	  char fieldInputFileName[100];
	  char fieldTableType[100];
	  nParameters=sscanf(&line[0],"%*s %*s %*s %s %s %g %g %d",
			     fieldTableType,fieldInputFileName,&fieldValue,&fieldValueFinal,&fieldNrOfSteps);
	  SetMagField(fieldTableType,fieldInputFileName,fieldValue*tesla);
	  G4cout<<"sr1DetectorConstruction.cc:  Tabulated magnetic field set to "<<fieldValue<<" T."<<G4endl;
	}
	else if (strcmp(tmpString2,"setOffsetOf2DField")==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 (sr1TabulatedField2D::GetInstance()==NULL) {G4cout<<"QQQQQQQQQQQQQQQQQQQQQQQQQQQ"<<G4endl;}
	  else { sr1TabulatedField2D::GetInstance()->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 (sr1TabulatedField3D::GetInstance()==NULL) {G4cout<<"WWWWWWWWWWWWWWWWWWWWWWWWWW"<<G4endl;}
	  else { sr1TabulatedField3D::GetInstance()->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<<G4endl<<"sr1DetectorConstruction.cc: ERROR: Field manager not found!"<<G4endl;
	    G4cout<<"Check that the magnetic field has been defined before setting"
		  <<" the accuracy parameters!"<<G4endl;
	    ReportGeometryProblem(line);
	  }
	  else if (propagMgr==NULL) {
	    G4cout<<"sr1DetectorConstruction.cc: ERROR: G4PropagatorInFile not found!"<<G4endl;
	    G4cout<<"Check that the magnetic field has been defined before setting"
		  <<" the accuracy parameters!"<<G4endl;
	    ReportGeometryProblem(line);
	  }
	  else {
	    // Now set the required parameter
	    char parameterName[100];
	    float parameterValue;
	    sscanf(&line[0],"%*s %*s %*s %s %g",parameterName,&parameterValue);
	    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<<"sr1DetectorConstruction.cc: ERROR: Unknown parameterName \""
		    <<parameterName<<"\" ."<<G4endl;
	      ReportGeometryProblem(line);
	    }
	  }
	}
	else if (strcmp(tmpString2,"printparameters")==0){    // Print the accuracy parameters for the magnetic field
	  G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager();
	  G4PropagatorInField* propagMgr = G4TransportationManager::GetTransportationManager()->GetPropagatorInField();
	  if (fieldMgr==NULL) {
	    G4cout<<"sr1DetectorConstruction: ERROR: Field manager not found!"<<G4endl;
	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
	  }
	  else { 
	    G4cout<<"GetDeltaIntersection()    = "<<fieldMgr->GetDeltaIntersection()/mm<<" mm"<<G4endl;
	    G4cout<<"GetDeltaOneStep()         = "<<fieldMgr->GetDeltaOneStep()/mm<<" mm"<<G4endl;
	    G4cout<<"GetMinimumEpsilonStep()   = "<<fieldMgr->GetMinimumEpsilonStep()<<G4endl;
	    G4cout<<"GetMaximumEpsilonStep()   = "<<fieldMgr->GetMaximumEpsilonStep()<<G4endl;
	  }

	  if (propagMgr==NULL) {
	    G4cout<<"sr1DetectorConstruction: ERROR: G4PropagatorInField not found!"<<G4endl;
	    G4cout<<"                          Can not print the accuracy parameters of the magnetic field."<<G4endl;
	  }
	  else {
	    G4cout<<"GetLargestAcceptableStep()= "<<propagMgr->GetLargestAcceptableStep()/mm<<" mm"<<G4endl;
	    G4cout<<"GetMaxLoopCount()         = "<<propagMgr->GetMaxLoopCount()<<G4endl;
	  }	    
	}
	else if (strcmp(tmpString2,"printFieldValueAtPoint")==0){    // Print the fieldvalue at the given point
	  float p0, p1, p2;
	  sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2);
	  double point[4]={p0,p1,p2,0};  double Bfi[3]={0,0,0};
	  G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->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.
	sr1EventAction* pointerToEventAction = sr1EventAction::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! sr1DetectorConstruction::Construct():   Logical Volume \""
		 << tmpString3 <<"\" not found!"<<G4endl<<"S T O P    F O R C E D"<<G4endl;
	  ReportGeometryProblem(line);
	}
	float minEnergy;
	sscanf(&line[0],"%*s %*s %*s %g",&minEnergy);
	if (minEnergy<1e-2) {
	  G4cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<G4endl;
	  G4cout<<"sr1DetectorConstruction: WARNING:  minEnergy too small!!  = "<<minEnergy<<G4endl;
	  G4cout<<"                                      ==>   minEnergy value ignored"<<G4endl;
	  G4cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<G4endl;
	}
	else {
	  G4UserLimits* User_lim = new G4UserLimits();
	  User_lim->SetUserMinEkine(minEnergy*eV);
	  pLogVol->SetUserLimits(User_lim);
	  G4cout<<"\nsr1DetectorConstruction.cc:  Energy limit in "<<tmpString2<<" set to "
		<<minEnergy<<" eV."<<G4endl;
	}
      }

      //! End of feature
      
     
      //   Set range cut for a given volume, if requested by the macro file
      else if (strcmp(tmpString1,"SetUserLimits")==0){
	G4LogicalVolume* pLogVol = FindLogicalVolume(tmpString2);
	if (pLogVol==NULL) {
	  G4cout << "ERROR! sr1DetectorConstruction::Construct():   Logical Volume \""
		 << tmpString3 <<"\" not found!"<<G4endl<<"S T O P    F O R C E D"<<G4endl;
	  ReportGeometryProblem(line);
	}
	float maxStep;
	sscanf(&line[0],"%*s %*s %*s %g",&maxStep);
	if (maxStep<0.000001*mm) {
	  G4cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<G4endl;
	  G4cout<<"sr1DetectorConstruction: WARNING:  maxStep very small!!  ="<<maxStep<<G4endl;
	  G4cout<<"                                      ==>   maxStep value ignored"<<G4endl;
	  G4cout<<" !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"<<G4endl;
          // cdel:   put here by mistake ?	  pLogVol->SetSensitiveDetector( aScintSD );
	}
	else {
	  pLogVol->SetUserLimits(new G4UserLimits(maxStep*mm));
	  G4cout<<"sr1DetectorConstruction.cc:  Range cut in "<<tmpString2<<" set to "
		<<maxStep<<" mm."<<G4endl;
	}
      }

      else if (strcmp(tmpString1,"typeofprocesses")==0){
	sr1Parameters* myParameters = sr1Parameters::GetInstance();
	myParameters -> SetMyTypeOfProcesses(tmpString2);
      }

      else if (strcmp(tmpString1,"storeOnlyEventsWithHits")==0){
	if (strcmp(tmpString2,"false")==0){ sr1Parameters::storeOnlyEventsWithHits = false; }
      } 

      else if (strcmp(tmpString1,"storeOnlyTheFirstTimeHit")==0){
	if (strcmp(tmpString2,"true")==0){ sr1Parameters::storeOnlyTheFirstTimeHit = true; }
      } 

      else if (strcmp(tmpString1,"signalSeparationTime")==0){
	float timeSeparation;
	sscanf(&line[0],"%*s %*s %g",&timeSeparation);
        sr1Parameters::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){ 
	sr1Parameters::includeMuoniumProcesses = false; 
	G4cout<<"\nIgnoring Muonium formation and other Mu processes.\n"<<G4endl;
	}
      }

/*!      else if (strcmp(tmpString1,"killAllPositrons")==0){
	if (strcmp(tmpString2,"true")==0){ sr1Parameters::killAllPositrons = true; }
      } 

      else if (strcmp(tmpString1,"G4GeneralParticleSource")==0){
	// This parameter is read in in the constructor of "sr1Parameters.cc", so ignore it here.
	//  if (strcmp(tmpString2,"true")==0){ sr1Parameters::boolG4GeneralParticleSource = true; }
      } 
*/
      
      //  Set the variables which the user does not want to store in the output root file.
      else if (strcmp(tmpString1,"rootOutput")==0){
	char storeIt[100];
	sscanf(&line[0],"%*s %*s %*s %s",storeIt);
	if (strcmp(storeIt,"off")==0) {
	  if (strcmp(tmpString2,"runID")==0)        {sr1RootOutput::store_runID = false;}
	  if (strcmp(tmpString2,"eventID")==0)      {sr1RootOutput::store_eventID = false;}
	  if (strcmp(tmpString2,"BFieldAtDecay")==0){sr1RootOutput::store_BFieldAtDecay = false;}
	  if (strcmp(tmpString2,"muIniPosX")==0)    {sr1RootOutput::store_muIniPosX = false;}
	  if (strcmp(tmpString2,"muIniPosY")==0)    {sr1RootOutput::store_muIniPosY = false;}
	  if (strcmp(tmpString2,"muIniPosZ")==0)    {sr1RootOutput::store_muIniPosZ = false;}
	  if (strcmp(tmpString2,"muIniMomX")==0)    {sr1RootOutput::store_muIniMomX = false;}
	  if (strcmp(tmpString2,"muIniMomY")==0)    {sr1RootOutput::store_muIniMomY = false;}
	  if (strcmp(tmpString2,"muIniMomZ")==0)    {sr1RootOutput::store_muIniMomZ = false;}
	  if (strcmp(tmpString2,"muIniPolX")==0)    {sr1RootOutput::store_muIniPolX = false;}
	  if (strcmp(tmpString2,"muIniPolY")==0)    {sr1RootOutput::store_muIniPolY = false;}
	  if (strcmp(tmpString2,"muIniPolZ")==0)    {sr1RootOutput::store_muIniPolZ = false;}
	  if (strcmp(tmpString2,"muDecayDetID")==0) {sr1RootOutput::store_muDecayDetID = false;}
	  if (strcmp(tmpString2,"muDecayPosX")==0)  {sr1RootOutput::store_muDecayPosX = false;}
	  if (strcmp(tmpString2,"muDecayPosY")==0)  {sr1RootOutput::store_muDecayPosY = false;}
	  if (strcmp(tmpString2,"muDecayPosZ")==0)  {sr1RootOutput::store_muDecayPosZ = false;}
	  if (strcmp(tmpString2,"muDecayTime")==0)  {sr1RootOutput::store_muDecayTime = false;}
	  if (strcmp(tmpString2,"muDecayPolX")==0)  {sr1RootOutput::store_muDecayPolX = false;}
	  if (strcmp(tmpString2,"muDecayPolY")==0)  {sr1RootOutput::store_muDecayPolY = false;}
	  if (strcmp(tmpString2,"muDecayPolZ")==0)  {sr1RootOutput::store_muDecayPolZ = false;}
	  if (strcmp(tmpString2,"muTargetTime")==0) {sr1RootOutput::store_muTargetTime = false;}
	  if (strcmp(tmpString2,"muTargetPosX")==0) {sr1RootOutput::store_muTargetPosX = false;}
	  if (strcmp(tmpString2,"muTargetPosY")==0) {sr1RootOutput::store_muTargetPosY = false;}
	  if (strcmp(tmpString2,"muTargetPosZ")==0) {sr1RootOutput::store_muTargetPosZ = false;}
	  if (strcmp(tmpString2,"muTargetPolX")==0) {sr1RootOutput::store_muTargetPolX = false;}
	  if (strcmp(tmpString2,"muTargetPolY")==0) {sr1RootOutput::store_muTargetPolY = false;}
	  if (strcmp(tmpString2,"muTargetPolZ")==0) {sr1RootOutput::store_muTargetPolZ = false;}
	  if (strcmp(tmpString2,"posIniMomX")==0)   {sr1RootOutput::store_posIniMomX = false;}
	  if (strcmp(tmpString2,"posIniMomY")==0)   {sr1RootOutput::store_posIniMomY = false;}
	  if (strcmp(tmpString2,"posIniMomZ")==0)   {sr1RootOutput::store_posIniMomZ = false;}
	  if (strcmp(tmpString2,"globalTime")==0)   {sr1RootOutput::store_globalTime = false;}
	  if (strcmp(tmpString2,"fieldValue")==0)   {sr1RootOutput::store_fieldValue = false;}
	  if (strcmp(tmpString2,"fieldNomVal")==0)  {sr1RootOutput::store_fieldNomVal = false;}
	  if (strcmp(tmpString2,"det_ID")==0)       {sr1RootOutput::store_det_ID = false;}
	  if (strcmp(tmpString2,"det_edep")==0)     {sr1RootOutput::store_det_edep = false;}
	  if (strcmp(tmpString2,"det_edep_el")==0)  {sr1RootOutput::store_det_edep_el = false;}
	  if (strcmp(tmpString2,"det_edep_pos")==0) {sr1RootOutput::store_det_edep_pos = false;}
	  if (strcmp(tmpString2,"det_edep_gam")==0) {sr1RootOutput::store_det_edep_gam = false;}
	  if (strcmp(tmpString2,"det_edep_mup")==0) {sr1RootOutput::store_det_edep_mup = false;}
	  if (strcmp(tmpString2,"det_nsteps")==0)   {sr1RootOutput::store_det_nsteps = false;}
	  if (strcmp(tmpString2,"det_length")==0)   {sr1RootOutput::store_det_length = false;}
	  if (strcmp(tmpString2,"det_start")==0)    {sr1RootOutput::store_det_start = false;}
	  if (strcmp(tmpString2,"det_end")==0)      {sr1RootOutput::store_det_end = false;}
	  if (strcmp(tmpString2,"det_x")==0)        {sr1RootOutput::store_det_x = false;}
	  if (strcmp(tmpString2,"det_y")==0)        {sr1RootOutput::store_det_y = false;}
	  if (strcmp(tmpString2,"det_z")==0)        {sr1RootOutput::store_det_z = false;}
	}
      }

      else ReportGeometryProblem(line);
      
    }
  }
  fclose(fSteeringFile);
  G4cout<< "sr1DetectorConstruction.cc:  pointerToWorldVolume="<<pointerToWorldVolume<<G4endl;
  return pointerToWorldVolume;
}


void sr1DetectorConstruction::DefineMaterials()
{
//--------- Material definitions ---------

  G4String symbol;         // a = mass of a mole;
  //G4double a, z, density;  // z = number of protons;  

  G4int natoms, ncomponents;
  G4double fractionmass, density;
  
  G4NistManager* man = G4NistManager::Instance();
  
  // Define elements and compounds needed for building materials
  
  // elements required for Stainless Steel
  G4Element* Cr = man->FindOrBuildElement("Cr");
  G4Element* Fe = man->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."<<G4endl;
  


//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......


#include "G4RunManager.hh"

void sr1DetectorConstruction::UpdateGeometry()
{
  G4RunManager::GetRunManager()->DefineWorldVolume(Construct());
}

void sr1DetectorConstruction::SetUniformMagField(G4double fieldValue)
{     // If we want a uniform field along z      

      G4FieldManager   *pFieldMgr;
      sr1MagneticField*  sr1UniformField= new sr1MagneticField();
          sr1UniformField->SetFieldValue(fieldValue);

      pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();
      G4Mag_SpinEqRhs *Mag_SpinEqRhs;
      G4MagIntegratorStepper* pStepper;
      Mag_SpinEqRhs = new G4Mag_SpinEqRhs(sr1UniformField);
      pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12);
      G4cout<< "DeltaStep "<<pFieldMgr->GetDeltaOneStep()/mm <<"mm" <<G4endl;
      //      G4cout<< "GetLargestAcceptableStep()="<<pFieldMgr->GetLargestAcceptableStep()/mm<<"mm"<<G4endl;      
      G4ChordFinder *pChordFinder = new G4ChordFinder(sr1UniformField,0.01*mm,pStepper);
      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->SetDetectorField(sr1UniformField);
      
      // store the pointer to the sr1MagneticField in the event action, which is needed for
      // the time dependent magnetic field. 
      sr1EventAction* pointerToEventAction = sr1EventAction::GetInstance();
      pointerToEventAction->SetPointerToMusrUniformField(sr1UniformField);

      //   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()="<<LargestAcceptableStep/mm<<"mm"<<G4endl;
}


void sr1DetectorConstruction::SetMagField(const char* fieldTableType, const char* inputFileName, G4double fieldValue )
{

//------------ Magnetic field from file ----------------------------
      G4FieldManager   *pFieldMgr;
      //      //      G4MagneticField* sr1Field = new sr1TabulatedField3D(inputFileName, fieldValue);
      // store the pointer to the sr1MagneticField in the sr1EventAction, which is needed for
      // the time dependent magnetic field. 
      sr1EventAction* pointerToEventAction = sr1EventAction::GetInstance();

      G4MagneticField* sr1Field = NULL;
      if (strcmp(fieldTableType,"3D")==0) {
	sr1TabulatedField3D* myMusrField = new sr1TabulatedField3D(inputFileName, fieldValue);
	sr1Field = myMusrField;
	pointerToEventAction->SetPointerToTabulatedField3D(myMusrField);
      }
      else if (strcmp(fieldTableType,"2D")==0) {
	sr1TabulatedField2D* myMusrField = new sr1TabulatedField2D(inputFileName, fieldValue, 1*cm, 0.00001);
	sr1Field = myMusrField;
	pointerToEventAction->SetPointerToTabulatedField2D(myMusrField);
      }
      else {
	char eMessage[200];
	sprintf(eMessage,"sr1DetectorConstruction.cc::SetMagField(): unknown type of the tabulated field \"%s\" .",fieldTableType);
	sr1ErrorMessage::GetInstance()->sr1Error(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};
      //        sr1Field->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(sr1Field);
      pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12);
      G4cout<< "DeltaStep "<<pFieldMgr->GetDeltaOneStep()/mm <<"mm" <<endl;
 
      G4ChordFinder *pChordFinder = new G4ChordFinder(sr1Field,0.01*mm,pStepper);
      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->SetDetectorField(sr1Field);
}



// void sr1DetectorConstruction::SetGaussianMagField(G4double fieldValue, G4double fieldSigma, G4double dummy) {
void sr1DetectorConstruction::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* sr1Field= new sr1GaussianField(fieldValue,fieldSigma);
      pFieldMgr=G4TransportationManager::GetTransportationManager()->GetFieldManager();

      G4Mag_SpinEqRhs *Mag_SpinEqRhs;
      G4MagIntegratorStepper* pStepper;

      Mag_SpinEqRhs = new G4Mag_SpinEqRhs(sr1Field);
      pStepper = new G4ClassicalRK4(Mag_SpinEqRhs,12);
      G4cout<< "DeltaStep "<<pFieldMgr->GetDeltaOneStep()/mm <<"mm" <<endl;
 
      G4ChordFinder *pChordFinder = new G4ChordFinder(sr1Field,0.01*mm,pStepper);
      pFieldMgr->SetChordFinder( pChordFinder );
      pFieldMgr->SetDetectorField(sr1Field);
}



void sr1DetectorConstruction::ReportGeometryProblem(char myString[501]) {
  G4cout<<"\nE R R O R  in sr1DetectorConstruction.cc:    "
	<<"Unknown keyword requested in \""<< parameterFileName <<"\" :"<<G4endl;
  G4cout<<"\""<<myString<<"\""<<G4endl;
  G4cout<<"S T O P     F O R C E D!"<<G4endl;
  exit(1);
}

void sr1DetectorConstruction::ReportProblemInStearingFile(char* myString) {
  G4cout<<"\nE R R O R  in sr1DetectorConstruction.cc:    "
	<<"Problem to interpret/execute the following line in \""<< parameterFileName <<"\" :"<<G4endl;
  G4cout<<"\""<<myString<<"\""<<G4endl;
}


G4Material* sr1DetectorConstruction::CharToMaterial(char myString[100]) {
  G4Material* Material=NULL;
//!  if (strcmp(myString,"VacMater")==0) { Material = G4Material::GetMaterial("ArgonGas"); }
//!  else {
    G4String materialName = myString;
    G4NistManager*   man = G4NistManager::Instance();
    Material = man->FindOrBuildMaterial(materialName);   
    //Material=G4Material::GetMaterial(materialName);
//!  }
  if (Material==NULL) {
    G4cout<<"\nE R R O R  in sr1DetectorConstruction.cc::CharToMaterial    "
	<<"Unknown material requested:"<<G4endl;
    G4cout<<myString<<G4endl;
    G4cout<<"S T O P     F O R C E D!"<<G4endl;
    exit(1);
  }
  return Material;
}


G4LogicalVolume* sr1DetectorConstruction::FindLogicalVolume(G4String LogicalVolumeName) {
  G4LogicalVolumeStore* pLogStore = G4LogicalVolumeStore::GetInstance();
  if (pLogStore==NULL) {
    G4cout<<"ERROR:  sr1DetectorConstruction.cc:    G4LogicalVolumeStore::GetInstance()  not found!"<<G4endl;
  }
  else {
   //	  for (G4int i=0; i<pLogStore->entries(); i++) {
    for (unsigned int i=0; i<pLogStore->size(); i++) {
      G4LogicalVolume* pLogVol=pLogStore->at(i);
      G4String iLogName=pLogVol->GetName(); 
      if (iLogName==LogicalVolumeName) {
	return pLogVol;
      }
    }
  }
  return NULL;
}

void sr1DetectorConstruction::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: sr1DetectorConstruction::SetColourOfLogicalVolume:   unknown colour requested: "<<colour<<G4endl;
      G4cout<<"    Ignoring and continuing"<<G4endl;
    }
  }
}



G4bool sr1DetectorConstruction::CheckIfStringContainsSubstring(char* string, char* substring) {
  int nString=strlen(string);
  int nSubString=strlen(substring);
  for (int i=0; i<(nString-nSubString+1); i++) {
    if ( memcmp( &string[i], substring, nSubString)==0 ) {
      return true;
    }
  }
  return false;
}