ShipFieldMaker Class Reference

Creates various magnetic fields and assigns them to geometry regions. More...

#include <ShipFieldMaker.h>

Inheritance diagram for ShipFieldMaker:

Collaboration diagram for ShipFieldMaker:

Classes
struct	fieldInfo
	Structure to hold volume name, field name and field scaling factor. More...
struct	transformInfo
	Structure to hold transformation information. More...

Public Types
typedef std::map< TString, TVirtualMagField * >	SFMap
	Typedef for a TString-TVirtualMagField* map.
typedef std::vector< std::string >	stringVect
	Typedef of a vector of strings.

Public Member Functions
	ShipFieldMaker (Bool_t verbose=kFALSE)
	Constructor.
virtual	~ShipFieldMaker ()
	Destructor.
virtual void	Construct ()
	Set-up all local and regional fields and assign them to volumes.
void	readInputFile (const std::string &inputFile)
	Read an input file to define fields and associated volumes.
void	defineUniform (const TString &name, const TVector3 &BVector)
	Define a uniform field.
void	defineConstant (const TString &name, const TVector2 &xRange, const TVector2 &yRange, const TVector2 &zRange, const TVector3 &BVector)
	Define a constant field.
void	defineBell (const TString &name, Double_t BPeak, Double_t zMiddle, Int_t orient=1, Double_t tubeR=500.0, Double_t xy=0.0, Double_t z=0.0, Double_t L=0.0)
	Define a Bell field.
void	defineFieldMap (const TString &name, const TString &mapFileName, const TVector3 &localCentre=TVector3(0.0, 0.0, 0.0), const TVector3 &localAngles=TVector3(0.0, 0.0, 0.0), Bool_t useSymmetry=kFALSE)
void	defineFieldMapCopy (const TString &name, const TString &mapNameToCopy, const TVector3 &translation, const TVector3 &eulerAngles=TVector3(0.0, 0.0, 0.0))
	Define a copy of a field map with a coordinate translation and optional rotation.
void	defineComposite (const TString &name, const TString &field1Name, const TString &field2Name, const TString &field3Name="", const TString &field4Name="")
	Define a composite field from up to four fields.
void	defineComposite (const TString &name, std::vector< TString > fieldNames)
	Define a composite field from a vector of field names.
void	defineGlobalField (const TString &field1Name, const TString &field2Name="", const TString &field3Name="", const TString &field4Name="")
	Define a composite field from up to four fields.
void	defineGlobalField (std::vector< TString > fieldNames)
	Define the Global field using a vector of field names.
void	defineRegionField (const TString &volName, const TString &fieldName, Double_t scale=1.0)
	Define a regional (local + global) field and volume pairing.
void	defineLocalField (const TString &volName, const TString &fieldName, Double_t scale=1.0)
	Define a localised field and volume pairing.
ShipCompField *	getGlobalField () const
	Get the global magnetic field.
SFMap	getAllFields () const
	Get the map storing volume names and their associated local magnetic fields.
TVirtualMagField *	getVolumeField (const TString &volName) const
	Get the magnetic field for the given volume.
Bool_t	gotField (const TString &name) const
	Check if we have a field stored with the given name name.
TVirtualMagField *	getField (const TString &name) const
	Get the field stored with the given name name.
void	plotXYField (const TVector3 &xAxis, const TVector3 &yAxis, const std::string &plotFile) const
	Plot the magnetic field in x-y plane.
void	plotZXField (const TVector3 &zAxis, const TVector3 &xAxis, const std::string &plotFile) const
	Plot the magnetic field in z-x plane.
void	plotZYField (const TVector3 &zAxis, const TVector3 &yAxis, const std::string &plotFile) const
	Plot the magnetic field in z-y plane.
void	plotField (Int_t type, const TVector3 &xAxis, const TVector3 &yAxis, const std::string &plotFile) const
	Plot the magnetic field in "x-y" plane.
void	generateFieldMap (TString fileName, const float step=2.5, const float xRange=179, const float yRange=317, const float zRange=1515.5, const float zShift=-4996)
	ClassDef (ShipFieldMaker, 1)
	Generate fieldMap csv file in the given region.

Protected Member Functions
void	defineUniform (const stringVect &inputLine)
	Define a uniform field based on information from the inputLine.
void	defineConstant (const stringVect &inputLine)
	Define a constant field based on information from the inputLine.
void	defineBell (const stringVect &inputLine)
	Define a Bell field based on information from the inputLine.
void	defineFieldMap (const stringVect &inputLine, Bool_t useSymmetry=kFALSE)
	Define a field map based on information from the inputLine.
void	defineFieldMapCopy (const stringVect &inputLine)
	Define a (translated) copy of a field map based on information from the inputLine.
void	defineComposite (const stringVect &inputLine)
	Define a composite field based on information from the inputLine.
void	defineGlobalField (const stringVect &inputLine)
	Define the global field based on information from the inputLine.
void	defineRegionField (const stringVect &inputLine)
	Define a regional (local+global) field based on the info from the inputLine.
void	setAllRegionFields ()
void	defineLocalField (const stringVect &inputLine)
	Define a local field only based on information from the inputLine.
void	checkLocalFieldMap (TVirtualMagField *&localField, const TString &volName, Double_t scale)
	Check if we have a local field map and store the volume global transformation.
void	setAllLocalFields ()
void	getTransformation (const TString &volName, transformInfo &theInfo)
	Get the transformation matrix for the volume position and orientation.
void	findNode (TGeoVolume *aVolume, const TString &volName)
	Update the current geometry node pointer that matches the required volume name.

Private Member Functions
stringVect	splitString (std::string &theString, std::string &splitter) const
	Split a string.

Private Attributes
ShipCompField *	globalField_
	The global magnetic field.
SFMap	theFields_
	The map storing all created magnetic fields.
std::vector< fieldInfo >	regionInfo_
	Vector of fieldInfo for regional fields.
std::vector< fieldInfo >	localInfo_
	Vector of fieldInfo for local fields.
Bool_t	verbose_
	Verbose boolean.
Double_t	Tesla_
	Double converting Tesla to kiloGauss (for VMC/FairRoot B field units)
TGeoNode *	theNode_
	The current volume node: used for finding volume transformations.
Bool_t	gotNode_
	Boolean to specify if we have found the volume node we need.

Detailed Description

Creates various magnetic fields and assigns them to geometry regions.

The internal units here are cm for distances and Tesla for fields. Geant4 units for distance are mm, B fields = 0.001 megaVolt*ns/mm^2 (1 Tesla). VMC units require cm and kGauss (0.1 Tesla). Internally, use cm and Tesla, so keep distances unchanged but multiply B fields by 10 (1 Tesla = 10 kGauss)

Author

John Back J.J.B.nosp@m.ack@.nosp@m.warwi.nosp@m.ck.a.nosp@m.c.uk

John Back J.J.B.nosp@m.ack@.nosp@m.warwi.nosp@m.ck.a.nosp@m.c.uk

This inherits from TG4VUserPostDetConstruction and overloads the Contruct() function so that the VMC (re)assigns the G4 magnetic fields to volumes using the input configuration file defined in this class constructor, together with the code below used in the addVMCFields function in python/geomGeant4.py:

geom = ROOT.TG4GeometryManager.Instance() geom.SetUserPostDetConstruction(fieldMaker) geom.ConstructSDandField()

Definition at line 35 of file ShipFieldMaker.h.

Member Typedef Documentation

SFMap

typedef std::map<TString, TVirtualMagField*> ShipFieldMaker::SFMap

Typedef for a TString-TVirtualMagField* map.

Definition at line 47 of file ShipFieldMaker.h.

stringVect

typedef std::vector<std::string> ShipFieldMaker::stringVect

Typedef of a vector of strings.

Definition at line 50 of file ShipFieldMaker.h.

Constructor & Destructor Documentation

ShipFieldMaker()

ShipFieldMaker::ShipFieldMaker

(

Bool_t

verbose = kFALSE

)

Constructor.

