Merge pull request #145 from gemc-eic/master

add cartesian_3d field support

SConstruct

@@ -93,7 +93,8 @@ field_sources = Split("""
 	fields/multipoleField.cc
 	fields/symmetries/dipole.cc
 	fields/symmetries/cylindrical.cc
-	fields/symmetries/phi-segmented.cc""")
+	fields/symmetries/phi-segmented.cc
+	fields/symmetries/cartesian_3d.cc""")
 env.Library(source = field_sources, target = "lib/gfields")
 
 

fields/asciiField.cc

@@ -192,15 +192,20 @@ void asciiField::loadFieldMap(gMappedField* map, double v)
 	// dipole field
 	if(map->symmetry == "dipole-x" || map->symmetry == "dipole-y" || map->symmetry == "dipole-z")
 		loadFieldMap_Dipole(map, v);
-
+	
 	// cylindrical field
-	if(map->symmetry == "cylindrical-x" || map->symmetry == "cylindrical-y" || map->symmetry == "cylindrical-z")
+	else if(map->symmetry == "cylindrical-x" || map->symmetry == "cylindrical-y" || map->symmetry == "cylindrical-z")
 		loadFieldMap_Cylindrical(map, v);
 
 	// phi-segmented field
 
-	if(map->symmetry == "phi-segmented")
+	else if(map->symmetry == "phi-segmented")
 		loadFieldMap_phiSegmented(map, v);
+
+	else if(map->symmetry == "cartesian_3D" || map->symmetry == "cartesian_3D_quadrant")
+		loadFieldMap_cartesian3d(map, v);	
+
+	else {cout << "can't recognize the field symmetry "<< map->symmetry << endl; exit(0);}
 }
 
 

fields/asciiField.h

@@ -19,7 +19,8 @@ class asciiField : public fieldFactory
 
 	void loadFieldMap_Dipole(gMappedField*, double);       // load dipole field map
 	void loadFieldMap_Cylindrical(gMappedField*, double);  // load cylindrical field map
-	void loadFieldMap_phiSegmented(gMappedField*, double); // load cylindrical field map
+	void loadFieldMap_phiSegmented(gMappedField*, double); // load phiSegmented field map
+	void loadFieldMap_cartesian3d(gMappedField*, double); // load cartesian3d field map	
 
 	static fieldFactory *createFieldFactory()
 	{

fields/fieldFactory.h

@@ -20,7 +20,9 @@ class fieldFactory
 		virtual void loadFieldMap(gMappedField*, double)             = 0; // load field map. This is a dispatcher function for the various types of fields below
 
 		virtual void loadFieldMap_Dipole(gMappedField*, double)      = 0; // load 1D dipole field depending on 2 coordinates (transverse and longitudinal)
-		virtual void loadFieldMap_Cylindrical(gMappedField*, double) = 0; // load phi-symmetric field depending on 2 coordinates (transverse and longitudinal)
+		virtual void loadFieldMap_Cylindrical(gMappedField*, double) = 0; // load cylindrical field depending on 2 coordinates (transverse and longitudinal)
+		virtual void loadFieldMap_phiSegmented(gMappedField*, double) = 0; // load phi-symmetric field depending on 3 coordinates (azimuthal,transverse and longitudinal)
+		virtual void loadFieldMap_cartesian3d(gMappedField*, double) = 0; // load cartesian_3D field depending on 3 coordinates (XX,YY,ZZ)				
 
 		string factoryType;
 

fields/mappedField.cc

@@ -46,7 +46,12 @@ void gMappedField::GetFieldValue(const double point[3], double *bField) const
 	// phi-segmented
 	} else 	if(symmetry == "phi-segmented") {
 		GetFieldValue_phiSegmented(rpoint, bField, FIRST_ONLY);
+	// cartesian_3D
+	} else 	if(symmetry == "cartesian_3D" || symmetry == "cartesian_3D_quadrant"){
+		GetFieldValue_cartesian3d(rpoint, bField, FIRST_ONLY);
 	}
+
+
 
 	if(verbosity == 99) FIRST_ONLY = 99;
 
@@ -89,13 +94,16 @@ void gMappedField::initializeMap()
 	if(symmetry == "dipole-x" || symmetry == "dipole-y" || symmetry == "dipole-z")
 	{
 		startMap = new double[2];
+		endMap = new double[2];		
 		cellSize = new double[2];
 		np       = new unsigned int[2];
 
 		np[0]       = getCoordinateWithName("longitudinal").np;
 		np[1]       = getCoordinateWithName("transverse").np;
 		startMap[0] = getCoordinateWithName("longitudinal").min;
 		startMap[1] = getCoordinateWithName("transverse").min;
+		endMap[0] = getCoordinateWithName("longitudinal").max;
+		endMap[1] = getCoordinateWithName("transverse").max;		
 		cellSize[0] = (getCoordinateWithName("longitudinal").max - startMap[0]) / (np[0] - 1);
 		cellSize[1] = (getCoordinateWithName("transverse").max   - startMap[1]) / (np[1] - 1);
 	}
@@ -105,13 +113,16 @@ void gMappedField::initializeMap()
 	if(symmetry == "cylindrical-x" || symmetry == "cylindrical-y" || symmetry == "cylindrical-z")
 	{
 		startMap = new double[2];
+		endMap = new double[2];				
 		cellSize = new double[2];
 		np       = new unsigned int[2];
 
 		np[0]       = getCoordinateWithName("transverse").np;
 		np[1]       = getCoordinateWithName("longitudinal").np;
 		startMap[0] = getCoordinateWithName("transverse").min;
 		startMap[1] = getCoordinateWithName("longitudinal").min;
+		endMap[0] = getCoordinateWithName("longitudinal").max;
+		endMap[1] = getCoordinateWithName("transverse").max;				
 		cellSize[0] = (getCoordinateWithName("transverse").max   - startMap[0]) / (np[0] - 1);
 		cellSize[1] = (getCoordinateWithName("longitudinal").max - startMap[1]) / (np[1] - 1);
 	}
@@ -121,6 +132,7 @@ void gMappedField::initializeMap()
 	if(symmetry == "phi-segmented")
 	{
 		startMap = new double[3];
+		endMap = new double[3];				
 		cellSize = new double[3];
 		np       = new unsigned int[3];
 
@@ -130,11 +142,36 @@ void gMappedField::initializeMap()
 		startMap[0] = getCoordinateWithName("azimuthal").min;
 		startMap[1] = getCoordinateWithName("transverse").min;
 		startMap[2] = getCoordinateWithName("longitudinal").min;
+		endMap[0] = getCoordinateWithName("azimuthal").max;
+		endMap[1] = getCoordinateWithName("transverse").max;				
+		endMap[2] = getCoordinateWithName("longitudinal").max;		
 		cellSize[0] = (getCoordinateWithName("azimuthal").max    - startMap[0]) / (np[0] - 1);
 		cellSize[1] = (getCoordinateWithName("transverse").max   - startMap[1]) / (np[1] - 1);
 		cellSize[2] = (getCoordinateWithName("longitudinal").max - startMap[2]) / (np[2] - 1);
 	}
 
+	//cartesian_3D
+	if(symmetry == "cartesian_3D" || symmetry == "cartesian_3D_quadrant")
+	{
+		startMap = new double[3];
+		endMap = new double[3];				
+		cellSize = new double[3];
+		np       = new unsigned int[3];
+
+		np[0]       = getCoordinateWithName("X").np;
+		np[1]       = getCoordinateWithName("Y").np;
+		np[2]       = getCoordinateWithName("Z").np;
+		startMap[0] = getCoordinateWithName("X").min;
+		startMap[1] = getCoordinateWithName("Y").min;
+		startMap[2] = getCoordinateWithName("Z").min;
+		endMap[0] = getCoordinateWithName("X").max;
+		endMap[1] = getCoordinateWithName("Y").max;				
+		endMap[2] = getCoordinateWithName("Z").max;				
+		cellSize[0] = (getCoordinateWithName("X").max    - startMap[0]) / (np[0] - 1);
+		cellSize[1] = (getCoordinateWithName("Y").max   - startMap[1]) / (np[1] - 1);
+		cellSize[2] = (getCoordinateWithName("Z").max - startMap[2]) / (np[2] - 1);
+	}	
+
 	// setting rotation sin and cosines
 	sinAlpha = sin(mapRotation[0]);
 	cosAlhpa = cos(mapRotation[0]);

fields/mappedField.h

@@ -119,6 +119,7 @@ class gMappedField : public G4MagneticField
 	// it will avoid the time spend in retrieving
 	// the infos in GetFieldValue
 	double *startMap;
+	double *endMap;	
 	double *cellSize;
 	unsigned int *np;
 	void initializeMap();
@@ -132,6 +133,7 @@ class gMappedField : public G4MagneticField
 	void GetFieldValue_Dipole( const double x[3], double *Bfield, int FIRST_ONLY) const;
 	void GetFieldValue_Cylindrical( const double x[3], double *Bfield, int FIRST_ONLY) const;
 	void GetFieldValue_phiSegmented( const double x[3], double *Bfield, int FIRST_ONLY) const;
+	void GetFieldValue_cartesian3d( const double x[3], double *Bfield, int FIRST_ONLY) const;
 
 
 	// we want to rotate the field (axes), not the point