3.12.2009 Kamil Sedlak

The reading of the field map files in the 3D Opera file format has been extended:
  1) The lenght unit can also be specified in cm instead of in meters only.
  2) Only one octant (in Kartesian coordinate system) of the field map can be
	specified, and the field can be extrapolated to other octants (this
	depends on the symmetry of the field).
  3) First line of the field map does not have to be empty (but it can be).
  4) Documentation has been modified accordingly.

doc/musrSim.tex

@@ -324,8 +324,29 @@ Three special volumes ``Target, M0, M1 and M2''.
 			the \emph{field normalisation factor} is 1. (Note that this default
                         normalisation is different from 2DBOpera and 2DBOperaXY options).
 			However, a different \emph{field normalisation factor} can be specified
-			in the field map using the keyword ``fieldNormalisation \emph{number}'' 
+			in the field map file using the keyword ``fieldNormalisation \emph{number}'' 
                         before the line started with 0.\\
+			The \emph{length unit} can be changed to 1\,cm by specifying ``[CENTIMETRE]''
+			after the ``0'' character in the field map file.\\
+			It is expected that the field map is defined in the full volume of the field.
+			Sometimes (due to the symmetry of the field), it is enough to define the
+			field in just one octant of the Kartesian coorinate system (e.g. for positive
+			$x$, $y$ and $z$).  In such cases, the user can specify this in the field map
+			file using the keyword ``symmetryType \emph{number}'', where the \emph{number}
+			specifies how the field should be extrapolated to other octants.  
+			The ``symmetryType 1'' case means that the planes of symmetry are (x,y) and
+			(x,z), i.e. if the field at point $(x,y,z)$ is $(F_x,F_y,F_z)$, the field in
+			different octants will look like this: 
+			$(-x,y,z) \rightarrow (F_x,F_y,F_z)$;
+			$(x,-y,z) \rightarrow (F_x,-F_y,F_z)$;
+			$(-x,-y,z) \rightarrow (F_x,-F_y,F_z)$;
+			$(x,y,-z) \rightarrow (F_x,F_y,-F_z)$;
+			$(-x,y,-z) \rightarrow (F_x,F_y,-F_z)$;
+			$(x,-y,-z) \rightarrow (F_x,-F_y,-F_z)$;
+			$(-x,-y,-z) \rightarrow (F_x,-F_y,-F_z)$. \\
+			Similar case is the ``symmetryType 2'', where the planes of symmetry are
+			(x,y) and (y,z).
+			These two symmetry types are realised in a the spin rotator oriented along the z axis.\\
 			Example of the beginning of the field map file:\\
 			2   2   55\\
 			1 X\\
@@ -336,6 +357,7 @@ Three special volumes ``Target, M0, M1 and M2''.
 			6 BZ\\
 			7 DUMMY\\
 			fieldNormalisation -22.5733634\\
+			symmetryType 2\\
 			0 [METRE]\\
 			-0.2  -0.2  -1.35  0.  0.     0. 0.\\
 			-0.2  -0.2  -1.30  0. -0.0002 0. 0.\\

include/musrTabulatedElementField.hh

@@ -83,6 +83,8 @@ private:
   double dx, dy, dz;
   double ffieldValue;
   double maximumWidth, maximumHeight, maximumLength;
+  int    symmetryType;  // this variable defines whether (and how) should be the field extended
+                        // in the case it is defined just in one octant of the Cartesian coordinates.
 
   void   Invert(const char* indexToInvert);
 

src/F04GlobalField.cc

@@ -320,7 +320,7 @@ void F04GlobalField::PrintFieldAtRequestedPoints() const {
     object->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);
-        printf ("   EM field value at (%.f, %.f, %.f) mm is:  B = (%0.10g, %0.10g, %0.10g) T,  E = (%0.10g, %0.10g, %0.10g) kV/mm\n",
+        printf ("  Field at (%.2f, %.2f, %.2f) mm is:  B = (%0.10g, %0.10g, %0.10g) T,  E = (%0.10g, %0.10g, %0.10g) kV/mm\n",
         point[0]/mm, point[1]/mm, point[2]/mm,
         Bfi[0]/tesla,         Bfi[1]/tesla,         Bfi[2]/tesla,
         Bfi[3]/(kilovolt/mm), Bfi[4]/(kilovolt/mm), Bfi[5]/(kilovolt/mm));

src/musrTabulatedElementField.cc

@@ -53,6 +53,8 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
     fieUnit= kilovolt/mm;
   }
   