Definition at line 37 of file ShipFieldMaker.cxx.

                                             :
    TG4VUserPostDetConstruction(),
    globalField_(0),
    theFields_(),
    regionInfo_(),
    localInfo_(),
    verbose_(verbose),
    Tesla_(10.0), // To convert T to kGauss for VMC/FairRoot
    theNode_(0),
    gotNode_(kFALSE)
{
}

~ShipFieldMaker()

ShipFieldMaker::~ShipFieldMaker ( )

virtual

Destructor.

Definition at line 50 of file ShipFieldMaker.cxx.

{
 
    // Delete the various magnetic fields
    SFMap::iterator iter;
    for (iter = theFields_.begin(); iter != theFields_.end(); ++iter) {
 
    delete iter->second;
 
    }
 
    theFields_.clear();
 
}

Member Function Documentation

checkLocalFieldMap()

void ShipFieldMaker::checkLocalFieldMap ( TVirtualMagField *&  localField, const TString &  volName, Double_t  scale  )

protected

Check if we have a local field map and store the volume global transformation.

Parameters

[in]	localField	The pointer (reference) to the field map (which may be updated)
[in]	volName	The name of the volume (which is used to find the transformation)
[in]	scale	The B field magnitude scaling factor

Definition at line 951 of file ShipFieldMaker.cxx.

                                    {
 
    // We assume that local field maps are stored using co-ordinates centred
    // on the volume. However, GetField() uses global co-ordinates, so we need
    // to find the global volume transformation (translation and rotation) so
    // that we can find the equivalent local co-ordinates. This also means that each
    // local volume needs to have its own lightweight field map copy, where we
    // reuse the field map data but just change the co-ordinate transformation info
 
    ShipBFieldMap* mapField = dynamic_cast<ShipBFieldMap*>(localField);
    if (mapField) {
 
    TString fieldName(mapField->GetName());
    TString localName(fieldName); localName += volName;
 
    if (verbose_) {
        std::cout<<"Checking local field map "<<fieldName
             <<" co-ord transformation for volume "<<volName<<std::endl;
    }
 
    // Check if we already have the local map to avoid duplication
    ShipBFieldMap* localMap = dynamic_cast<ShipBFieldMap*>(this->getField(localName));
    
    if (!localMap && volName.Length() > 0) {
 
        // Get the volume and its associate global transformation
        TString volName1(volName); volName1 += "_1";
 
        transformInfo theInfo;
        this->getTransformation(volName1, theInfo);
        
        // The original field map may have defined its own translation and rotation.
        // Apply this before the volume global transformation
        Double_t origX0 = mapField->GetXOffset();
        Double_t origY0 = mapField->GetYOffset();
        Double_t origZ0 = mapField->GetZOffset();
        TGeoTranslation origTrans("origTrans", origX0, origY0, origZ0);
 
        Double_t origPhi = mapField->GetPhi();
        Double_t origTheta = mapField->GetTheta();
        Double_t origPsi = mapField->GetPsi();
        TGeoRotation origRot("origRot", origPhi, origTheta, origPsi);
 
        TGeoCombiTrans origComb(origTrans, origRot);
        if (verbose_) {
        std::cout<<"The original field map transformation:"<<std::endl;
        origComb.Print();
        }
 
        TGeoTranslation newTrans("newTrans", theInfo.x0_, theInfo.y0_, theInfo.z0_);
        TGeoRotation newRot("newRot", theInfo.phi_, theInfo.theta_, theInfo.psi_);
        
        TGeoCombiTrans newComb(newTrans, newRot);
 
        if (verbose_) {
        std::cout<<"The volume transformation:"<<std::endl;
        newComb.Print();
        }
 
        // The full transformation
        newComb = newComb*origComb;
 
        if (verbose_) {
        std::cout<<"The combined transformation:"<<std::endl;
        newComb.Print();
        }
 
        // Update transformation info
        const Double_t* newTransArray = newComb.GetTranslation();
        theInfo.x0_ = newTransArray[0];
        theInfo.y0_ = newTransArray[1];
        theInfo.z0_ = newTransArray[2];
        
        const TGeoRotation* fullRot = newComb.GetRotation();
        if (fullRot) {
        fullRot->GetAngles(theInfo.phi_, theInfo.theta_, theInfo.psi_);
        }
 
        // Create a lightweight copy, reusing the map data but just updating
        // the global transformation
        if (verbose_) {
        std::cout<<"Creating field map copy for "<<localName
             <<" based on "<<mapField->GetName()
             <<": x0 = "<<theInfo.x0_<<", y0 = "<<theInfo.y0_
             <<", z0 = "<<theInfo.z0_<<", phi = "<<theInfo.phi_
             <<", theta = "<<theInfo.theta_
             <<", psi = "<<theInfo.psi_<<", scale = "<<scale
             <<" and symmetry = "<<mapField->HasSymmetry()<<std::endl;
        }
 
        localMap = new ShipBFieldMap(localName.Data(), *mapField, 
                     theInfo.x0_, theInfo.y0_, theInfo.z0_,
                     theInfo.phi_, theInfo.theta_, theInfo.psi_, scale);
        // Keep track that we have created this field pointer
        theFields_[localName] = localMap;
 
    }
    
    // Set the localField pointer to use the (new or already existing) localMap pointer
    localField = localMap;
    
    }
 
}

ClassDef()

ShipFieldMaker::ClassDef	(	ShipFieldMaker	,
		1
	)

Generate fieldMap csv file in the given region.

ClassDef for ROOT

Construct()

void ShipFieldMaker::Construct ( )

virtual

Set-up all local and regional fields and assign them to volumes.

Definition at line 65 of file ShipFieldMaker.cxx.

{
 
    // Assign volumes with their regional (local + global) or local B fields
    this->setAllRegionFields();
    this->setAllLocalFields();
 
}

defineBell() [1/2]

void ShipFieldMaker::defineBell ( const stringVect &  inputLine )

protected

Define a Bell field based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 320 of file ShipFieldMaker.cxx.

{
 
    size_t nWords = inputLine.size();
 
    // Expecting a line such as:
    // Bell Name BPeak zMiddle orientationInt tubeRadius
 
    if (nWords == 6 || nWords == 9) {
 
    TString name(inputLine[1].c_str());
 
    Double_t BPeak = std::atof(inputLine[2].c_str());
    Double_t zMiddle = std::atof(inputLine[3].c_str()); // cm
 
    Int_t orient = std::atoi(inputLine[4].c_str());
    Double_t tubeR = std::atof(inputLine[5].c_str()); // cm
 
    Double_t xy(0.0), z(0.0), L(0.0);
 
    if (nWords == 9) {
        // Specify "target" dimensions (cm)
        xy = std::atof(inputLine[6].c_str());
        z  = std::atof(inputLine[7].c_str());
        L  = std::atof(inputLine[8].c_str());
    }
 
    this->defineBell(name, BPeak, zMiddle, orient, tubeR, xy, z, L);
 
    } else {
 
    std::cout<<"Expecting 6 or 9 words for the definition of the Bell field: "
         <<"Bell Name BPeak zMiddle orientationInt tubeRadius [targetXY targetZ0 targetL]"<<std::endl;
 
    }
 
}

defineBell() [2/2]

void ShipFieldMaker::defineBell	(	const TString &	name,
		Double_t	BPeak,
		Double_t	zMiddle,
		Int_t	orient = 1,
		Double_t	tubeR = 500.0,
		Double_t	xy = 0.0,
		Double_t	z = 0.0,
		Double_t	L = 0.0
	)

Define a Bell field.

Parameters

[in]	name	The name of the field
[in]	BPeak	The peak B field magnitude (Tesla)
[in]	zMiddle	The middle z global position of the magnet (cm)
[in]	orient	Orientation flag: 1 => Bx = 0 (default), 0 => By = 0
[in]	tubeR	The largest inner radius of the tube ellipse (cm); default = 500 cm
[in]	xy	Optional target xy radius region (cm)
[in]	z	Optional target start z global position (cm)
[in]	L	Optional target region length (cm)

Definition at line 358 of file ShipFieldMaker.cxx.

