composite_calorimeter/src/CCalMagneticField.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCalMagneticField.cc
0028 // Description: User Field class implementation.
0029 ///////////////////////////////////////////////////////////////////////////////
0030 #include <fstream>
0031
0032 #include "CCalMagneticField.hh"
0033 #include "CCalutils.hh"
0034 #include "G4SystemOfUnits.hh"
0035 #include "G4FieldManager.hh"
0036
0037 //#define ddebug
0038 //#define debug
0039
0040 //Constructor and destructor:
0041
0042 CCalMagneticField::CCalMagneticField(const G4String &filename) :
0043   fval(0), pos(0), slope(0), intercept(0)
0044 {
0045 #ifdef debug
0046   fVerbosity = 1;
0047 #else
0048   fVerbosity = 0;
0049 #endif
0050
0051   //Let's open the file
0052   G4cout << " ==> Opening file " << filename << " to read magnetic field..."
0053          << G4endl;
0054   G4String pathName = std::getenv("CCAL_GLOBALPATH");
0055   std::ifstream is;
0056   G4bool ok = openGeomFile(is, pathName, filename);
0057
0058   if (ok) {
0059     findDO(is, G4String("FLDM"));
0060     is >> fval >> npts >> xoff;
0061
0062     if (fVerbosity)
0063       G4cout << "Field value " << fval << " # points " << npts <<
0064         " offset in x "
0065              << xoff*mm << G4endl;
0066
0067     if (npts > 0) {
0068       pos       = new G4double[npts];
0069       slope     = new G4double[npts];
0070       intercept = new G4double[npts];
0071
0072       for (G4int i = 0; i < npts; i++) {
0073         is >> pos[i] >> slope[i] >> intercept[i];
0074         if (fVerbosity)
0075           G4cout << tab << "Position " << i << " " << pos[i] << " Slope "
0076                  << slope[i] << " Intercept " << intercept[i] << G4endl;
0077       }
0078     }
0079
0080     ///////////////////////////////////////////////////////////////
0081     // Close the file
0082     G4cout << " ==> Closing file " << filename << G4endl;
0083     is.close();
0084   }
0085 }
0086
0087
0088 CCalMagneticField::~CCalMagneticField() {
0089   if (pos)
0090     delete[] pos;
0091   if (slope)
0092     delete[] slope;
0093   if (intercept)
0094     delete[] intercept;
0095 }
0096
0097
0098 // Member functions
0099
0100 void CCalMagneticField::MagneticField(const G4double x[3], G4double B[3]) const
0101 {
0102   G4int i=0;
0103   for (i=0; i<2; i++) {
0104     B[i]   = 0*kilogauss;
0105   }
0106
0107   G4double m1=0;
0108   G4double c1=0;
0109   G4double xnew = x[0]/mm + xoff;
0110   if (npts > 0) {
0111     for (i=0; i<npts; i++) {
0112       if (xnew > pos[i]*mm) {
0113         m1 = slope[i];
0114         c1 = intercept[i];
0115       }
0116     }
0117   }
0118   G4double scor = c1 + m*xnew;
0119   if (scor < 0.) scor = 0.;
0120   if (scor > 1.) scor = 1.0;
0121
0122   B[2] = scor*fval*kilogauss;
0123   if (fVerbosity)
0124     {
0125
0126       G4cout << "Field at x: " << x[0]/mm << "mm (" << xnew << ") = " <<
0127         B[2]/tesla
0128              << "T (m = " << m1 << ", c = " <<
0129         c1 << ", scale = " << scor << ")"
0130              << G4endl;
0131     }
0132 }
0133
0134
0135 CLHEP::Hep3Vector CCalMagneticField::
0136 MagneticField(const CLHEP::Hep3Vector point) const {
0137
0138   G4double x[3],B[3];
0139   CLHEP::Hep3Vector v;
0140
0141   x[0] = point.x();
0142   x[1] = point.y();
0143   x[2] = point.z();
0144   CCalMagneticField::MagneticField(x, B);
0145   v.setX(B[0]);
0146   v.setY(B[1]);
0147   v.setZ(B[2]);
0148   return v;
0149 }
0150
0151
0152 void CCalMagneticField::GetFieldValue(const G4double x[3], G4double* B) const {
0153   CCalMagneticField::MagneticField(x, B);
0154 }
0155