+  symmetryType=0;
+
   fldDim = (fieldTableType[0]=='3') ? 3:2;  
   if (fldDim==2) ny=1;
 
@@ -67,7 +69,9 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
 
   char buffer[256];
   char tmpString1[100]="Unset";
+  char tmpString2[100]="Unset";
   float fvalue;
+  int   ivalue;
   G4bool boolMinimaAndMaximaDefinedInTheFile = false;
   if (fldType=='E') G4cout<<"      Electric field ";
   if (fldType=='B') G4cout<<"      Magnetic field ";
@@ -112,8 +116,13 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
     G4cout << "3D, field-map file format from OPERA (Kamil)" << G4endl;
     G4cout << " ---> Assumed order (7 col.): x, y, z, "<<fldType<<"x, "<<fldType<<"y, "<<fldType<<"z, Dummy"<<G4endl;
 
-    file.getline(buffer,256); // Skip the first empty line of the file
+    //  Skip the first few empty lines of the file (if any) and read in the nx, ny and nz of the field map.
-    file >> nx >> ny >> nz; // Note dodgy order
+    int tmp_nx=-1; int tmp_ny=-1; int tmp_nz=-1; int nVarReadIn =-1;
+    do {
+      file.getline(buffer,256);
+      nVarReadIn = sscanf(&buffer[0],"%d %d %d",&tmp_nx,&tmp_ny,&tmp_nz);
+      nx = tmp_nx; ny = tmp_ny; nz = tmp_nz;
+    } while ((nx<=0)||(ny<=0)||(nz<=0)||(nVarReadIn<3));
     
     // Ignore other header information    
     // The first line whose second character is '0' is considered to
@@ -123,8 +132,18 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
       sscanf(&buffer[0],"%s %g",tmpString1,&fvalue);
       if (strcmp(tmpString1,"fieldNormalisation")==0) {
 	fieldNormalisation = fvalue;
-	G4cout << "DEBUG:  musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
+	//	G4cout << "DEBUG:  musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
       }
+      else if (strcmp(tmpString1,"symmetryType")==0) {
+	sscanf(&buffer[0],"%s %d",tmpString1,&ivalue);
+	symmetryType = ivalue;
+	if (symmetryType!=0) {
+	  G4cout<<"      The field should be extrapolated to different octants \n";
+	  G4cout<<"          of the Cartesian coord. system according to symmetryType = "<<symmetryType<<G4endl;
+	}
+      }
+      sscanf(&buffer[0],"%s %s",tmpString1,tmpString2);
+      if (strcmp(tmpString2,"[CENTIMETRE]")==0) {lenUnit=1*cm;}
     } while ( buffer[1]!='0');
   }
   else if ((strcmp(fieldTableType,"2DE")==0)||(strcmp(fieldTableType,"2DB")==0)||(strcmp(fieldTableType,"2DEf")==0)) {
@@ -169,7 +188,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
       sscanf(&buffer[0],"%s %g",tmpString1,&fvalue);
       if (strcmp(tmpString1,"fieldNormalisation")==0) {
 	fieldNormalisation = fvalue;
-	G4cout << "DEBUG:  musrTabulatedElementField: fieldNormalisation set to "<<fieldNormalisation<<G4endl;
+	G4cout << "                       fieldNormalisation set to "<<fieldNormalisation<<G4endl;
       }
     } while ( buffer[1]!='0');
   }
@@ -228,7 +247,9 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
   // Read in the data
   double xval,yval,zval,bx,by,bz;
   double permeability; // Not used in this example.