{
    
    if (!this->gotField(name)) {
    
    ShipBellField* theField = new ShipBellField(name.Data(), BPeak*Tesla_, zMiddle, orient, tubeR);
 
    // Set additional parameters if we have a non-zero target length
    if (fabs(L) > 0.0) {
        theField->IncludeTarget(xy, z, L);
    }
 
    theFields_[name] = theField;
    
    } else {
 
    if (verbose_) {
        std::cout<<"We already have a Bell field with the name "
             <<name.Data()<<std::endl;
    }
    }
}

defineComposite() [1/3]

void ShipFieldMaker::defineComposite ( const stringVect &  inputLine )

protected

Define a composite field based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 557 of file ShipFieldMaker.cxx.

{
 
    size_t nWords = inputLine.size();
 
    // Expecting a line such as:
    // Composite Name Field1 Field2 ... FieldN
 
    if (nWords > 2) {
 
    TString name(inputLine[1].c_str());
 
    std::vector<TString> fieldNames;
    for (size_t i = 2; i < nWords; i++) {
 
        TString aName(inputLine[i].c_str());
        fieldNames.push_back(aName);
 
    }
 
    this->defineComposite(name, fieldNames);
 
    } else {
 
    std::cout<<"Expecting at least 3 words for the composite definition: "
         <<"Composite Name Field1 Field2 ... FieldN"<<std::endl;
 
    }
 
}

defineComposite() [2/3]

void ShipFieldMaker::defineComposite	(	const TString &	name,
		const TString &	field1Name,
		const TString &	field2Name,
		const TString &	field3Name = "",
		const TString &	field4Name = ""
	)

Define a composite field from up to four fields.

Parameters

[in]	name	The name of the composite field
[in]	field1Name	The name of the first field
[in]	field2Name	The name of the second field
[in]	field3Name	The name of the third field (optional)
[in]	field4Name	The name of the fourth field (optional)

Definition at line 588 of file ShipFieldMaker.cxx.

{
 
    std::vector<TString> fieldNames;
    fieldNames.push_back(field1Name);
    fieldNames.push_back(field2Name);
 
    if (field3Name.Length() > 0) {
    fieldNames.push_back(field3Name);
    }
    if (field4Name.Length() > 0) {
    fieldNames.push_back(field4Name);
    }
 
    this->defineComposite(name, fieldNames);
 
}

defineComposite() [3/3]

void ShipFieldMaker::defineComposite	(	const TString &	name,
		std::vector< TString >	fieldNames
	)

Define a composite field from a vector of field names.

Parameters

[in]	name	The name of the composite field
[in]	fieldNames	Vector of all of the field names to be combined

Definition at line 608 of file ShipFieldMaker.cxx.

{
 
    // Check if the field is already in the map
    if (!this->gotField(name)) {
 
    // Loop over the list of fields and add them to the composite
    std::vector<TVirtualMagField*> vectFields;
    
    std::vector<TString>::iterator iter;
    for (iter = fieldNames.begin(); iter != fieldNames.end(); ++iter) {
 
        TString aName = *iter;
        TVirtualMagField* aField = this->getField(aName);
        if (aField) {
        if (verbose_) {std::cout<<"Adding field "<<aName<<std::endl;}
        vectFields.push_back(aField);
        }
 
        ShipCompField* composite = new ShipCompField(name.Data(), vectFields);
        theFields_[name] = composite;
 
    }
 
    } else {
 
    std::cout<<"We already have a composite field with the name "
         <<name.Data()<<std::endl;
    }
 
}

defineConstant() [1/2]

void ShipFieldMaker::defineConstant ( const stringVect &  inputLine )

protected

Define a constant field based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 247 of file ShipFieldMaker.cxx.

{
 
    size_t nWords = inputLine.size();
 
    // Expecting a line such as:
    // Constant Name xMin xMax yMin yMax zMin zMax Bx By Bz
 
    if (nWords == 11) {
 
    TString name(inputLine[1].c_str());
 
    Double_t xMin = std::atof(inputLine[2].c_str());
    Double_t xMax = std::atof(inputLine[3].c_str());
 
    Double_t yMin = std::atof(inputLine[4].c_str());
    Double_t yMax = std::atof(inputLine[5].c_str());
 
    Double_t zMin = std::atof(inputLine[6].c_str());
    Double_t zMax = std::atof(inputLine[7].c_str());
 
    // Input field in Tesla_, interface needs kGauss units
    Double_t Bx = std::atof(inputLine[8].c_str());
    Double_t By = std::atof(inputLine[9].c_str());
    Double_t Bz = std::atof(inputLine[10].c_str());
    
    const TVector2 xRange(xMin, xMax);
    const TVector2 yRange(yMin, yMax);
    const TVector2 zRange(zMin, zMax);
    const TVector3 BVector(Bx, By, Bz);
 
    this->defineConstant(name, xRange, yRange, zRange, BVector);
 
    } else {
 
    std::cout<<"Expecting 11 words for the definition of the constant field: "
         <<"Constant Name xMin xMax yMin yMax zMin zMax Bx By Bz"<<std::endl;
 
    }
 
}

defineConstant() [2/2]

void ShipFieldMaker::defineConstant	(	const TString &	name,
		const TVector2 &	xRange,
		const TVector2 &	yRange,
		const TVector2 &	zRange,
		const TVector3 &	BVector
	)

Define a constant field.

Parameters

[in]	name	The name of the field
[in]	xRange	The x range as a TVector2(xMin, xMax)
[in]	yRange	The y range as a TVector2(yMin, yMax)
[in]	zRange	The z range as a TVector2(zMin, zMax)
[in]	BVector	The vector of B field components (Bx, By, Bz) in Tesla

Definition at line 289 of file ShipFieldMaker.cxx.

{
 
    // Check if the field is already in the map
    if (!this->gotField(name)) {
 
    Double_t xMin = xRange.X();
    Double_t xMax = xRange.Y();
    Double_t yMin = yRange.X();
    Double_t yMax = yRange.Y();
    Double_t zMin = zRange.X();
    Double_t zMax = zRange.Y();
 
    Double_t Bx = BVector.X()*Tesla_;
    Double_t By = BVector.Y()*Tesla_;
    Double_t Bz = BVector.Z()*Tesla_;
 
    ShipConstField* theField = new ShipConstField(name.Data(), xMin, xMax, yMin, yMax,
                              zMin, zMax, Bx, By, Bz);
    theFields_[name] = theField;
 
    }  else {
 
    if (verbose_) {
        std::cout<<"We already have a constant field with the name "
             <<name.Data()<<std::endl;
    }
    }
}

defineFieldMap() [1/2]

void ShipFieldMaker::defineFieldMap ( const stringVect &  inputLine, Bool_t  useSymmetry = kFALSE  )

protected

Define a field map based on information from the inputLine.

Parameters

[in]	inputLine	The space separated input line
[in]	useSymmetry	Boolean to specify if we have x-y quadrant symmetry

Definition at line 383 of file ShipFieldMaker.cxx.

{
 
    size_t nWords = inputLine.size();
 
    // Expecting the line:
    // FieldMap Name mapFileName [x0 y0 z0] [phi theta psi]
 
    if (nWords == 3 || nWords == 6 || nWords == 9) {
 
    const TString name(inputLine[1].c_str());
    const TString mapFileName(inputLine[2].c_str());
    
    Double_t x0(0.0), y0(0.0), z0(0.0);
    Double_t phi(0.0), theta(0.0), psi(0.0);
 
    if (nWords > 5) {
        x0 = std::atof(inputLine[3].c_str());
        y0 = std::atof(inputLine[4].c_str());
        z0 = std::atof(inputLine[5].c_str());
    }
 
    if (nWords == 9) {
        phi = std::atof(inputLine[6].c_str());
        theta = std::atof(inputLine[7].c_str());
        psi = std::atof(inputLine[8].c_str());
    }
 
    const TVector3 localCentre(x0, y0, z0);
    const TVector3 localAngles(phi, theta, psi);
 
    this->defineFieldMap(name, mapFileName, localCentre, localAngles, useSymmetry);
 
    } else {
 
    std::cout<<"Expecting 3, 6 or 9 words for the definition of the field map: "
         <<"(Sym)FieldMap Name mapFileName [x0 y0 z0] [[phi theta psi]]"<<std::endl;
 
    }
}

