Problem 1085

Summary: G4QuadrupoleMagField
Product: Geant4 Reporter: b.riese
Component: geometry/magneticfieldAssignee: John Apostolakis <John.Apostolakis>
Status: RESOLVED FIXED    
Severity: normal CC: b.riese
Priority: P3    
Version: 9.2   
Hardware: PC   
OS: Linux   

Description b.riese 2009-10-13 16:26:22 CEST
Hi,
in my G4 simulation I use quadrupoles that are not in the origin of my coordinate system (0,0,0). 
G4 uses the global coordinate system with these formulas:

     Bx = B[0] = fGradient*X ,
     By = B[1] = fGradient*Y ,
     Bz = B[2] = 0 .

Hence there's just a deflection no focussing or defocussing effect.

How can I fix this?

Thanks for your help,
Björn Riese.
Comment 1 b.riese 2009-10-19 14:49:09 CEST
Here in addition the definition of my quadrupole field:

QuadField = new G4QuadrupoleMagField(fGradient);
       		fEquation = new G4Mag_UsualEqRhs(QuadField);
		fieldMgr = new G4FieldManager();
		fieldMgr->SetMinimumEpsilonStep( minEps );
		fieldMgr->SetMaximumEpsilonStep( maxEps );
		fieldMgr->SetDeltaOneStep(valueD1step); 
		fieldMgr->SetDeltaIntersection(valueDintersection); 
		fieldMgr->SetDetectorField(QuadField);

      		fStepper = new G4ClassicalRK4(fEquation);
      		fChordFinder = new G4ChordFinder(QuadField,fMinStep,fStepper);
		fieldMgr->SetChordFinder( fChordFinder );
		
       		Ptr->SetFieldManager(fieldMgr,true);

Thanks and best regards,
Björn Riese.
Comment 2 John Apostolakis 2009-10-21 17:42:54 CEST
Indeed the Quadrople class does not cater for being positioned at a different location.

It could be revised to add this as an option, and it will be a simple and useful extension.

For an immediate solution I have written to you, suggesting some code how to do this.  It is a simple extension of this class.

As this is a proposed revision or new requirement, I will close this problem report.

Best regards,
John Apostolakis
Comment 3 b.riese 2009-11-30 16:06:59 CET
Hi,
if you're using the methods
G4VPhysicalVolume->GetObjectTranslation();
G4VPhysicalVolume->GetObjectRotation();
as parameter value the quadrupole classes should work also in rotated volumes.

Best regards,
Björn Riese.

------------------- "G4QuadrupoleMagField.hh": --------------------------

#ifndef G4QUADRUPOLEMAGFIELD_HH
#define G4QUADRUPOLEMAGFIELD_HH

#include "G4MagneticField.hh"
#include "G4ThreeVector.hh"
#include "G4RotationMatrix.hh"

class G4QuadrupoleMagField : public G4MagneticField
{
  public: // with description

    G4QuadrupoleMagField(G4double pGradient, G4ThreeVector     
                   pOrigin, G4RotationMatrix* pMatrix);
    G4QuadrupoleMagField(G4double pGradient);
   ~G4QuadrupoleMagField();

    void GetFieldValue(const G4double yTrack[],
                             G4double B[]     ) const;
  private:

    G4double fGradient;
    G4ThreeVector fOrigin;
    G4RotationMatrix* fMatrix;
};

#endif


------------------- "G4QuadrupoleMagField.cc": --------------------------

#include "G4QuadrupoleMagField.hh"

G4QuadrupoleMagField::G4QuadrupoleMagField(G4double pGradient)
{
   fGradient = pGradient ;
}

/////////////////////////////////////////////////////////////////////////

G4QuadrupoleMagField::G4QuadrupoleMagField(G4double pGradient, G4ThreeVector pOrigin, G4RotationMatrix* pMatrix)
{
   fGradient    = pGradient ;
   fOrigin      = pOrigin ;
   fMatrix      = pMatrix ;
}

/////////////////////////////////////////////////////////////////////////

G4QuadrupoleMagField::~G4QuadrupoleMagField()
{
}

////////////////////////////////////////////////////////////////////////

void G4QuadrupoleMagField::GetFieldValue( const G4double y[7],
                                                G4double B[3]  ) const  
{

   G4ThreeVector r_global = G4ThreeVector(
		y[0] - fOrigin.x(), 
		y[1] - fOrigin.y(), 
		y[2] - fOrigin.z());

   G4ThreeVector r_local = G4ThreeVector(
   		fMatrix->colX() * r_global,
   		fMatrix->colY() * r_global,
   		fMatrix->colZ() * r_global);

   G4ThreeVector B_local = G4ThreeVector(
   		fGradient * r_local.y(),
		fGradient * r_local.x(),
		0);

   G4ThreeVector B_global = G4ThreeVector(
   		fMatrix->inverse().rowX() * B_local,
   		fMatrix->inverse().rowY() * B_local,
   		fMatrix->inverse().rowZ() * B_local);

   B[0] = B_global.x() ;
   B[1] = B_global.y() ;
   B[2] = B_global.z() ;

}
Comment 4 John Apostolakis 2012-12-06 14:22:42 CET
The improvement adding an origin is incorporated in Geant4 release 9.6  (November 30th, 2012.)