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

@@ -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);
 

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;
+    if (strcmp(fieldTableType,"3DBQuadVrankovic")==0) {
-    G4cout << " ---> Assumed order (7 col.): x, y, z, "<<fldType<<"x, "<<fldType<<"y, "<<fldType<<"z, Dummy"<<G4endl;
+      G4cout<<"3D field of a quadrupole defined by Vjeran Vrankovic."<<G4endl;
+    }
+    else {
+      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;
+    }
 
     //  Skip the first few empty lines of the file (if any) and read in the nx, ny and nz of the field map.
     int tmp_nx=-1; int tmp_ny=-1; int tmp_nz=-1; int nVarReadIn =-1;
@@ -269,6 +276,9 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
 	  file >> xval >> yval >> zval >> bx >> by >> bz >> permeability;
 	  //	  G4cout<< xval <<" "<< yval <<" "<< zval <<" "<< bx <<" "<< by <<" "<< bz <<G4endl;
 	}
+	else if (strcmp(fieldTableType,"3DBQuadVrankovic")==0) {
+	  file >> 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<<"."<<G4endl;
   if (symmetryType!=0) {
-    if ((minimumx!=0)||(minimumy!=0)||(minimumz!=0)) {
+    if ((minimumx!=0)||(minimumy!=0)||( (minimumz!=0)&&(symmetryType!=101) )) {
       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;
@@ -410,7 +420,9 @@ musrTabulatedElementField::musrTabulatedElementField( const char* filename, cons
     else maximumLength = dz;
   }
   else {
-    if (symmetryType>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 (="<<symmetryType<<")"<<G4endl;
     }
@@ -492,6 +510,11 @@ void musrTabulatedElementField::addFieldValue3D(const G4double point[4],
     int yindex = static_cast<int>(ydindex);
     int zindex = static_cast<int>(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 (="<<symmetryType<<")"<<G4endl;
       }