defineFieldMap() [2/2]

void ShipFieldMaker::defineFieldMap	(	const TString &	name,
		const TString &	mapFileName,
		const TVector3 &	localCentre = TVector3(0.0, 0.0, 0.0),
		const TVector3 &	localAngles = TVector3(0.0, 0.0, 0.0),
		Bool_t	useSymmetry = kFALSE
	)

Parameters

[in]	name	The name of the field
[in]	mapFileName	The location of the field map ROOT file relative to VMCWORKDIR
[in]	localCentre	The TVector3(x,y,z) offset shift applied to all field map coordinates
[in]	localAngles	The TVector3(phi, theta, psi) Euler rotation applied to all map coords
[in]	useSymmetry	Boolean to specify if the map has quadrant symmetry (default = false)

Definition at line 424 of file ShipFieldMaker.cxx.

{
    // Check if the field is already in the map
    if (!this->gotField(name)) {
 
    // Add the VMCWORKDIR prefix to this map file location
    std::string fullFileName = getenv("VMCWORKDIR");
    fullFileName += "/"; fullFileName += mapFileName.Data();
 
    if (verbose_) {
        if (useSymmetry) {
        std::cout<<"Creating symmetric field map called "<<name.Data()<<" using "<<fullFileName<<std::endl;
        } else {
        std::cout<<"Creating field map called "<<name.Data()<<" using "<<fullFileName<<std::endl;
        }
    }
 
    // Since field maps use floating point precision we convert any double parameter
    // values to floats, i.e. we can't simply pass TVector3 objects
    Float_t x0 = localCentre.X();
    Float_t y0 = localCentre.Y();
    Float_t z0 = localCentre.Z();
 
    Float_t phi = localAngles.X();
    Float_t theta = localAngles.Y();
    Float_t psi = localAngles.Z();
 
    Float_t scale(1.0);
 
    ShipBFieldMap* mapField = new ShipBFieldMap(name.Data(), fullFileName, x0, y0, z0, 
                            phi, theta, psi, scale, useSymmetry);
    theFields_[name] = mapField;
 
    } else {
 
    if (verbose_) {
        std::cout<<"We already have a field map with the name "
             <<name.Data()<<std::endl;      
    }
 
    }
 
}

defineFieldMapCopy() [1/2]

void ShipFieldMaker::defineFieldMapCopy ( const stringVect &  inputLine )

protected

Define a (translated) copy of a field map based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 470 of file ShipFieldMaker.cxx.

{
 
    // Define a (translated) copy of a field map based in the input file line
 
    size_t nWords = inputLine.size();
 
    // Expecting the line:
    // CopyMap Name MapNameToCopy x0 y0 z0 [phi theta psi]
 
    if (nWords == 6 || nWords == 9) {
 
    const TString name(inputLine[1].c_str());
        
    // We want to try to copy and transpose an already existing field map
    const TString mapNameToCopy(inputLine[2].c_str());
 
    Double_t x0 = std::atof(inputLine[3].c_str());
    Double_t y0 = std::atof(inputLine[4].c_str());
    Double_t z0 = std::atof(inputLine[5].c_str());
 
    Double_t phi(0.0), theta(0.0), psi(0.0), scale(1.0);
 
    if (nWords == 9) {
        phi = std::atof(inputLine[6].c_str());
        theta = std::atof(inputLine[7].c_str());
        psi = std::atof(inputLine[8].c_str());
    }
 
    const TVector3 translation(x0, y0, z0);
    const TVector3 eulerAngles(phi, theta, psi);
 
    this->defineFieldMapCopy(name, mapNameToCopy, translation, eulerAngles);
 
    } else {
 
    std::cout<<"Expecting 6 or 9 words for the copy of a field map: "
         <<"CopyMap Name MapNameToCopy x0 y0 z0 [phi theta psi]"<<std::endl;
 
    }
 
}

defineFieldMapCopy() [2/2]

void ShipFieldMaker::defineFieldMapCopy	(	const TString &	name,
		const TString &	mapNameToCopy,
		const TVector3 &	translation,
		const TVector3 &	eulerAngles = TVector3(0.0, 0.0, 0.0)
	)

Define a copy of a field map with a coordinate translation and optional rotation.

Parameters

[in]	name	The name of the field
[in]	mapNameToCopy	The name of the field map that is to be copied
[in]	translation	The TVector3(x,y,z) coordinate translation
[in]	eulerAngles	The TVector3(phi, theta, psi) Euler angle rotation

Definition at line 513 of file ShipFieldMaker.cxx.

{
 
    // Check if the field is already in the map
    if (!this->gotField(name)) {
 
    ShipBFieldMap* fieldToCopy = 
        dynamic_cast<ShipBFieldMap*>(this->getField(mapNameToCopy));
 
    if (fieldToCopy) {
        
        if (verbose_) {
        std::cout<<"Creating map field copy "<<name.Data()
             <<" based on "<<mapNameToCopy.Data()<<std::endl;
        }
 
        // Since field maps use floating point precision we convert any double parameter
        // values to floats, i.e. we can't simply pass TVector3 objects
        Float_t x0 = translation.X();
        Float_t y0 = translation.Y();
        Float_t z0 = translation.Z();
 
        Float_t phi = eulerAngles.X();
        Float_t theta = eulerAngles.Y();
        Float_t psi = eulerAngles.Z();
        
        Float_t scale(1.0);
 
        ShipBFieldMap* copiedMap = new ShipBFieldMap(name.Data(), *fieldToCopy,
                             x0, y0, z0, phi, theta, psi, scale);
        theFields_[name] = copiedMap;           
        
    }
 
    } else {
 
    std::cout<<"We already have a copied field map with the name "
         <<name.Data()<<std::endl;
    }
 
}

defineGlobalField() [1/3]

void ShipFieldMaker::defineGlobalField ( const stringVect &  inputLine )

protected

Define the global field based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 640 of file ShipFieldMaker.cxx.

{
 
    // Define the global field using the input file line
 
    size_t nWords = inputLine.size();
 
    // Expecting the line:
    // Global Field1 ... FieldN
 
    if (nWords > 1) {
 
    TString name("Global");
 
    std::vector<TString> fieldNames;
    for (size_t i = 1; i < nWords; i++) {
 
        TString aName(inputLine[i].c_str());
        fieldNames.push_back(aName);
 
    }
 
    this->defineGlobalField(fieldNames);
 
    } else {
 
    std::cout<<"Expecting at least two words for the global field definition: "
         <<"Global Field1 ... FieldN"<<std::endl;
 
    }
 
}

defineGlobalField() [2/3]

void ShipFieldMaker::defineGlobalField	(	const TString &	field1Name,
		const TString &	field2Name = "",
		const TString &	field3Name = "",
		const TString &	field4Name = ""
	)

Define a composite field from up to four fields.

Parameters

[in]	field1Name	The name of the first field
[in]	field2Name	The name of the second field (optional)
[in]	field3Name	The name of the third field (optional)
[in]	field4Name	The name of the fourth field (optional)

Definition at line 673 of file ShipFieldMaker.cxx.

{
 
    std::vector<TString> fieldNames;
    fieldNames.push_back(field1Name);
 
    if (field2Name.Length() > 0) {
    fieldNames.push_back(field2Name);
    }
    if (field3Name.Length() > 0) {
    fieldNames.push_back(field3Name);
    }
    if (field4Name.Length() > 0) {
    fieldNames.push_back(field4Name);
    }
 
    this->defineGlobalField(fieldNames);
 
}

defineGlobalField() [3/3]

void ShipFieldMaker::defineGlobalField

(

std::vector< TString >

fieldNames

)

Define the Global field using a vector of field names.

Parameters

[in]

fieldNames

Vector of all of the field names to be combined

Definition at line 694 of file ShipFieldMaker.cxx.

