From 751641e1668041097d3a4f82a41f9f23635ab15e Mon Sep 17 00:00:00 2001 From: Kamil Sedlak Date: Thu, 30 Sep 2010 12:24:55 +0000 Subject: [PATCH] 30.9.2010 Kamil Sedlak - added special geometry of GPDmHolder for the m-counter light guide at GPD - added specail case of 3DBQuadVrankovic field map for the quadrupoles (format defined by Vjeran Vrankovic) --- src/musrDetectorConstruction.cc | 15 ++++++++++++ src/musrTabulatedElementField.cc | 39 ++++++++++++++++++++++++++++---- 2 files changed, 49 insertions(+), 5 deletions(-) diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 1455ae8..d866237 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -397,6 +397,21 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() { G4ThreeVector zTransA5(0,0,111.25*mm); solid = new G4SubtractionSolid(solidName, solidA123, solidGPDBoxA5, yRotA12, zTransA5); } + else if (strcmp(tmpString2,"GPDmHolder")==0){ + // Light guide that holds m0 in its position + 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",sensitiveDet,&volumeID); + solidName+=name; + G4Box* solidGPDcutOut = new G4Box ("solidGPDcutOut",x1*mm,x1*mm,(x3+0.01)*mm); + G4Box* solidGPDmHolder = new G4Box ("solidGPDmHolder",sqrt(2)*x1*mm,x2*mm,x3*mm); + // G4RotationMatrix* rot = new G4RotationMatrix(45*deg,0,0); + G4RotationMatrix* rot = new G4RotationMatrix(G4ThreeVector(0,0,1),45*deg); + // G4RotationMatrix* rot = new G4RotationMatrix(); + G4ThreeVector trans(0,x2*mm,0); + solid = new G4SubtractionSolid(solidName,solidGPDmHolder,solidGPDcutOut,rot,trans); + // solid = new G4SubtractionSolid(solidName,solidGPDcutOut,solidGPDmHolder,rot,trans); + } else ReportGeometryProblem(line); diff --git a/src/musrTabulatedElementField.cc b/src/musrTabulatedElementField.cc index 157ed4c..cee4a21 100644 --- a/src/musrTabulatedElementField.cc +++ b/src/musrTabulatedElementField.cc @@ -38,6 +38,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons // 2DBOperaXY = 2D_OperaXY .... 2D field, magnetic, Opera format with (Kamil-like), X and Y components // 2DE ... 2D field , electric, Toni-like (2DE WAS TESTED) // 2DB ... 2D field , magnetic, Toni-like + // 3DBQuadVrankovic ... 3D field of a quadrupole magnet asymmetric in z but with the symmetric octants in (x,y) plane G4double lenUnit = 1*m; // length unit of the field map grid coordinates G4double fieldNormalisation = 1.; // Normalisation factor by which the the field map has to be multiplied @@ -110,13 +111,19 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons file.ignore(256, '\n'); } while (file.peek() == '%'); } - else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)||(strcmp(fieldTableType,"3DEOpera")==0)) { + else if ((strcmp(fieldTableType,"3D")==0)||(strcmp(fieldTableType,"3DBOpera")==0)||(strcmp(fieldTableType,"3DEOpera")==0) + ||(strcmp(fieldTableType,"3DBQuadVrankovic")==0)) { // OPERA format of the input file: // Read table dimensions lenUnit = 1*m; fieldNormalisation = 1.; - G4cout << "3D, field-map file format from OPERA (Kamil)" << G4endl; - G4cout << " ---> Assumed order (7 col.): x, y, z, "< Assumed order (7 col.): x, y, z, "<> xval >> yval >> zval >> bx >> by >> bz >> permeability; // G4cout<< xval <<" "<< yval <<" "<< zval <<" "<< bx <<" "<< by <<" "<< bz <> xval >> yval >> zval >> bx >> by >> bz; + } else if ((strcmp(fieldTableType,"2DE")==0)||(strcmp(fieldTableType,"2DB")==0)||(strcmp(fieldTableType,"2DEf")==0)) { file >> xval >> zval >> bx >> bz; } @@ -336,7 +346,7 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons file.close(); // G4cout<<"."< UNDERSTAND THE PROBLEM BEFORE PROCEEDING FURTHER!"<0) {maximumWidth = 2*dx; maximumHeight = 2*dy; maximumLength = 2*dz;} + if (symmetryType==101) {maximumWidth = 2*dx; maximumHeight = 2*dx; maximumLength = 2*dz;} // "maximumHeight = 2*dx" is correct here! + // maximumLength = 2*dz ==> 1*dz would be fine for symmetric (in z) case, but 2*dz is safe for z going eg. from -80 to 20 cm. + else if (symmetryType>0) {maximumWidth = 2*dx; maximumHeight = 2*dy; maximumLength = 2*dz;} else {maximumWidth = dx; maximumHeight = dy; maximumLength = dz;} } } @@ -441,6 +453,8 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4], double y_sign = 1.; double z_sign = 1.; + G4bool xySwap=false; // used only for Quadrupoles, which have additional symmetry line at 45 degrees. + if (symmetryType!=0) { x_sign = (x>0) ? 1.:-1.; y_sign = (y>0) ? 1.:-1.; @@ -450,6 +464,10 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4], case (2): x=fabs(x); y=fabs(y); z=fabs(z); break; + case (101): + x=fabs(x); y=fabs(y); + if (y>x) {xySwap=true; double xtmp=x; x=y; y=xtmp;} + break; default : G4cout<<" musrTabulatedElementField: unsupported symmetryType (="<(ydindex); int zindex = static_cast(zdindex); + // The following is not neede (it should be treated by the code several lines above). + // if (symmetryType==101) { + // if (yindex>(ny-2)) return; // field not defined in this region for the quadrupole + // } + //cks The following check is necessary - even though xindex and zindex should never be out of range, // it may happen (due to some rounding error ?). It is better to leave the check here. if ((xindex<0)||(xindex>(nx-2))) { @@ -552,6 +575,12 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4], case (2): B[0]*=x_sign; B[2]*=z_sign; break; + case (101): + if (xySwap) { + double B0tmp= B[0]; B[0]=B[1]; B[1]=B0tmp; + } + B[0]*=y_sign; B[1]*=x_sign; + break; default : G4cout<<" musrTabulatedElementField: unsupported symmetryType (="<