+  //  G4cout<<"      ";
   for (ix=0; ix<nx; ix++) {
+    //    G4cout<<"#"; G4cout.flush();
     for (iy=0; iy<ny; iy++) {
       for (iz=0; iz<nz; iz++) {
 	if ((strcmp(fieldTableType,"3DE")==0)||(strcmp(fieldTableType,"3DB")==0)) {
@@ -282,6 +303,16 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
     }
   }
   file.close();
+  //  G4cout<<"."<<G4endl;
+  if (symmetryType!=0) {
+    if ((minimumx!=0)||(minimumy!=0)||(minimumz!=0)) {
+      G4cout << "      musrTabulatedElementField::musrTabulatedElementField: Symmetrical extrapolation of the field"<<G4endl;
+      G4cout << "               required (symmetryType="<<symmetryType<<"), but the field does not start at point (0,0,0)"<<G4endl;
+      G4cout << "               STRANGE -> UNDERSTAND THE PROBLEM BEFORE PROCEEDING FURTHER!"<<G4endl;
+      G4cout << "                 =====>   S T O P " << G4endl;
+      musrErrorMessage::GetInstance()->musrError(FATAL,"musrTabulatedElementField: non-zero symmetryType for field that does not start at (0,0,0)",false);
+    }
+  }
 
 
   if (!boolMinimaAndMaximaDefinedInTheFile) {
@@ -347,9 +378,8 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
     else maximumLength = dz;
   }
   else {
-    maximumWidth = dx;
+    if (symmetryType>0) {maximumWidth = 2*dx; maximumHeight = 2*dy; maximumLength = 2*dz;}
-    maximumHeight = dy;
+    else                {maximumWidth =   dx; maximumHeight =   dy; maximumLength =   dz;}
-    maximumLength = dz;
   }
 }
 
@@ -375,13 +405,38 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
   double y = local.y();
   double z = local.z();
 
-  //  G4cout<<"Global points= "<<point[0]<<", "<<point[1]<<", "<<point[2]<<",  Local point= "<<x<<", "<<y<<", "<<z<<G4endl;
+  double x_sign = 1.;
+  double y_sign = 1.;
+  double z_sign = 1.;
+  
+  if (symmetryType!=0) {
+    x_sign = (x>0) ? 1.:-1.;
+    y_sign = (y>0) ? 1.:-1.;
+    z_sign = (z>0) ? 1.:-1.;
+    switch(symmetryType) {
+    case (1):
+    case (2):
+      x=fabs(x); y=fabs(y); z=fabs(z);
+      break;
+    default :
+      G4cout<<" musrTabulatedElementField: unsupported symmetryType (="<<symmetryType<<")"<<G4endl;
+    }
+  }
+
+  // G4cout<<"Global points= "<<point[0]<<", "<<point[1]<<", "<<point[2]<<",  Local point= "<<x<<", "<<y<<", "<<z<<G4endl;
+
+  // 26.11.2009  K.S.
+  // Check whether any of the axis is symmetrically or asymmetrically extended to negative values:
+  //  if (boolSymmetricOrAntisymmetricX||boolSymmetricOrAntisymmetricY||boolSymmetricOrAntisymmetricZ) {
+  //    if ( boolSymmetricOrAntisymmetricX && (x<0) ) {x=-x; x_sign=-1;}
+  //    
+  //  }
 
 
   // Check that the point is within the defined region 
-  if ( x>minimumx && x<maximumx &&
+  if ( x>=minimumx && x<maximumx &&
-       y>minimumy && y<maximumy && 
+       y>=minimumy && y<maximumy && 
-       z>minimumz && z<maximumz ) {
+       z>=minimumz && z<maximumz ) {
 
     // Position of given point within region, normalized to the range
     // [0,1]
@@ -456,6 +511,20 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
     B[1] *= ffieldValue;
     B[2] *= ffieldValue;
 
+    // G4cout<<"Pole1 = " << B[0]<<", "<<B[1]<<", "<<B[2]<<G4endl;
+    if (symmetryType!=0) {
+      switch(symmetryType) {
+      case (1): 
+	B[1]*=y_sign;  B[2]*=z_sign;
+	break;
+      case (2): 
+	B[0]*=x_sign;  B[2]*=z_sign;
+	break;
+      default :
+	G4cout<<" musrTabulatedElementField: unsupported symmetryType (="<<symmetryType<<")"<<G4endl;
+      }
+    }
+
     G4ThreeVector finalField(B[0],B[1],B[2]);
     finalField = global2local.Inverse().TransformAxis(finalField);