{
 
    if (globalField_) {
    if (verbose_) {
        std::cout<<"Deleting already existing Global field"<<std::endl;
    }
    delete globalField_;
    }
 
    if (verbose_) {std::cout<<"Setting the Global field"<<std::endl;}
 
    std::vector<TVirtualMagField*> vectFields;
 
    std::vector<TString>::iterator iter;
    for (iter = fieldNames.begin(); iter != fieldNames.end(); ++iter) {
 
    TString aName = *iter;
    TVirtualMagField* aField = this->getField(aName);
    
    if (aField) {
        if (verbose_) {std::cout<<"Adding field "<<aName<<" to Global"<<std::endl;}
        vectFields.push_back(aField);
    } else {
        std::cout<<"Could not find the field "<<aName<<std::endl;
    }
 
    }
 
    TString name("Global");
    globalField_ = new ShipCompField(name.Data(), vectFields);
 
    // Set this as the global field in the virtual MC
    if (gMC) {
    gMC->SetMagField(globalField_);
    } else {
    std::cout<<"The virtual MC pointer gMC is null! The global field can't be used by Geant4 but will work for track fitting and track extrapolation"<<std::endl;
    }
 
}

defineLocalField() [1/2]

void ShipFieldMaker::defineLocalField ( const stringVect &  inputLine )

protected

Define a local field only based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 865 of file ShipFieldMaker.cxx.

{
 
    // Set only the local field for the region using info from inputLine
    size_t nWords = inputLine.size();
 
    // Expecting the line:
    // Local VolName FieldName [FieldMapScaleFactor]
 
    if (nWords == 3 || nWords == 4) {
    
    TString volName(inputLine[1].c_str());
    TString fieldName(inputLine[2].c_str());
 
    Double_t scale(1.0);
    if (nWords == 4) {scale = std::atof(inputLine[3].c_str());}
 
    this->defineLocalField(volName, fieldName, scale);
 
    } else {
 
    std::cout<<"Expecting 3 or 4 words for the local field definition: "
         <<"Local VolName LocalFieldName [FieldMapScale]"<<std::endl;
 
    }
 
}

defineLocalField() [2/2]

void ShipFieldMaker::defineLocalField	(	const TString &	volName,
		const TString &	fieldName,
		Double_t	scale = 1.0
	)

Define a localised field and volume pairing.

Parameters

[in]	volName	The name of the volume
[in]	fieldName	The name of the local field for the volume
[in]	scale	Optional scale factor for field maps (default = 1.0)

Definition at line 893 of file ShipFieldMaker.cxx.

{
 
    if (verbose_) {
    std::cout<<"ShipFieldMaker::defineLocalField for volume "
         <<volName.Data()<<" and field "<<fieldName.Data()
         <<" with scale = "<<scale<<std::endl;
    }
    
    fieldInfo theInfo(volName, fieldName,scale);
    localInfo_.push_back(theInfo);
 
}

defineRegionField() [1/2]

void ShipFieldMaker::defineRegionField ( const stringVect &  inputLine )

protected

Define a regional (local+global) field based on the info from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 735 of file ShipFieldMaker.cxx.

{
 
    // Set the local + global field for the region using info from inputLine
    size_t nWords = inputLine.size();
 
    // Expecting the line:
    // Region VolName FieldName [FieldMapScaleFactor]
 
    if (nWords == 3 || nWords == 4) {
    
    TString volName(inputLine[1].c_str());
    TString fieldName(inputLine[2].c_str());
 
    Double_t scale(1.0);
    if (nWords == 4) {scale = std::atof(inputLine[3].c_str());}
 
    this->defineRegionField(volName, fieldName, scale);
 
    } else {
 
    std::cout<<"Expecting 3 or 4 words for the region (local + global) field definition: "
         <<"Region VolName LocalFieldName [LocalFieldMapScale]"<<std::endl;
 
    }
 
}

defineRegionField() [2/2]

void ShipFieldMaker::defineRegionField	(	const TString &	volName,
		const TString &	fieldName,
		Double_t	scale = 1.0
	)

Define a regional (local + global) field and volume pairing.

Parameters

[in]	volName	The name of the volume
[in]	fieldName	The name of the field for the volume
[in]	scale	Optional scale factor for field maps (default = 1.0)

Definition at line 763 of file ShipFieldMaker.cxx.

{
 
    if (verbose_) {
    std::cout<<"ShipFieldMaker::defineRegionField for volume "
         <<volName.Data()<<" and field "<<fieldName.Data()
         <<" with scale = "<<scale<<std::endl;
    }
 
    fieldInfo theInfo(volName, fieldName, scale);
    regionInfo_.push_back(theInfo);
 
}

defineUniform() [1/2]

void ShipFieldMaker::defineUniform ( const stringVect &  inputLine )

protected

Define a uniform field based on information from the inputLine.

Parameters

[in]

inputLine

The space separated input line

Definition at line 194 of file ShipFieldMaker.cxx.

{
 
    size_t nWords = inputLine.size();
 
    // Expecting a line such as:
    // Uniform Name Bx By Bz
 
    if (nWords == 5) {
 
    const TString name(inputLine[1].c_str());
 
    Double_t Bx = std::atof(inputLine[2].c_str());
    Double_t By = std::atof(inputLine[3].c_str());
    Double_t Bz = std::atof(inputLine[4].c_str());
    const TVector3 BVector(Bx, By, Bz);
 
    this->defineUniform(name, BVector);
 
    } else {
 
    std::cout<<"Expecting 5 words for the definition of the uniform field: "
         <<"Uniform Name Bx By Bz"<<std::endl;
 
    }
 
}

defineUniform() [2/2]

void ShipFieldMaker::defineUniform	(	const TString &	name,
		const TVector3 &	BVector
	)

Define a uniform field.

Parameters

[in]	name	The name of the field
[in]	BVector	The vector of B field components (Bx, By, Bz) in Tesla

Definition at line 222 of file ShipFieldMaker.cxx.

{
 
    // Check if the field is already in the map
    if (!this->gotField(name)) {
 
    if (verbose_) {std::cout<<"Creating uniform field for "<<name.Data()<<std::endl;}
 
    Double_t Bx = BVector.X()*Tesla_;
    Double_t By = BVector.Y()*Tesla_;
    Double_t Bz = BVector.Z()*Tesla_;
 
    TGeoUniformMagField* uField = new TGeoUniformMagField(Bx, By, Bz);
    theFields_[name] = uField;
 
    }  else {
 
    if (verbose_) {
        std::cout<<"We already have a uniform field with the name "
             <<name.Data()<<std::endl;
    }
    }
}

findNode()

void ShipFieldMaker::findNode ( TGeoVolume *  aVolume, const TString &  volName  )

protected

Update the current geometry node pointer that matches the required volume name.

Parameters

[in]	aVolume	The current volume, whose nodes are looked at
[in]	volName	The required volume name we want to find

Definition at line 1115 of file ShipFieldMaker.cxx.

                                                                         {
 
    // Find the geometry node that matches the required volume name
 
    // Immediately exit the function if we have already found the volume
    if (gotNode_) {return;}
 
    if (aVolume) {
 
    TObjArray* volNodes = aVolume->GetNodes();
 
    if (volNodes) {
 
        // Loop over the geometry nodes
        int nNodes = volNodes->GetEntries();
        for (int i = 0; i < nNodes; i++) {
 
        TGeoNode* node = dynamic_cast<TGeoNode*>(volNodes->At(i));
 
        if (node) {
 
            const TString nodeName(node->GetName());
            if (!nodeName.CompareTo(volName, TString::kExact)) {
 
            // We have a match. The node has the transformation we need
            theNode_ = node;
            gotNode_ = kTRUE;
            break;
 
            } else if (node->GetNodes()) {
 
            // We have sub-volumes. Recursively call this function
            this->findNode(node->GetVolume(), volName);
 
            }
 
        }
 
        }
 
    }
 
    }
 
}

generateFieldMap()

void ShipFieldMaker::generateFieldMap	(	TString	fileName,
		const float	step = 2.5,
		const float	xRange = 179,
		const float	yRange = 317,
		const float	zRange = 1515.5,
		const float	zShift = -4996
	)

