diff --git a/SConstruct b/SConstruct index 7389c281..c1b8d09f 100644 --- a/SConstruct +++ b/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") diff --git a/fields/asciiField.cc b/fields/asciiField.cc index 12a48832..6ca4fcfb 100644 --- a/fields/asciiField.cc +++ b/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);} } diff --git a/fields/asciiField.h b/fields/asciiField.h index 75cf5bbf..c02dd82f 100644 --- a/fields/asciiField.h +++ b/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() { diff --git a/fields/fieldFactory.h b/fields/fieldFactory.h index 3f1bc841..73cd2e47 100644 --- a/fields/fieldFactory.h +++ b/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; diff --git a/fields/mappedField.cc b/fields/mappedField.cc index 6a55b6fd..ac0bcf26 100644 --- a/fields/mappedField.cc +++ b/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,6 +94,7 @@ 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]; @@ -96,6 +102,8 @@ void gMappedField::initializeMap() 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,6 +113,7 @@ 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]; @@ -112,6 +121,8 @@ void gMappedField::initializeMap() 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]); diff --git a/fields/mappedField.h b/fields/mappedField.h index af0c2d89..87c65ea3 100644 --- a/fields/mappedField.h +++ b/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 diff --git a/fields/symmetries/cartesian_3d.cc b/fields/symmetries/cartesian_3d.cc new file mode 100644 index 00000000..44f9e816 --- /dev/null +++ b/fields/symmetries/cartesian_3d.cc @@ -0,0 +1,285 @@ +// gemc headers +#include "asciiField.h" +#include "utils.h" + + +// 3D field in cartesian coordinates. Field itself is in cartesian coordinates. +// Dependent on 3 cartesian coordinates (Y, X, Z) +// The values can be loaded from the map in any order as long as their speed is ordered. +// The values are indexed as B1_3D[X][Y][Z],B2_3D[X][Y][Z],B3_3D[X][Y][Z] +// The field is three dimensional, ordered in the class as B1=Bx, B2=By, B3=Bz + +//symmetry "cartesian_3D" is for full 3D map +//symmetry "cartesian_3D_quadrant" is for 3D map covering only the 1st quadrant where x>=0 && y>=0 + +// from fieldFactory: load field map +void asciiField::loadFieldMap_cartesian3d(gMappedField* map, double verbosity) +{ + setlocale(LC_NUMERIC, "en_US"); + + int np_1 = map->getCoordinateWithSpeed(0).np; + int np_2 = map->getCoordinateWithSpeed(1).np; + int np_3 = map->getCoordinateWithSpeed(2).np; + double min1 = map->getCoordinateWithSpeed(0).min; + double min2 = map->getCoordinateWithSpeed(1).min; + double min3 = map->getCoordinateWithSpeed(2).min; + double cell1 = (map->getCoordinateWithSpeed(0).max - min1)/(np_1 - 1); + double cell2 = (map->getCoordinateWithSpeed(1).max - min2)/(np_2 - 1); + double cell3 = (map->getCoordinateWithSpeed(2).max - min3)/(np_3 - 1); + + // Allocate memory. [AZI][TRANSVERSE][LONGI] + // as initialized in the map + map->B1_3D = new double**[map->np[0]]; + map->B2_3D = new double**[map->np[0]]; + map->B3_3D = new double**[map->np[0]]; + for (unsigned i = 0; i < map->np[0]; ++i) + { + map->B1_3D[i] = new double*[map->np[1]]; + map->B2_3D[i] = new double*[map->np[1]]; + map->B3_3D[i] = new double*[map->np[1]]; + for (unsigned j = 0; j < map->np[1]; ++j) + { + map->B1_3D[i][j] = new double[map->np[2]]; + map->B2_3D[i][j] = new double[map->np[2]]; + map->B3_3D[i][j] = new double[map->np[2]]; + + } + } + + double unit1 = get_number("1*" + map->getCoordinateWithSpeed(0).unit); + double unit2 = get_number("1*" + map->getCoordinateWithSpeed(1).unit); + double unit3 = get_number("1*" + map->getCoordinateWithSpeed(2).unit); + + double scale = map->scaleFactor*get_number("1*" + map->unit); + + double d1, d2, d3, b1, b2, b3; + + // progress bar + int barWidth = 50; + + // using fscanf instead of c++ to read file, this is a lot faster + char ctmp[100]; + string tmp; + FILE *fp = fopen (map->identifier.c_str(), "r"); + + // ignoring header + while(tmp != "") + { + if(fscanf(fp, "%s", ctmp) != 0) + tmp = string(ctmp); + } + + // now reading map + // values as read from map + for(int i1 = 0; i1 0.001) { + cout << " !! Error: coordinate index wrong. Map point should be " << min1 + i1*cell1 + << " but it's " << d1 << " instead." << endl; + } + // checking map consistency for second coordinate + if( (min2 + i2*cell2 - d2)/d2 > 0.001) { + cout << " !! Error: coordinate index wrong. Map point should be " << min2 + i2*cell2 + << " but it's " << d2 << " instead." << endl; + } + + // checking map consistency for third coordinate + if( (min3 + i3*cell3 - d3)/d3 > 0.001) { + cout << " !! Error: coordinate index wrong. Map point should be " << min2 + i2*cell2 + << " but it's " << d2 << " instead." << endl; + } + + + // calculating index + // this + unsigned t1 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + unsigned t2 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + unsigned t3 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + + + // The values are indexed as B1_3D[X][Y][Z] + if( map->getCoordinateWithSpeed(0).name == "X" + && map->getCoordinateWithSpeed(1).name == "Y" + && map->getCoordinateWithSpeed(2).name == "Z" ) { + map->B1_3D[t1][t2][t3] = b1; + map->B2_3D[t1][t2][t3] = b2; + map->B3_3D[t1][t2][t3] = b3; + } else if( map->getCoordinateWithSpeed(0).name == "X" + && map->getCoordinateWithSpeed(1).name == "Z" + && map->getCoordinateWithSpeed(2).name == "Y" ) { + t1 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + t2 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + t3 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + + map->B1_3D[t1][t3][t2] = b1; + map->B2_3D[t1][t3][t2] = b2; + map->B3_3D[t1][t3][t2] = b3; + } else if( map->getCoordinateWithSpeed(0).name == "Z" + && map->getCoordinateWithSpeed(1).name == "X" + && map->getCoordinateWithSpeed(2).name == "Y" ) { + t1 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + t2 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + t3 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + + map->B1_3D[t3][t1][t2] = b1; + map->B2_3D[t3][t1][t2] = b2; + map->B3_3D[t3][t1][t2] = b3; + } else if( map->getCoordinateWithSpeed(0).name == "Z" + && map->getCoordinateWithSpeed(1).name == "Y" + && map->getCoordinateWithSpeed(2).name == "X" ) { + t1 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + t2 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + t3 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + + map->B1_3D[t3][t2][t1] = b1; + map->B2_3D[t3][t2][t1] = b2; + map->B3_3D[t3][t2][t1] = b3; + } else if( map->getCoordinateWithSpeed(0).name == "Y" + && map->getCoordinateWithSpeed(1).name == "Z" + && map->getCoordinateWithSpeed(2).name == "X" ) { + t1 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + t2 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + t3 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + + map->B1_3D[t2][t3][t1] = b1; + map->B2_3D[t2][t3][t1] = b2; + map->B3_3D[t2][t3][t1] = b3; + } else if( map->getCoordinateWithSpeed(0).name == "Y" + && map->getCoordinateWithSpeed(1).name == "X" + && map->getCoordinateWithSpeed(2).name == "Z" ) { + t1 = (unsigned) floor( ( d2 - min2 + cell2/2 ) / ( cell2 ) ) ; + t2 = (unsigned) floor( ( d1 - min1 + cell1/2 ) / ( cell1 ) ) ; + t3 = (unsigned) floor( ( d3 - min3 + cell3/2 ) / ( cell3 ) ) ; + + map->B1_3D[t2][t1][t3] = b1; + map->B2_3D[t2][t1][t3] = b2; + map->B3_3D[t2][t1][t3] = b3; + } + + if(verbosity>4 && verbosity != 99) { + cout << " Loading Map: coordinates (" << d1 << ", " << d2 << ", " << d3 << ") values: (" << b1 << ", " << b2 << ", " << b3 << ")"; + cout << ", indexes (t1, t2, t3) = (" << t1 << ", " << t2 << ", " << t3 << ") " ; + cout << ", array sizes = (" << np_1 << ", " << np_2 << ", " << np_3 << ") " << endl; + } + + + } + } + } + + cout << " ["; + double progress = (double)i1/(double)np_1; + int pos = progress*barWidth; + for (int i = 0; i < barWidth; ++i) + { + if (i < pos) cout << "="; + else if (i == pos) cout << ">"; + else cout << " "; + } + cout << "] " << int(progress * 100.0) << " %\r"; + cout.flush(); + } + fclose(fp); + + cout << endl; +} + +// from mappedField: GetFieldValue +void gMappedField::GetFieldValue_cartesian3d( const double x[3], double *Bfield, int FIRST_ONLY) const +{ + + double xx = x[0]; + double yy = x[1]; + double zz = x[2]; + + double XX = 0; + double YY = 0; + double ZZ = 0; + + if(symmetry == "cartesian_3D"){ + XX = xx; YY = yy; ZZ = zz; + }else if (symmetry == "cartesian_3D_quadrant"){ + if (xx>=0 && yy>=0) { XX = xx; YY = yy; ZZ = zz;} + if (xx>=0 && yy<0) { XX = -yy; YY = xx; ZZ = zz;} + if (xx<0 && yy<0) { XX = -xx; YY = -yy; ZZ = zz;} + if (xx<0 && yy>=0) { XX = yy; YY = -xx; ZZ = zz;} + + if (XX<0 || YY <0) {cout << "xx " << xx << " yy " << yy << " zz " << zz << " XX " << XX << " YY " << YY << " ZZ " << ZZ << endl; return;} + } + + // map indexes, bottom of the cell + // 0 to np-1 + unsigned int IXX = floor( ( XX - startMap[0] ) / cellSize[0] ); + unsigned int IYY = floor( ( YY - startMap[1] ) / cellSize[1] ); + unsigned int IZZ = floor( ( ZZ - startMap[2] ) / cellSize[2] ); + + // outside map, returning no field + if (XX < startMap[0] || YY < startMap[1] || ZZ < startMap[2]) return; + if (XX >= endMap[0] || YY >= endMap[1] || ZZ >= endMap[2]) return; + + // no interpolation + if(interpolation == "none") + { + // checking if the point is closer to the top of the cell + if( fabs( startMap[0] + IXX*cellSize[0] - XX) > fabs( startMap[0] + (IXX+1)*cellSize[0] - XX) ) IXX++; + if( fabs( startMap[0] + IYY*cellSize[0] - YY) > fabs( startMap[0] + (IYY+1)*cellSize[0] - YY) ) IYY++; + if( fabs( startMap[0] + IZZ*cellSize[0] - ZZ) > fabs( startMap[0] + (IZZ+1)*cellSize[0] - ZZ) ) IZZ++; + + Bfield[0] = B1_3D[IXX][IYY][IZZ]; + Bfield[1] = B2_3D[IXX][IYY][IZZ]; + Bfield[2] = B3_3D[IXX][IYY][IZZ]; + } + else if (interpolation == "linear") + { + // relative positions within cell + double XXR = (XX - (startMap[0] + IXX*cellSize[0])) / cellSize[0]; + double YYR = (YY - (startMap[1] + IYY*cellSize[1])) / cellSize[1]; + double ZZR = (ZZ - (startMap[2] + IZZ*cellSize[2])) / cellSize[2]; + + // field component interpolation + double B1 = B1_3D[IXX][IYY][IZZ] * (1.0 - XXR) + B1_3D[IXX][IYY][IZZ] * XXR; + double B2 = B2_3D[IXX][IYY][IZZ] * (1.0 - YYR) + B2_3D[IXX][IYY][IZZ] * YYR; + double B3 = B3_3D[IXX][IYY][IZZ] * (1.0 - ZZR) + B3_3D[IXX][IYY][IZZ] * ZZR; + + Bfield[0] = B1; + Bfield[1] = B2; + Bfield[2] = B3; + } + else + { + cout << " !! Unkown field interpolation method >" << interpolation << "<" << endl; + return; + } + + // we don't worry about computer speed + // if verbosity is set this high + // so we can output units as well + if(verbosity>3 && FIRST_ONLY != 99) + { + cout << " > Track position in magnetic field map, with displacement and rotations (x,y,z)/cm:" + << "(" << x[0]/cm << ", " + << x[1]/cm << ", " + << x[2]/cm << ") cm, " << endl; + cout << " 3D XYZ: "; + cout << "loc. pos. = (" << x[0]/cm << ", " << x[1]/cm << ", " << x[2]/cm << ") cm, "; + cout << "XX=" << XX/cm << "cm, YY=" << YY/cm << "cm, ZZ=" << ZZ/cm << ", "; + cout << "IXX=" << IXX << " "; + cout << "IYY=" << IYY << " "; + cout << "IZZ=" << IZZ << ", "; + cout << "B = (" << Bfield[0]/gauss << ", " << Bfield[1]/gauss << ", " << Bfield[2]/gauss << ") gauss " << endl; + } +}