Definition at line 1383 of file ShipFieldMaker.cxx.

                                                                                                                                                       {
        std::ofstream myfile;
        myfile.open (fileName);
        int xSteps = ceil(xRange/step) + 1;   //field map has X range from 0 to xMax
        int ySteps = ceil(yRange/step) + 1;   //from 0 up to yMax
        int zSteps = ceil(zRange*2./step) + 1;//from -zMax up to zMax
        Double_t position[3] = {0.0, 0.0, 0.0};
        myfile<<xSteps<<"    "<<ySteps<<"    "<<zSteps<<"    "<<step<<"    "<<xRange<<"    "<<yRange<<"    "<<zRange<<std::endl;
        for (int i =0 ; i < xSteps; i++){
                for (int k=0;  k<ySteps;k++){
                        for (int m=0; m<zSteps; m++){
                                Double_t B[3] = {0.0, 0.0, 0.0};
                                Double_t x = step*i;
                                Double_t y = step * k;
                                Double_t z = m * step - zRange + zShift;
                                position[0] = x;
                                position[1] = y;
                                position[2] = z;
                                Bool_t inside(kFALSE);
                                TGeoNode* theNode = gGeoManager->FindNode(position[0], position[1], position[2]);
                                if (theNode) {
                                        TGeoVolume* theVol = theNode->GetVolume();
                                        if (theVol) {
                                                TVirtualMagField* theField = dynamic_cast<TVirtualMagField*>(theVol->GetField());
                                                if (theField) {  
                                                    theField->Field(position, B);
                                                    inside = kTRUE;  
                                                }
                                        }
                                }
                                if (inside == kFALSE && globalField_) {
                                        globalField_->Field(position, B);
                                }
                                myfile<<x<<"    "<<y<<"    "<<z<<"    "<<B[0]/Tesla_<<"    "<<B[1]/Tesla_<<"    "<<B[2]/Tesla_<<std::endl;
                        }
                }
        }
        myfile.close();
}

getAllFields()

SFMap ShipFieldMaker::getAllFields ( ) const

inline

Get the map storing volume names and their associated local magnetic fields.

Returns

the map of volume names and their corresponding magnetic field pointers

Definition at line 201 of file ShipFieldMaker.h.

{return theFields_;}

getField()

TVirtualMagField * ShipFieldMaker::getField

(

const TString &

name

)

const

Get the field stored with the given name name.

Parameters

[in]

name

The name of the field

Returns

the pointer to the magnetic field object

Definition at line 1201 of file ShipFieldMaker.cxx.

{
  
    TVirtualMagField* theField(0);
 
    // Iterate over the internal map and see if we have a match 
    SFMap::const_iterator iter;
    for (iter = theFields_.begin(); iter != theFields_.end(); ++iter) {
   
    TString key = iter->first;
    TVirtualMagField* BField = iter->second;
 
    // Check that we have the key already stored
    if (!key.CompareTo(name, TString::kExact)) {
        theField = BField;
        break;
    }
 
    }
 
    return theField;
 
}

getGlobalField()

ShipCompField * ShipFieldMaker::getGlobalField ( ) const

inline

Get the global magnetic field.

Returns

the global magnetic field pointer

Definition at line 195 of file ShipFieldMaker.h.

{return globalField_;}

getTransformation()

void ShipFieldMaker::getTransformation ( const TString &  volName, transformInfo &  theInfo  )

protected

Get the transformation matrix for the volume position and orientation.

Parameters

[in]	volName	The name of the volume
[in]	theInfo	The transformation information structure

Definition at line 1057 of file ShipFieldMaker.cxx.

                                                                                     {
 
    // Find the geometry node that matches the volume name. We need to search
    // the geometry hierarchy recursively until we get a match. Note that nodes
    // have "_1" appended to the volume name.
 
    theInfo.x0_ = 0.0; theInfo.y0_ = 0.0; theInfo.z0_ = 0.0;
    theInfo.phi_ = 0.0; theInfo.theta_ = 0.0; theInfo.psi_ = 0.0;
 
    TGeoMatrix* theMatrix(0);
 
    TGeoVolume* topVolume = gGeoManager->GetTopVolume();
    if (!topVolume) {
    std::cout<<"Can't find the top volume in ShipFieldMaker::getTransformation"<<std::endl;
    return;
    }
 
    if (verbose_) {std::cout<<"Finding transformation for "<<volName<<std::endl;}
 
    // Find the geometry node that matches the required name
    theNode_ = 0;
    gotNode_ = kFALSE;
    this->findNode(topVolume, volName);
 
    if (theNode_) {theMatrix = theNode_->GetMatrix();}
 
    // Retrieve the translation and rotation information
    if (theMatrix) {
 
    // Translation displacement components
    const Double_t* theTrans = theMatrix->GetTranslation();
    theInfo.x0_ = theTrans[0];
    theInfo.y0_ = theTrans[1];
    theInfo.z0_ = theTrans[2];
 
    // Euler rotation angles. First check if we have a combined translation
    // and rotation, then check if we just have a pure rotation
    if (theMatrix->IsCombi()) {
 
        TGeoCombiTrans* theCombi = dynamic_cast<TGeoCombiTrans*>(theMatrix);
        if (theCombi) {     
        TGeoRotation* combRotn = theCombi->GetRotation();
        if (combRotn) {
            combRotn->GetAngles(theInfo.phi_, theInfo.theta_, theInfo.psi_);
        }
        }
 
    } else if (theMatrix->IsRotation()) {
 
        TGeoRotation* theRotn = dynamic_cast<TGeoRotation*>(theMatrix);
        if (theRotn) {
        theRotn->GetAngles(theInfo.phi_, theInfo.theta_, theInfo.psi_);
        }
    }
    }
 
}

getVolumeField()

TVirtualMagField * ShipFieldMaker::getVolumeField

(

const TString &

volName

)

const

Get the magnetic field for the given volume.

Parameters

[in]

volName

The name of the TGeo volume

Returns

the pointer of the local magnetic field for the volume

Definition at line 1161 of file ShipFieldMaker.cxx.

{
 
    TVirtualMagField* theField(0);
    TGeoVolume* theVol(0);
    if (gGeoManager) {theVol = gGeoManager->FindVolumeFast(volName.Data());}
 
    if (theVol) {
    // Need to cast the TObject* to a TVirtualMagField*
    theField = dynamic_cast<TVirtualMagField*>(theVol->GetField());
    }
 
    return theField;
 
}

gotField()

Bool_t ShipFieldMaker::gotField

(

const TString &

name

)

const

Check if we have a field stored with the given name name.

Parameters

[in]

name

The name of the field

Returns

a boolean to say if we already have the field stored in the internal map

Definition at line 1177 of file ShipFieldMaker.cxx.

{
 
    Bool_t result(kFALSE);
 
    // Iterate over the internal map and see if we have a match 
    SFMap::const_iterator iter;
    for (iter = theFields_.begin(); iter != theFields_.end(); ++iter) {
   
    TString key = iter->first;
    TVirtualMagField* theField = iter->second;
 
    // Check that we have the key already stored and the pointer is not null
    if (!key.CompareTo(name, TString::kExact) && theField) {
        result = kTRUE;
        break;
    }
 
    }
 
    return result;
 
}

plotField()

void ShipFieldMaker::plotField	(	Int_t	type,
		const TVector3 &	xAxis,
		const TVector3 &	yAxis,
		const std::string &	plotFile
	)		const

Plot the magnetic field in "x-y" plane.

Parameters

[in]	type	The co-ordinate type: 0 = x-y, 1 = z-x and 2 = z-y
[in]	xAxis	Three vector specifying the min, max and bin width of the x axis
[in]	yAxis	Three vector specifying the min, max and bin width of the y axis
[in]	plotFile	The name of the output file containing the plot of the magnetic field

Definition at line 1243 of file ShipFieldMaker.cxx.

{
 
    std::cout<<"ShipFieldMaker plotField "<<plotFile<<": type = "<<type<<std::endl;
 
    Double_t xMin = xAxis(0);
    Double_t xMax = xAxis(1);
    Double_t dx   = xAxis(2);
    Int_t Nx(0);
    if (dx > 0.0) {Nx = static_cast<Int_t>(((xMax - xMin)/dx) + 0.5);}
    
    Double_t yMin = yAxis(0);
    Double_t yMax = yAxis(1);
    Double_t dy   = yAxis(2);
    Int_t Ny(0);
    if (dy > 0.0) {Ny = static_cast<Int_t>(((yMax - yMin)/dy) + 0.5);}
 
    // Create a 2d histogram
    const int nhistograms = 4; //x,y,z,and magnitude
    const int ncoordinates = 3; //x,y,z
    
    TH2D theHist[nhistograms]; 
    std::string titles[nhistograms] = {"Bx (T)","By (T)","Bz (T)","B (T)"};
    for (int icomponent = 0; icomponent< nhistograms; icomponent++){
      theHist[icomponent] = TH2D(Form("theHist[%i]",icomponent), titles[icomponent].data(), Nx, xMin, xMax, Ny, yMin, yMax);
      theHist[icomponent].SetDirectory(0);
      if (type == 0) {
    // x-y
    theHist[icomponent].SetXTitle("x (cm)"); 
    theHist[icomponent].SetYTitle("y (cm)");
     } 
      else if (type == 1) {
    // z-x
    theHist[icomponent].SetXTitle("z (cm)"); 
    theHist[icomponent].SetYTitle("x (cm)");
      } 
      else if (type == 2) {
    // z-y
    theHist[icomponent].SetXTitle("z (cm)"); 
    theHist[icomponent].SetYTitle("y (cm)");
      }
    }
    // Get list of volumes (to check for local fields)
    TObjArray* theVolumes = gGeoManager->GetListOfVolumes();
    Int_t nVol(0);
    if (theVolumes) {nVol = theVolumes->GetSize();}
 
    // Loop over "x" axis
    for (Int_t ix = 0; ix < Nx; ix++) {
 
    Double_t x = dx*ix + xMin;
 
    // Loop over "y" axis
    for (Int_t iy = 0; iy < Ny; iy++) {
 
        Double_t y = dy*iy + yMin;
 
        // Initialise the B field array to zero
        Double_t B[ncoordinates] = {0.0, 0.0, 0.0};
 
        // Initialise the position array to zero
        Double_t position[ncoordinates] = {0.0, 0.0, 0.0};
        if (type == 0) {
        // x-y
            position[0] = x, position[1] = y;
        } else if (type == 1) {
        // z-x
        position[0] = y, position[2] = x;
        } else if (type == 2) {
        // z-y
        position[1] = y; position[2] = x;
        }
 
        // Find out if the point is inside one of the geometry volumes
        Bool_t inside(kFALSE);
 
        // Find the geoemtry node (volume path)
        TGeoNode* theNode = gGeoManager->FindNode(position[0], position[1], position[2]);
        
        if (theNode) {
        
        // Get the volume
        TGeoVolume* theVol = theNode->GetVolume();
 
        if (theVol) {
 
            // Get the magnetic field
            TVirtualMagField* theField = dynamic_cast<TVirtualMagField*>(theVol->GetField());
 
            if (theField) {
 
            // Find the "local" field inside the volume (using global co-ords)
            theField->Field(position, B);
            inside = kTRUE;
 
            } // theField
 
        } // volume
 
        } // node
 
        // If no local volumes found, use global field if it exists
        if (inside == kFALSE && globalField_) {
        globalField_->Field(position, B);
        }
            
        // Divide by the Tesla_ factor, since we want to plot Tesla_ not kGauss (VMC/FairRoot units)
        for (int icomponent = 0; icomponent<ncoordinates; icomponent++){
          theHist[icomponent].Fill(x,y, B[icomponent]/Tesla_);
        }
        Double_t BMag = sqrt(B[0]*B[0] + B[1]*B[1] + B[2]*B[2])/Tesla_;
        theHist[3].Fill(x, y, BMag);
 
    } // "y" axis
 
    } // "x" axis
 
    Bool_t wasBatch = gROOT->IsBatch();
    // Disable pop-up windows
    gROOT->SetBatch(kTRUE);
    TCanvas theCanvas("theCanvas", "", 900, 700);
    gROOT->SetStyle("Plain");
    gStyle->SetOptStat(0);
    // Set colour style
    gStyle->SetPalette(kBird);
    theCanvas.UseCurrentStyle();
 
    theCanvas.Divide(2,2);
    for (int icomponent = 0; icomponent < nhistograms; icomponent++){
        theCanvas.cd(icomponent+1);
        theHist[icomponent].Draw("colz");
    }
    theCanvas.Print(plotFile.c_str());
 
    // Reset the batch boolean
    if (wasBatch == kFALSE) {gROOT->SetBatch(kFALSE);}
 
}

plotXYField()

void ShipFieldMaker::plotXYField	(	const TVector3 &	xAxis,
		const TVector3 &	yAxis,
		const std::string &	plotFile
	)		const

Plot the magnetic field in x-y plane.

Parameters

[in]	xAxis	Three vector specifying the min, max and bin width of the x axis
[in]	yAxis	Three vector specifying the min, max and bin width of the y axis
[in]	plotFile	The name of the output file containing the plot of the magnetic field

Definition at line 1225 of file ShipFieldMaker.cxx.

{    
    this->plotField(0, xAxis, yAxis, plotFile);
}

plotZXField()

void ShipFieldMaker::plotZXField	(	const TVector3 &	zAxis,
		const TVector3 &	xAxis,
		const std::string &	plotFile
	)		const

Plot the magnetic field in z-x plane.

Parameters

[in]	zAxis	Three vector specifying the min, max and bin width of the z axis
[in]	xAxis	Three vector specifying the min, max and bin width of the x axis
[in]	plotFile	The name of the output file containing the plot of the magnetic field

Definition at line 1231 of file ShipFieldMaker.cxx.

{    
    this->plotField(1, zAxis, xAxis, plotFile);
}

plotZYField()

void ShipFieldMaker::plotZYField	(	const TVector3 &	zAxis,
		const TVector3 &	yAxis,
		const std::string &	plotFile
	)		const

Plot the magnetic field in z-y plane.

Parameters

[in]	zAxis	Three vector specifying the min, max and bin width of the z axis
[in]	yAxis	Three vector specifying the min, max and bin width of the y axis
[in]	plotFile	The name of the output file containing the plot of the magnetic field

Definition at line 1237 of file ShipFieldMaker.cxx.

{    
    this->plotField(2, zAxis, yAxis, plotFile);
}

readInputFile()

void ShipFieldMaker::readInputFile

(

const std::string &

inputFile

)

Read an input file to define fields and associated volumes.

Parameters

[in]

inputFile

The input text file

Definition at line 74 of file ShipFieldMaker.cxx.

{
 
    // Check that we have a non-empty string
    if (inputFile.size() < 1) {
    std::cerr<<"Skipping ShipFieldMaker::readInputFile(): file name is empty"<<std::endl;
    return;
    }
 
    std::string fullFileName = getenv("VMCWORKDIR");
    fullFileName += "/"; fullFileName += inputFile.c_str();
 
    std::cout<<"ShipFieldMaker::makeFields inputFile = "<<fullFileName<<std::endl;
 
    std::ifstream getData(fullFileName);
    std::string whiteSpace(" ");
 
    // Loop while the input file is readable
    while (getData.good()) {
    
    if (getData.peek() == '\n') {
        
        // Finish reading line
        char c;
        getData.get(c);
        
        // Stop while loop if we have reached the end of the file
        if (getData.eof()) {break;}
 
    } else if (getData.peek() == '#') {
        
        // Skip comment line
        getData.ignore(1000, '\n');
        getData.putback('\n');
 
        // Stop while loop if we have reached the end of the file
        if (getData.eof()) {break;}
 
    } else {
 
        // Read data line
        std::string line("");
        std::getline(getData, line);
 
        // Stop while loop if we have reached the end of the file
        if (getData.eof()) {break;}
        
        // Split up the line according to white spaces
        std::vector<std::string> lineVect = this->splitString(line, whiteSpace);
 
        size_t nWords = lineVect.size();
 
        // Check to see if we have at least one keyword at the start of the line
        if (nWords > 1) {
 
        TString keyWord(lineVect[0].c_str());
        keyWord.ToLower();
 
        if (!keyWord.CompareTo("uniform")) {
 
            // Define the uniform magnetic field
            this->defineUniform(lineVect);
 
        } else if (!keyWord.CompareTo("constant")) {
 
            // Define a uniform field with an x,y,z boundary
            this->defineConstant(lineVect);
 
        } else if (!keyWord.CompareTo("bell")) {
 
            // Define the Bell-shaped field
            this->defineBell(lineVect);
 
        } else if (!keyWord.CompareTo("fieldmap")) {
 
            // Define the field map
            this->defineFieldMap(lineVect);
 
        } else if (!keyWord.CompareTo("symfieldmap")) {
 
            // Define the symmetric field map
            this->defineFieldMap(lineVect, kTRUE);
  
        } else if (!keyWord.CompareTo("copymap")) {
 
            // Copy (& translate) the field map
            this->defineFieldMapCopy(lineVect);
 
        } else if (!keyWord.CompareTo("composite")) {
 
            // Define the composite field
            this->defineComposite(lineVect);
 
        } else if (!keyWord.CompareTo("global")) {
 
            // Define which fields are global
            this->defineGlobalField(lineVect);
 
        } else if (!keyWord.CompareTo("region")) {
 
            // Define the local and global fields for the given volume
            this->defineRegionField(lineVect);
 
        } else if (!keyWord.CompareTo("local")) {
 
            // Define the field for the given volume as the local one only
            this->defineLocalField(lineVect);
 
        }
 
        }
 
    }
 
    }
 
    getData.close();
 
}

setAllLocalFields()

void ShipFieldMaker::setAllLocalFields ( )

protected

Definition at line 907 of file ShipFieldMaker.cxx.

{
 
    // Loop over all entries in the localInfo_ vector and assign fields to their volumes
    std::vector<fieldInfo>::iterator localIter;
 
    for (localIter = localInfo_.begin(); localIter != localInfo_.end(); ++localIter) {
 
    fieldInfo theInfo = *localIter;
    const TString volName = theInfo.volName_;
    TString fieldName = theInfo.fieldName_;
    Double_t scale = theInfo.scale_;
 
    TGeoVolume* theVol(0);
    if (gGeoManager) {theVol = gGeoManager->FindVolumeFast(volName.Data());}
 
    if (theVol) {       
 
        TVirtualMagField* localField = this->getField(fieldName);
        
        if (localField) {
 
        this->checkLocalFieldMap(localField, volName, scale);
        theVol->SetField(localField);
 
        if (verbose_) {
            std::cout<<"Setting local field "<<localField->GetName()
                 <<" for volume "<<volName<<std::endl;
        }
 
        } else {
 
        std::cout<<"Could not find the field "<<fieldName.Data()<<std::endl;
        }
 
    } else {
 
        std::cout<<"Could not find the volume "<<volName.Data()<<std::endl;
    }
 
    } // localIter loop
 
}

setAllRegionFields()

void ShipFieldMaker::setAllRegionFields ( )

protected

Definition at line 777 of file ShipFieldMaker.cxx.

{
 
    // Loop over all entries in the regionInfo_ vector and assign fields to their volumes
    std::vector<fieldInfo>::iterator regionIter;
 
    for (regionIter = regionInfo_.begin(); regionIter != regionInfo_.end(); ++regionIter) {
 
    fieldInfo theInfo = *regionIter;
    const TString volName = theInfo.volName_;
    TString fieldName = theInfo.fieldName_;
    Double_t scale = theInfo.scale_;
 
    // Find the volume
    TGeoVolume* theVol(0);
    if (gGeoManager) {theVol = gGeoManager->FindVolumeFast(volName.Data());}
 
    if (theVol) {       
 
        // Find the local field
        TVirtualMagField* localField = this->getField(fieldName);
 
        if (localField) {
 
        // Check local field maps know about their associated volume position and orientation.
        // This will update the localField pointer if required
        this->checkLocalFieldMap(localField, volName, scale);
 
        // Reset the fieldName to use the name from the localField pointer, since this
        // could have changed for a local field map, for example
        fieldName = localField->GetName();
 
        // See if we have already combined this local field with the global field
        if (globalField_ && fieldName.Length() > 0) {
 
            TString lgName(fieldName); lgName += "Global";
            TVirtualMagField* lgField = this->getField(lgName);
 
            if (!lgField) {
 
            // Create the combined local + global field and store in the internal map.
            // Other volumes that use the same combined field will use the stored pointer
            if (verbose_) {
                std::cout<<"Creating the combined field "<<lgName.Data()<<", with local field "
                     <<fieldName.Data()<<" with the global field for volume "
                     <<volName.Data()<<std::endl;
            }
 
            ShipCompField* combField = new ShipCompField(lgName.Data(), localField, globalField_);
            theFields_[lgName] = combField;
            theVol->SetField(combField);
 
            } else {
 
            if (verbose_) {
                std::cout<<"Setting the field "<<lgName.Data()
                     <<" for volume "<<volName.Data()<<std::endl;
            }
            theVol->SetField(lgField);
 
            }
 
        } else {
 
            if (verbose_) {
            std::cout<<"There is no global field defined. Just setting the local field"<<std::endl;
            }
            theVol->SetField(localField);
 
        }
 
        } else {
 
        std::cout<<"Could not find the local field "<<fieldName.Data()<<std::endl;
 
        }
 
 
    } else {
 
        std::cout<<"Could not find the volume "<<volName<<std::endl;
 
    }
 
    } // regionIter loop
 
}

splitString()

ShipFieldMaker::stringVect ShipFieldMaker::splitString ( std::string &  theString, std::string &  splitter  ) const

private

Split a string.

Parameters

[in]	theString	The string to be split up
[in]	splitted	The delimiter that will be used to split up the string

Returns

a vector of the delimiter-separated strings

Definition at line 1423 of file ShipFieldMaker.cxx.

                                                          {
 
    // Code from STLplus
    stringVect result;
 
    if (!theString.empty() && !splitter.empty()) {
 
    for (std::string::size_type offset = 0;;) {
 
        std::string::size_type found = theString.find(splitter, offset);
 
        if (found != std::string::npos) {
        std::string tmpString = theString.substr(offset, found-offset);
        if (tmpString.size() > 0) {result.push_back(tmpString);}
        offset = found + splitter.size();
        } else {
        std::string tmpString = theString.substr(offset, theString.size()-offset);
        if (tmpString.size() > 0) {result.push_back(tmpString);}
        break;
        }
    }
    }
    
    return result;
 
}

Member Data Documentation

globalField_

ShipCompField* ShipFieldMaker::globalField_

private

The global magnetic field.

Definition at line 374 of file ShipFieldMaker.h.

gotNode_

Bool_t ShipFieldMaker::gotNode_

private

Boolean to specify if we have found the volume node we need.

Definition at line 395 of file ShipFieldMaker.h.

localInfo_

std::vector<fieldInfo> ShipFieldMaker::localInfo_

private

Vector of fieldInfo for local fields.

Definition at line 383 of file ShipFieldMaker.h.

regionInfo_

std::vector<fieldInfo> ShipFieldMaker::regionInfo_

private

Vector of fieldInfo for regional fields.

Definition at line 380 of file ShipFieldMaker.h.

Tesla_

Double_t ShipFieldMaker::Tesla_

private

Double converting Tesla to kiloGauss (for VMC/FairRoot B field units)

Definition at line 389 of file ShipFieldMaker.h.

theFields_

SFMap ShipFieldMaker::theFields_

private

The map storing all created magnetic fields.

Definition at line 377 of file ShipFieldMaker.h.

theNode_

TGeoNode* ShipFieldMaker::theNode_

private

The current volume node: used for finding volume transformations.

Definition at line 392 of file ShipFieldMaker.h.

verbose_

Bool_t ShipFieldMaker::verbose_

private

Verbose boolean.

Definition at line 386 of file ShipFieldMaker.h.

---

The documentation for this class was generated from the following files:

• field/ShipFieldMaker.h
• field/ShipFieldMaker.cxx