Changeset 1340 for trunk/source/geometry


Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (15 years ago)
Author:
garnier
Message:

update ti head

Location:
trunk/source/geometry
Files:
39 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/geometry/magneticfield/History

    r1315 r1340  
    1 $Id: History,v 1.146 2009/11/12 15:05:49 japost Exp $
     1$Id: History,v 1.152 2010/09/10 16:02:44 japost Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19Sep 10th,  2010 J.Apostolakis - field-V09-03-03
     20-----------------------------
     21- Revised constructor of G4MagErrorStepper to add number of State variables
     22- Corrected MagErrorStepper to copy State Variable to output
     23- Enable G4ConstRK4 to copy remaining State Variables (must integrate 6)
     24
     25Sep 10th,  2010 J.Apostolakis - field-V09-03-02
     26-----------------------------
     27- Fixed passing of time in G4NystromRK4
     28
     29Jul 21st,  2010 T.Nikitina - field-V09-03-01
     30--------------------------
     31- Fixed cases of memory corruption in G4RKG3_Stepper.
     32- Fixed case of unused array data member in G4ExactHelixStepper.
     33- Removed useless code never executed in G4ConstRK4::Stepper().
     34- Fixed initialization in G4NystromRK4 constructor.
     35
     36Jul 14th,  2010 G.Cosmo - field-V09-03-00
     37-----------------------
     38- Added dummy initialization of members in constructors for G4CashKarpRKF45,
     39  G4ConstRK4, G4EqEMFieldWithEDM, G4EqEMFieldWithSpin, G4ExactHelixStepper,
     40  G4FieldTrack, G4MagHelicalStepper, G4MagInt_Driver, G4Mag_EqRhs,
     41  G4Mag_SpinEqRhs, G4Mag_UsualEqRhs, G4NystromRK4, G4RKG3_Stepper,
     42  G4UniformElectricField, G4UniformElectricField.
     43
     44Nov 12th,  2009 J.Apostolakis - field-V09-02-09
     45-----------------------------
     46- G4MagIntegratorDriver:  activate check on integration error for spin.
    1947
    2048Nov 12th,  2009 J.Apostolakis - field-V09-02-08
    21 -----------------------------
     49-----------------------------  (fix only)
    2250- G4Nystrom:  Corrected interface method getField: array now has explicit dimension[4]
    2351                (Problem found by gcc 4.3 - it checked indices used in inline method! )

  • trunk/source/geometry/magneticfield/include/G4CashKarpRKF45.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4CashKarpRKF45.hh,v 1.11 2008/01/11 15:23:54 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4CashKarpRKF45.hh,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    8888  private:
    8989
    90    // G4int fNumberOfVariables ;  // Already kept in G4MagIntegratorStepper
    91    G4double *ak2, *ak3, *ak4, *ak5, *ak6, *ak7, *yTemp, *yIn;  // scratch space
     90    G4double *ak2, *ak3, *ak4, *ak5, *ak6, *ak7, *yTemp, *yIn;
     91      // scratch space
    9292
    93   // for DistChord calculations
     93    G4double fLastStepLength;
     94    G4double *fLastInitialVector, *fLastFinalVector,
     95             *fLastDyDx, *fMidVector, *fMidError;
     96      // for DistChord calculations
    9497
    95   G4double fLastStepLength;
    96   G4double *fLastInitialVector, *fLastFinalVector,
    97            *fLastDyDx, *fMidVector, *fMidError;
    98   G4CashKarpRKF45* fAuxStepper;
     98    G4CashKarpRKF45* fAuxStepper;
     99
    99100};
    100101

  • trunk/source/geometry/magneticfield/include/G4ConstRK4.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ConstRK4.hh,v 1.2 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ConstRK4.hh,v 1.3 2010/09/10 15:50:17 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5454   public:  // with description
    5555
    56      G4ConstRK4(G4Mag_EqRhs *EquationMotion, G4int numberOfVariables = 8);
     56    G4ConstRK4(G4Mag_EqRhs *EquationMotion, G4int numberOfStateVariables=8);
    5757    ~G4ConstRK4();
    5858

  • trunk/source/geometry/magneticfield/include/G4EqEMFieldWithEDM.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithEDM.hh,v 1.1 2009/11/04 23:51:42 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithEDM.hh,v 1.2 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    7878    G4double fElectroMagCof ;
    7979    G4double fMassCof;
    80 
    8180    G4double omegac;
    8281    G4double anomaly;
    8382    G4double eta;
    84     G4double ParticleCharge;
     83    G4double pcharge;
    8584    G4double E;
    8685    G4double gamma;

  • trunk/source/geometry/magneticfield/include/G4EqEMFieldWithSpin.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithSpin.hh,v 1.3 2008/11/14 13:37:09 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithSpin.hh,v 1.4 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    7373    G4double fElectroMagCof ;
    7474    G4double fMassCof;
    75 
    7675    G4double omegac;
    7776    G4double anomaly;
    78     G4double ParticleCharge;
     77    G4double pcharge;
    7978    G4double E;
    8079    G4double gamma;

  • trunk/source/geometry/magneticfield/include/G4ExactHelixStepper.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4ExactHelixStepper.hh,v 1.5 2007/05/18 12:50:31 tnikitin Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4ExactHelixStepper.hh,v 1.7 2010/07/21 13:45:37 tnikitin Exp $
     27// GEANT4 tag $Name: field-V09-03-03 $
    2828//
    2929//
    … …  
    8282   
    8383  private:
    84     G4ThreeVector    fBfieldValue;   //  Initial value of field at last step
    85     G4ThreeVector    yInitialEHS,  yFinalEHS; 
    86     G4ThreeVector    pInitial;
    87     G4double         fLastStepSize;  // Length of last step
    88     G4double         fYInSav[7];     // Starting state of  x, p, ...
    89      // Values saved for calculating mid-point for chord
    9084
    91      G4Mag_EqRhs*  fPtrMagEqOfMot;
     85    G4ThreeVector    fBfieldValue;
     86      //  Initial value of field at last step
     87    G4Mag_EqRhs*  fPtrMagEqOfMot;
    9288};
    9389

  • trunk/source/geometry/magneticfield/include/G4MagErrorStepper.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagErrorStepper.hh,v 1.11 2006/06/29 18:22:50 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagErrorStepper.hh,v 1.12 2010/09/10 15:52:53 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5252  public:  // with description
    5353
    54     // G4MagErrorStepper(G4Mag_EqRhs *EqRhs, G4int numberOfVariables);
    55     G4MagErrorStepper(G4EquationOfMotion *EqRhs, G4int numberOfVariables);
     54    G4MagErrorStepper(G4EquationOfMotion *EqRhs, G4int numberOfVariables, G4int numStateVariables=12);
    5655    virtual ~G4MagErrorStepper();
    5756 

  • trunk/source/geometry/magneticfield/include/G4MagErrorStepper.icc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagErrorStepper.icc,v 1.13 2006/06/29 18:22:52 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagErrorStepper.icc,v 1.14 2010/09/10 15:52:53 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    3232inline
    3333G4MagErrorStepper::G4MagErrorStepper(G4EquationOfMotion *EquationRhs,
    34                                      G4int numberOfVariables)
    35      : G4MagIntegratorStepper(EquationRhs,numberOfVariables)
     34                                     G4int numberOfVariables,
     35                                     G4int numStateVariables)
     36     : G4MagIntegratorStepper(EquationRhs,numberOfVariables,numStateVariables)
    3637  {
    3738      G4int nvar = std::max(this->GetNumberOfVariables(), 8);

  • trunk/source/geometry/magneticfield/include/G4MagIntegratorDriver.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagIntegratorDriver.hh,v 1.20 2007/05/10 10:10:05 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagIntegratorDriver.hh,v 1.21 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    241241     G4int  fNoTotalSteps, fNoBadSteps, fNoSmallSteps, fNoInitialSmallSteps;
    242242     G4double fDyerr_max, fDyerr_mx2;
    243      G4double fDyerrPos_smTot, fDyerrVel_smTot, fDyerrPos_lgTot, fDyerrVel_lgTot;
    244      G4double fSumH_sm,   fSumH_lg;
     243     G4double fDyerrPos_smTot, fDyerrPos_lgTot, fDyerrVel_lgTot;
     244     G4double fSumH_sm, fSumH_lg;
    245245        // Step Statistics
    246246

  • trunk/source/geometry/magneticfield/include/G4Mag_SpinEqRhs.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_SpinEqRhs.hh,v 1.11 2006/06/29 18:23:09 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_SpinEqRhs.hh,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    7575     G4double omegac;
    7676     G4double anomaly;
    77      G4double ParticleCharge;
    78 
     77     G4double pcharge;
    7978     G4double E;
    8079     G4double gamma;

  • trunk/source/geometry/magneticfield/include/G4NystromRK4.hh

    r1337 r1340  
    2424// ********************************************************************
    2525//
    26 // $Id: G4NystromRK4.hh,v 1.3 2009/11/12 15:01:36 japost Exp $
    27 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     26// $Id: G4NystromRK4.hh,v 1.4 2010/07/14 10:00:36 gcosmo Exp $
     27// GEANT4 tag $Name: field-V09-03-03 $
    2828//
    2929// class G4NystromRK4
    … …  
    123123  G4double dz = P[2]-m_fldPosition[2];
    124124
    125   if((dx*dx+dy*dy+dz*dz) > m_magdistance2) {
    126 
     125  if((dx*dx+dy*dy+dz*dz) > m_magdistance2)
     126  {
    127127    m_fldPosition[0] = P[0];
    128128    m_fldPosition[1] = P[1];

  • trunk/source/geometry/magneticfield/include/G4RKG3_Stepper.hh

    r1337 r1340  
    2626//
    2727//
    28 // $Id: G4RKG3_Stepper.hh,v 1.13 2007/05/18 12:44:02 tnikitin Exp $
    29 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     28// $Id: G4RKG3_Stepper.hh,v 1.14 2010/07/23 14:13:09 tnikitin Exp $
     29// GEANT4 tag $Name: field-V09-03-03 $
    3030//
    3131//
    … …  
    6969    G4double  DistChord() const ;
    7070 
    71     void StepNoErr( const G4double tIn[7],
    72                     const G4double dydx[7],
     71    void StepNoErr( const G4double tIn[8],
     72                    const G4double dydx[6],
    7373                          G4double Step,
    74                           G4double tOut[7],
     74                          G4double tOut[8],
    7575                          G4double B[3] );
    7676      // Integrator RK Stepper from G3 with only two field evaluation per
    … …  
    7979      // B[3] is magnetic field which is passed from substep to substep.
    8080
    81     void StepWithEst( const G4double  tIn[7],
    82                       const G4double dydx[7],
     81    void StepWithEst( const G4double  tIn[8],
     82                      const G4double dydx[6],
    8383                            G4double Step,
    84                             G4double tOut[7],
     84                            G4double tOut[8],
    8585                            G4double& alpha2,
    8686                            G4double& beta2,

  • trunk/source/geometry/magneticfield/src/G4CashKarpRKF45.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4CashKarpRKF45.cc,v 1.15 2008/01/11 18:11:44 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4CashKarpRKF45.cc,v 1.16 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// The Cash-Karp Runge-Kutta-Fehlberg 4/5 method is an embedded fourth
    … …  
    4949                                 G4int noIntegrationVariables,
    5050                                 G4bool primary)
    51   : G4MagIntegratorStepper(EqRhs, noIntegrationVariables)
    52 {
    53   // unsigned int noVariables= std::max(numberOfVariables,8); // For Time .. 7+1
     51  : G4MagIntegratorStepper(EqRhs, noIntegrationVariables),
     52    fLastStepLength(0.), fAuxStepper(0)
     53{
    5454  const G4int numberOfVariables = noIntegrationVariables;
    5555
    … …  
    6969  fMidVector = new G4double[numberOfVariables];
    7070  fMidError =  new G4double[numberOfVariables];
    71   fAuxStepper = 0;   
    72   if( primary )
    73       fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
    74 
     71  if( primary )
     72  {
     73    fAuxStepper = new G4CashKarpRKF45(EqRhs, numberOfVariables, !primary);
     74  }
    7575}
    7676

  • trunk/source/geometry/magneticfield/src/G4ConstRK4.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ConstRK4.cc,v 1.2 2008/10/29 14:17:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ConstRK4.cc,v 1.5 2010/09/10 15:51:10 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    3838//////////////////////////////////////////////////////////////////
    3939//
    40 // Constructor sets the number of variables (default = 8)
    41 
    42 G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numberOfVariables)
    43   : G4MagErrorStepper(EqRhs, numberOfVariables)
    44 {
    45   if(numberOfVariables !=8 )
     40// Constructor sets the number of *State* variables (default = 8)
     41//   The number of variables integrated is always 6
     42
     43G4ConstRK4::G4ConstRK4(G4Mag_EqRhs* EqRhs, G4int numStateVariables)
     44  : G4MagErrorStepper(EqRhs, 6, numStateVariables)
     45{
     46  // const G4int numberOfVariables= 6;
     47  if( numStateVariables < 8 )
    4648  {
     49    G4cerr << "ERROR in G4ConstRK4::G4ConstRK4 "
     50           << "   The number of State variables at least 8 " << G4endl;
     51    G4cerr << "   Instead it is  numStateVariables= " << numStateVariables << G4endl;
    4752    G4Exception("G4ConstRK4::G4ConstRK4()", "InvalidSetup", FatalException,
    48                 "Valid only for number of variables=8. Use another Stepper!");
     53        "Valid only for number of state variables of 8 or more. Use another Stepper!");
    4954  }
    50   else
    51   {
    52     fEq=EqRhs;
    53     yMiddle= new G4double[8];
    54     dydxMid= new G4double[8];
    55     yInitial= new G4double[8];
    56     yOneStep= new G4double[8];
    57 
    58     dydxm = new G4double[8];
    59     dydxt = new G4double[8];
    60     yt    = new G4double[8];
    61     Field[0]=0.;Field[1]=0.;Field[2]=0.;
    62   }
     55
     56  fEq = EqRhs;
     57  yMiddle  = new G4double[8];
     58  dydxMid  = new G4double[8];
     59  yInitial = new G4double[8];
     60  yOneStep = new G4double[8];
     61
     62  dydxm = new G4double[8];
     63  dydxt = new G4double[8];
     64  yt    = new G4double[8];
     65  Field[0]=0.; Field[1]=0.; Field[2]=0.;
    6366}
    6467
    … …  
    150153                           G4double yError [] )
    151154{
    152    const G4int nvar = 8 ;
    153    const G4int maxvar= 8;
     155   const G4int nvar = 6;  // number of variables integrated
     156   const G4int maxvar= GetNumberOfStateVariables();
     157
     158   // Correction for Richardson extrapolation
     159   G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
    154160
    155161   G4int i;
    156 
    157    // Correction for Richardson extrapolation
    158    //
    159    G4double  correction = 1. / ( (1 << IntegratorOrder()) -1 );
    160162   
    161163   // Saving yInput because yInput and yOutput can be aliases for same array
    162 
    163    for (i=0;i<nvar;i++) { yInitial[i]=yInput[i]; }
    164 
    165    yInitial[7]= yInput[7];  // Copy the time in case...even if not really needed
    166    yMiddle[7] = yInput[7];  // Copy the time from initial value
    167    yOneStep[7] = yInput[7]; // As it contributes to final value of yOutput ?
    168    yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
    169    for (i=nvar;i<maxvar;i++) { yOutput[i]=yInput[i]; }
     164   for (i=0;    i<maxvar; i++) { yInitial[i]= yInput[i]; }
     165 
     166   // Must copy the part of the state *not* integrated to the output
     167   for (i=nvar; i<maxvar; i++) { yOutput[i]=  yInput[i]; }
     168
     169   // yInitial[7]= yInput[7];  //  The time is typically needed
     170   yMiddle[7]  = yInput[7];   // Copy the time from initial value
     171   yOneStep[7] = yInput[7];   // As it contributes to final value of yOutput ?
     172   // yOutput[7] = yInput[7];  // -> dumb stepper does it too for RK4
    170173   yError[7] = 0.0;         
    171174

  • trunk/source/geometry/magneticfield/src/G4EqEMFieldWithEDM.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithEDM.cc,v 1.2 2009/11/06 22:31:54 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithEDM.cc,v 1.3 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    4343
    4444G4EqEMFieldWithEDM::G4EqEMFieldWithEDM(G4ElectroMagneticField *emField )
    45       : G4EquationOfMotion( emField )
     45      : G4EquationOfMotion( emField ), fElectroMagCof(0.), fMassCof(0.),
     46        omegac(0.), anomaly(0.0011659208), eta(0.), pcharge(0.), E(0.),
     47        gamma(0.), beta(0.)
    4648{
    47   anomaly = 0.0011659208;
    48   eta = 0.;
    4949}
    5050
    … …  
    6363   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6464
    65    ParticleCharge = particleCharge;
     65   pcharge = particleCharge;
    6666
    6767   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    145145
    146146   G4ThreeVector dSpin
    147      = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
     147     = pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
    148148                               // from Jackson
    149149                               // -uce*Spin.cross(u.cross(EField)) )
    150150                               // but this form has one less operation
    151                                - uce*(u*(Spin*EField) - EField*(Spin*u))
    152                                +eta/2.*(Spin.cross(EField) - ude*(Spin.cross(u))
    153                                         // +Spin.cross(u.cross(Bfield))
    154                                         + (u*(Spin*BField) - BField*(Spin*u))
    155                                        )
    156                              );
    157                              
     151                       - uce*(u*(Spin*EField) - EField*(Spin*u))
     152                       + eta/2.*(Spin.cross(EField) - ude*(Spin.cross(u))
     153                               // +Spin.cross(u.cross(Bfield))
     154                       + (u*(Spin*BField) - BField*(Spin*u)) ) );
     155     
    158156   dydx[ 9] = dSpin.x();
    159157   dydx[10] = dSpin.y();

  • trunk/source/geometry/magneticfield/src/G4EqEMFieldWithSpin.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4EqEMFieldWithSpin.cc,v 1.8 2009/11/06 22:31:35 gum Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4EqEMFieldWithSpin.cc,v 1.9 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    3333//  30.08.2007 Chris Gong, Peter Gumplinger
    3434//  14.02.2009 Kevin Lynch
    35 //  06.11.2009 Hiromi Iinuma see:
    36 //  http://hypernews.slac.stanford.edu/HyperNews/geant4/get/emfields/161.html
     35//  06.11.2009 Hiromi Iinuma
    3736//
    3837// -------------------------------------------------------------------
    … …  
    4443
    4544G4EqEMFieldWithSpin::G4EqEMFieldWithSpin(G4ElectroMagneticField *emField )
    46   : G4EquationOfMotion( emField )
     45  : G4EquationOfMotion( emField ), fElectroMagCof(0.), fMassCof(0.),
     46    omegac(0.), anomaly(0.0011659208), pcharge(0.), E(0.), gamma(0.), beta(0.)
    4747{
    48   anomaly = 0.0011659208;
    4948}
    5049
    … …  
    6362   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6463
    65    ParticleCharge = particleCharge;
     64   pcharge = particleCharge;
    6665
    6766   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    135134
    136135   G4ThreeVector dSpin
    137      = ParticleCharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
     136     = pcharge*omegac*( ucb*(Spin.cross(BField))-udb*(Spin.cross(u))
    138137                               // from Jackson
    139138                               // -uce*Spin.cross(u.cross(EField)) );
    140139                               // but this form has one less operation
    141                                - uce*(u*(Spin*EField) - EField*(Spin*u)) );
     140                      - uce*(u*(Spin*EField) - EField*(Spin*u)) );
    142141
    143142   dydx[ 9] = dSpin.x();

  • trunk/source/geometry/magneticfield/src/G4ExactHelixStepper.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ExactHelixStepper.cc,v 1.9 2008/10/29 14:34:35 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ExactHelixStepper.cc,v 1.11 2010/07/21 13:46:01 tnikitin Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//  Helix a-la-Explicity Euler: x_1 = x_0 + helix(h)
    … …  
    4343#include "G4LineSection.hh"
    4444
    45 
    4645G4ExactHelixStepper::G4ExactHelixStepper(G4Mag_EqRhs *EqRhs)
    4746  : G4MagHelicalStepper(EqRhs),
    48     fBfieldValue(DBL_MAX, DBL_MAX, DBL_MAX), yInitialEHS(DBL_MAX), yFinalEHS(-DBL_MAX)
    49    
     47    fBfieldValue(DBL_MAX, DBL_MAX, DBL_MAX),
     48    fPtrMagEqOfMot(EqRhs)
    5049{
    51    const G4int nvar = 6 ;
    52    G4int i;
    53    for(i=0;i<nvar;i++)  {
    54      fYInSav[i]= DBL_MAX;
    55    }
    56     fPtrMagEqOfMot=EqRhs;
     50  ;
    5751}
    5852
    … …  
    6155void
    6256G4ExactHelixStepper::Stepper( const G4double yInput[],
    63                               const G4double*,
    64                                     G4double hstep,
    65                                     G4double yOut[],
    66                                     G4double yErr[]      )
     57                              const G4double*,
     58                                    G4double hstep,
     59                                    G4double yOut[],
     60                                    G4double yErr[]      )
    6761{ 
    68    const G4int nvar = 6 ;
     62   const G4int nvar = 6;
    6963
    7064   G4int i;
    71    // G4double      yTemp[7], yIn[7] ;
    7265   G4ThreeVector Bfld_value;
    7366
    74    for(i=0;i<nvar;i++)  {
    75       // yIn[i]=     yInput[i];
    76       fYInSav[i]= yInput[i];
    77    }
    78 
    79    MagFieldEvaluate(yInput, Bfld_value) ;       
    80      
    81    // DumbStepper(yIn, Bfld_value, hstep, yTemp);
     67   MagFieldEvaluate(yInput, Bfld_value);       
    8268   AdvanceHelix(yInput, Bfld_value, hstep, yOut);
    8369
    84    // We are assuming a constant field: helix is exact.
    85    for(i=0;i<nvar;i++) {
    86      yErr[i] = 0.0 ;
    87    }
     70  // We are assuming a constant field: helix is exact
     71  //
     72  for(i=0;i<nvar;i++)
     73  {
     74    yErr[i] = 0.0 ;
     75  }
    8876
    89     yInitialEHS = G4ThreeVector( yInput[0],   yInput[1],   yInput[2]);
    90     yFinalEHS   = G4ThreeVector( yOut[0],  yOut[1],  yOut[2]);
    91     fBfieldValue=Bfld_value;
    92    
     77  fBfieldValue=Bfld_value;
    9378}
    9479
    9580void
    9681G4ExactHelixStepper::DumbStepper( const G4double  yIn[],
    97                                    G4ThreeVector   Bfld,
    98                                    G4double  h,
    99                                    G4double  yOut[])
     82                                        G4ThreeVector   Bfld,
     83                                        G4double  h,
     84                                        G4double  yOut[])
    10085{
    10186  // Assuming a constant field: solution is a helix
     87
    10288  AdvanceHelix(yIn, Bfld, h, yOut);
    10389
    104   G4Exception("G4ExactHelixStepper::DumbStepper should not be called.",
    105               "EHS:NoDumbStepper", FatalException, "Stepper must do all the work." );
     90  G4Exception("G4ExactHelixStepper::DumbStepper",
     91              "EHS:NoDumbStepper", FatalException,
     92              "Should not be called. Stepper must do all the work." );
    10693} 
    10794
    … …  
    10996// ---------------------------------------------------------------------------
    11097
    111 G4double G4ExactHelixStepper::DistChord()   const
     98G4double
     99G4ExactHelixStepper::DistChord() const
    112100{
    113101  // Implementation : must check whether h/R >  pi  !!
    114102  //   If( h/R <  pi)   DistChord=h/2*std::tan(Ang_curve/4)
    115103  //   Else             DistChord=R_helix
    116   //
     104
    117105  G4double distChord;
    118106  G4double Ang_curve=GetAngCurve();
    119107
    120      
    121          if(Ang_curve<=pi){
    122            distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
    123          }
    124          else
    125          if(Ang_curve<twopi){
    126            distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
    127          }
    128          else{
    129           distChord=2.*GetRadHelix(); 
    130          }
    131 
    132    
     108  if (Ang_curve<=pi)
     109  {
     110    distChord=GetRadHelix()*(1-std::cos(0.5*Ang_curve));
     111  }
     112  else if(Ang_curve<twopi)
     113  {
     114    distChord=GetRadHelix()*(1+std::cos(0.5*(twopi-Ang_curve)));
     115  }
     116  else
     117  {
     118    distChord=2.*GetRadHelix(); 
     119  }
    133120
    134121  return distChord;
    135  
    136122}   
    137123

  • trunk/source/geometry/magneticfield/src/G4FieldTrack.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4FieldTrack.cc,v 1.14 2007/10/03 15:34:42 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4FieldTrack.cc,v 1.15 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3636     const G4double *SixV = SixVec.SixVector;
    3737     os << " ( ";
    38      os << " X= " << SixV[0] << " " << SixV[1] << " " << SixV[2] << " ";  // Position
    39      os << " V= " << SixV[3] << " " << SixV[4] << " " << SixV[5] << " ";  // Momentum
    40      os << " v2= " << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag();     // mom magnitude
     38     os << " X= " << SixV[0] << " " << SixV[1] << " "
     39                  << SixV[2] << " ";  // Position
     40     os << " V= " << SixV[3] << " " << SixV[4] << " "
     41                  << SixV[5] << " ";  // Momentum
     42     os << " v2= "
     43        << G4ThreeVector(SixV[3], SixV[4], SixV[5]).mag(); // mom magnitude
    4144     os << " mdm= " << SixVec.fMomentumDir.mag();
    4245     os << " l= " << SixVec.GetCurveLength();
    … …  
    5760   fRestMass_c2(restMass_c2),
    5861   fLabTimeOfFlight(LaboratoryTimeOfFlight),
    59    // fProperTimeOfFlight(0.0),
     62   fProperTimeOfFlight(0.),
    6063   // fMomentumDir(pMomentumDirection),
    6164   fChargeState(  charge, magnetic_dipole_moment )
    … …  
    6366  G4double momentum  = std::sqrt(kineticEnergy*kineticEnergy
    6467                            +2.0*restMass_c2*kineticEnergy);
    65 
    6668  G4ThreeVector pMomentum= momentum * pMomentumDirection;
    6769  SetCurvePnt( pPosition, pMomentum, curve_length );
    68   // Sets momentum direction as well.
     70    // Sets momentum direction as well.
    6971
    70   // Set the momentum direction again - keeping value from argument exactly
    7172  fMomentumDir=pMomentumDirection;
     73    // Set the momentum direction again - keeping value from argument exactly
    7274
    7375  InitialiseSpin( Spin );
    74 
    75   // fpChargeState = new G4ChargeState( charge, magnetic_dipole_moment );
    7676}
    7777
    … …  
    107107  else         Spin= *pSpin;
    108108  InitialiseSpin( Spin );
    109 
    110   // fpChargeState = new G4ChargeState( DBL_MAX );  //  charge not yet set !!
    111109}
    112110
    113111G4FieldTrack::G4FieldTrack( char )                  //  Nothing is set !!
    114   : fRestMass_c2(0.0), fLabTimeOfFlight(0.0),
    115    fChargeState( DBL_MAX ) //  charge not set
     112  : fKineticEnergy(0.), fRestMass_c2(0.), fLabTimeOfFlight(0.),
     113    fProperTimeOfFlight(0.), fChargeState( DBL_MAX )
    116114{
    117115  G4ThreeVector Zero(0.0, 0.0, 0.0);
    118116  SetCurvePnt( Zero, Zero, 0.0 );
    119117  InitialiseSpin( Zero );
    120   // fpChargeState = new G4ChargeState( DBL_MAX );   
    121118}
    122119

  • trunk/source/geometry/magneticfield/src/G4MagHelicalStepper.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagHelicalStepper.cc,v 1.23 2007/09/05 12:20:17 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagHelicalStepper.cc,v 1.24 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// --------------------------------------------------------------------
    … …  
    4545
    4646G4MagHelicalStepper::G4MagHelicalStepper(G4Mag_EqRhs *EqRhs)
    47    : G4MagIntegratorStepper(EqRhs, 6)  // integrate over 6 variables only !!
     47   : G4MagIntegratorStepper(EqRhs, 6), // integrate over 6 variables only !!
    4848                                       // position & velocity
    49 {
    50   fPtrMagEqOfMot = EqRhs;
     49     fPtrMagEqOfMot(EqRhs), fAngCurve(0.), frCurve(0.), frHelix(0.)
     50{
    5151}
    5252

  • trunk/source/geometry/magneticfield/src/G4MagIntegratorDriver.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4MagIntegratorDriver.cc,v 1.53 2009/11/05 22:31:43 japost Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4MagIntegratorDriver.cc,v 1.57 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    7878    fNoInitialSmallSteps(0), fDyerr_max(0.0), fDyerr_mx2(0.0),
    7979    fDyerrPos_smTot(0.0), fDyerrPos_lgTot(0.0), fDyerrVel_lgTot(0.0),
    80     fSumH_sm(0.0),   fSumH_lg(0.0),
     80    fSumH_sm(0.0), fSumH_lg(0.0),
    8181    fStatisticsVerboseLevel(statisticsVerbose)
    8282{ 
    … …  
    208208
    209209#ifdef G4DEBUG_FIELD
     210    G4double xSubStepStart= x;
    210211    for (i=0;i<nvar;i++)  { ySubStepStart[i] = y[i]; }
    211212    yFldTrkStart.LoadFromArray(y, fNoIntegrationVariables);
    … …  
    227228      lastStepSucceeded= (hdid == h);   
    228229#ifdef G4DEBUG_FIELD
    229       if (dbg>2)
    230       {
    231         PrintStatus( ySubStepStart, xSubStart, y, x, h,  nstp); // Only
     230      if (dbg>2) {
     231        PrintStatus( ySubStepStart, xSubStepStart, y, x, h,  nstp); // Only
    232232      }
    233233#endif
    … …  
    291291    {
    292292      if( nstp==nStpPr )  { G4cout << "***** Many steps ****" << G4endl; }
     293      G4cout << "MagIntDrv: " ;
    293294      G4cout << "hdid="  << std::setw(12) << hdid  << " "
    294              << "hnext=" << std::setw(12) << hnext << " " << G4endl;
     295             << "hnext=" << std::setw(12) << hnext << " "
     296             << "hstep=" << std::setw(12) << hstep << " (requested) "
     297             << G4endl;
    295298      PrintStatus( ystart, x1, y, x, h, (nstp==nStpPr) ? -nstp: nstp);
    296299    }
    … …  
    369372        lastStep = true;
    370373#ifdef G4DEBUG_FIELD
    371         if (dbg)
     374        if (dbg>2)
    372375        {
     376          int prec= G4cout.precision(12);
    373377          G4cout << "Warning: G4MagIntegratorDriver::AccurateAdvance"
    374378                 << G4endl
    375379                 << "  Integration step 'h' became "
    376                  << h << " due to roundoff " << G4endl
     380                 << h << " due to roundoff. " << G4endl
     381                 << " Calculated as difference of x2= "<< x2 << " and x=" << x
    377382                 << "  Forcing termination of advance." << G4endl;
     383          G4cout.precision(prec);
    378384        }         
    379385#endif
    … …  
    581587                 /  ( sqr(y[9]) + sqr(y[10]) + sqr(y[11]) );
    582588      errspin_sq *= inv_eps_vel_sq;
    583     }
     589      errmax_sq = std::max( errmax_sq, errspin_sq );
     590   }
    584591
    585592    if ( errmax_sq <= 1.0 )  { break; } // Step succeeded.

  • trunk/source/geometry/magneticfield/src/G4Mag_EqRhs.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_EqRhs.cc,v 1.11 2006/06/29 18:24:36 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_EqRhs.cc,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//  This is the standard right-hand side for equation of motion 
    … …  
    4848//
    4949G4Mag_EqRhs::G4Mag_EqRhs( G4MagneticField *magField )
    50    : G4EquationOfMotion(magField)
     50   : G4EquationOfMotion(magField), fCof_val(0.)
    5151{
    5252}

  • trunk/source/geometry/magneticfield/src/G4Mag_SpinEqRhs.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_SpinEqRhs.cc,v 1.15 2009/03/25 15:29:02 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_SpinEqRhs.cc,v 1.16 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// This is the standard right-hand side for equation of motion.
    … …  
    4343
    4444G4Mag_SpinEqRhs::G4Mag_SpinEqRhs( G4MagneticField* MagField )
    45   : G4Mag_EqRhs( MagField )
     45  : G4Mag_EqRhs( MagField ), omegac(0.), anomaly(0.0011659208),
     46    pcharge(0.), E(0.), gamma(0.), beta(0.)
    4647{
    47    anomaly = 0.0011659208;
    4848}
    4949
    … …  
    6262   omegac = 0.105658387*GeV/particleMass * 2.837374841e-3*(rad/cm/kilogauss);
    6363
    64    ParticleCharge = particleCharge;
     64   pcharge = particleCharge;
    6565
    6666   E = std::sqrt(sqr(MomentumXc)+sqr(particleMass));
    … …  
    101101   G4ThreeVector dSpin;
    102102
    103    dSpin = ParticleCharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
     103   dSpin = pcharge*omegac*(ucb*(Spin.cross(BField))-udb*(Spin.cross(u)));
    104104
    105105   dydx[ 9] = dSpin.x();

  • trunk/source/geometry/magneticfield/src/G4Mag_UsualEqRhs.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4Mag_UsualEqRhs.cc,v 1.12 2006/06/29 18:24:42 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4Mag_UsualEqRhs.cc,v 1.13 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    4040
    4141G4Mag_UsualEqRhs::G4Mag_UsualEqRhs( G4MagneticField* MagField )
    42   : G4Mag_EqRhs( MagField ) {}
     42  : G4Mag_EqRhs( MagField ), fInvCurrentMomentumXc(0.) {}
    4343
    4444G4Mag_UsualEqRhs::~G4Mag_UsualEqRhs() {}

  • trunk/source/geometry/magneticfield/src/G4NystromRK4.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4NystromRK4.cc,v 1.6 2009/12/16 17:49:57 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4NystromRK4.cc,v 1.9 2010/09/10 15:42:09 japost Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// History:
    … …  
    4949    m_cachedMom( false )
    5050{
    51   m_fldPosition[0]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
    52   m_fldPosition[1]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
    53   m_fldPosition[2]  = m_iPoint[1] = m_fPoint[2] = 9.9999999e+99 ;
     51  m_fldPosition[0]  = m_iPoint[0] = m_fPoint[0] = m_mPoint[0] = 9.9999999e+99 ;
     52  m_fldPosition[1]  = m_iPoint[1] = m_fPoint[1] = m_mPoint[1] = 9.9999999e+99 ;
     53  m_fldPosition[2]  = m_iPoint[2] = m_fPoint[2] = m_mPoint[2] = 9.9999999e+99 ;
    5454  m_fldPosition[3]  = -9.9999999e+99;
     55  m_lastField[0] = m_lastField[1] = m_lastField[2] = 0.0;
    5556
    5657  m_magdistance2 = distanceConstField*distanceConstField;
    … …  
    196197G4NystromRK4::ComputeRightHandSide(const G4double P[],G4double dPdS[])
    197198{
    198   getField(P);
     199  G4double P4vec[4]= { P[0], P[1], P[2], P[7] }; // Time is P[7]
     200  getField(P4vec);
    199201  m_mom   = std::sqrt(P[3]*P[3]+P[4]*P[4]+P[5]*P[5])     ;
    200202  m_imom  = 1./m_mom                                ;

  • trunk/source/geometry/magneticfield/src/G4RKG3_Stepper.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4RKG3_Stepper.cc,v 1.15 2007/08/21 10:17:41 tnikitin Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4RKG3_Stepper.cc,v 1.17 2010/07/23 14:13:49 tnikitin Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030// -------------------------------------------------------------------
    … …  
    3535
    3636G4RKG3_Stepper::G4RKG3_Stepper(G4Mag_EqRhs *EqRhs)
    37   : G4MagIntegratorStepper(EqRhs,6)
    38 {
    39  
    40   fPtrMagEqOfMot=EqRhs;
    41 
     37  : G4MagIntegratorStepper(EqRhs,6), hStep(0.), fPtrMagEqOfMot(EqRhs)
     38{
    4239}
    4340
    … …  
    4643}
    4744
    48 void G4RKG3_Stepper::Stepper(  const G4double  yInput[7],
    49                                const G4double dydx[7],
     45void G4RKG3_Stepper::Stepper(  const G4double yInput[8],
     46                               const G4double dydx[6],
    5047                                     G4double Step,
    51                                      G4double yOut[7],
     48                                     G4double yOut[8],
    5249                                     G4double yErr[])
    5350{
    … …  
    5754   G4double  by15 = 1. / 15. ; // was  0.066666666 ;
    5855
    59    G4double yTemp[7], dydxTemp[6], yIn[7] ;
     56   G4double yTemp[8], dydxTemp[6], yIn[8] ;
    6057   //  Saving yInput because yInput and yOut can be aliases for same array
    6158   for(i=0;i<nvar;i++) yIn[i]=yInput[i];
    62 
     59   yIn[6] = yInput[6];
     60   yIn[7] = yInput[7];
    6361   G4double h = Step * 0.5;
    6462   hStep=Step;
    … …  
    8280
    8381   h *= 2 ;
    84    StepNoErr(yIn,dydx,h,yTemp,B); // ,beTemp2) ;
     82   StepNoErr(yIn,dydx,h,yTemp,B);
    8583   for(i=0;i<nvar;i++)
    8684   {
    … …  
    127125// is passed from substep to substep.
    128126
    129 void G4RKG3_Stepper::StepNoErr(const G4double tIn[7],
    130                                const G4double dydx[7],
     127void G4RKG3_Stepper::StepNoErr(const G4double tIn[8],
     128                               const G4double dydx[6],
    131129                                     G4double Step,
    132                                      G4double tOut[7],
     130                                     G4double tOut[8],
    133131                                     G4double B[3]      )     // const
    134132   
    … …  
    137135  //  Copy and edit the routine above, to delete alpha2, beta2, ...
    138136   G4double K1[7],K2[7],K3[7],K4[7] ;
    139    G4double tTemp[7], yderiv[6] ;
     137   G4double tTemp[8], yderiv[6] ;
    140138
    141139  // Need Momentum value to give correct values to the coefficients in equation
    … …  
    188186      tOut[i+3] = tIn[i+3] + mom*(K1[i] + 2*K2[i] + 2*K3[i] +K4[i])*c3 ;
    189187   }
    190  
     188   tOut[6] = tIn[6];
     189   tOut[7] = tIn[7];
    191190   // NormaliseTangentVector( tOut );
    192191 

  • trunk/source/geometry/magneticfield/src/G4UniformElectricField.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UniformElectricField.cc,v 1.12 2006/06/29 18:24:56 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UniformElectricField.cc,v 1.13 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5252                                               G4double vPhi    )
    5353{
    54    if(vField >= 0 &&
    55       vTheta >= 0 && vTheta <= pi &&
    56       vPhi >= 0 && vPhi <= twopi)
    57    {
    58       fFieldComponents[0] = 0.0;
    59       fFieldComponents[1] = 0.0;
    60       fFieldComponents[2] = 0.0;
    61       fFieldComponents[3] = vField*std::sin(vTheta)*std::cos(vPhi) ;
    62       fFieldComponents[4] = vField*std::sin(vTheta)*std::sin(vPhi) ;
    63       fFieldComponents[5] = vField*std::cos(vTheta) ;
    64    }
    65    else
     54   if ( (vField<0) || (vTheta<0) || (vTheta>pi) || (vPhi<0) || (vPhi>twopi) )
    6655   {
    6756      G4Exception("G4UniformElectricField::G4UniformElectricField()",
    6857                  "WrongArgumentValue", FatalException, "Invalid parameters.");
    6958   }
     59
     60   fFieldComponents[0] = 0.0;
     61   fFieldComponents[1] = 0.0;
     62   fFieldComponents[2] = 0.0;
     63   fFieldComponents[3] = vField*std::sin(vTheta)*std::cos(vPhi) ;
     64   fFieldComponents[4] = vField*std::sin(vTheta)*std::sin(vPhi) ;
     65   fFieldComponents[5] = vField*std::cos(vTheta) ;
    7066}
    7167

  • trunk/source/geometry/magneticfield/src/G4UniformMagField.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4UniformMagField.cc,v 1.11 2006/06/29 18:24:58 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4UniformMagField.cc,v 1.12 2010/07/14 10:00:36 gcosmo Exp $
     28// GEANT4 tag $Name: field-V09-03-03 $
    2929//
    3030//
    … …  
    5656                                     G4double vPhi    )
    5757{
    58    if(vField >= 0 &&
    59       vTheta >= 0 && vTheta <= pi &&
    60       vPhi >= 0 && vPhi <= twopi)
    61    {
    62       fFieldComponents[0] = vField*std::sin(vTheta)*std::cos(vPhi) ;
    63       fFieldComponents[1] = vField*std::sin(vTheta)*std::sin(vPhi) ;
    64       fFieldComponents[2] = vField*std::cos(vTheta) ;
    65    }
    66    else
     58   if ( (vField<0) || (vTheta<0) || (vTheta>pi) || (vPhi<0) || (vPhi>twopi) )
    6759   {
    6860      G4Exception("G4UniformMagField::G4UniformMagField()",
    6961                  "WrongArgumentValue", FatalException, "Invalid parameters.") ;
    7062   }
     63   fFieldComponents[0] = vField*std::sin(vTheta)*std::cos(vPhi) ;
     64   fFieldComponents[1] = vField*std::sin(vTheta)*std::sin(vPhi) ;
     65   fFieldComponents[2] = vField*std::cos(vTheta) ;
    7166}
    7267

  • trunk/source/geometry/management/GNUmakefile

    r831 r1340  
    1 # $Id: GNUmakefile,v 1.5 2004/11/12 09:05:05 gcosmo Exp $
     1# $Id: GNUmakefile,v 1.6 2010/10/19 10:18:01 gcosmo Exp $
    22# -----------------------------------------------------------------------
    33# GNUmakefile for geometry/management library.  Gabriele Cosmo, 16/11/96.
    … …  
    1212include $(G4INSTALL)/config/architecture.gmk
    1313
    14 CPPFLAGS += -DG4GEOM_ALLOC_EXPORT
     14CPPFLAGS += -DG4ALLOC_EXPORT
    1515CPPFLAGS += -I$(G4BASE)/graphics_reps/include \
    1616            -I$(G4BASE)/intercoms/include \

  • trunk/source/geometry/management/History

    r1315 r1340  
    1 $Id: History,v 1.139 2010/02/09 16:50:44 gcosmo Exp $
     1$Id: History,v 1.145 2010/10/19 15:20:33 gcosmo Exp $
    22-------------------------------------------------------------------
    33
    … …  
    1717     * Reverse chronological order (last date on top), please *
    1818     ----------------------------------------------------------
     19
     20October 19, 2010  G.Cosmo                  geommng-V09-03-05
     21- Added Clone() virtual method to G4VSolid, to return a pointer to a
     22  dynamically allocated copy of the solid. To be used by Geant4-MT.
     23- Added copy-ctr and assignment operator to G4ReflectedSolid, now also
     24  implementing the Clone() method.
     25
     26September 6, 2010  G.Cosmo                 geommng-V09-03-04
     27- Fixed false-positive null forward case in G4SmartVoxelHeader.
     28
     29July 16, 2010  G.Cosmo                     geommng-V09-03-03
     30- Fixed compilation errors for code within G4GEOMETRY_VOXELDEBUG flag
     31  in G4GeometryManager and G4SmartVoxelHeader.
     32- Added printout of voxel limits in G4SmartVoxelHeader::Buildnodes() and
     33  removed protection on node width previously introduced.
     34
     35July 13, 2010  G.Cosmo                     geommng-V09-03-02
     36- G4SmartVoxelHeader: added protection against rare cases of zero computed
     37  node width in BuildNodes().
     38
     39July 2, 2010  G.Cosmo                      geommng-V09-03-01
     40- G4LogicalSurface: added missing initialisation for 'theTransRadSurface'
     41  data member. Made virtual destructor and constructors not inline.
     42- G4ErrorTarget: added missing default initialisation of a data member.
     43- G4SmartVoxelHeader: added return statement after exceptions in methods
     44  BuildVoxelsWithinLimits() and BuildNodes(); fixed case of invalid iterator
     45  in deletion of test slices in BuildVoxelsWithinLimits(); added protection
     46  for cases of invalide number of nodes in RefineNodes(), and added return
     47  statement after exceptions.
     48- G4GeometryManager: restore cout precision after printing voxel statistics.
     49- G4LogicalVolume: added return statement after exception in GetMass().
     50- G4ReflectedSolid: moved reserve() statement for 'vertices' within null-check
     51  condition in CalculateExtent() ...
    1952
    2053February 9, 2010  G.Cosmo                  geommng-V09-03-00

  • trunk/source/geometry/management/include/G4LogicalSurface.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4LogicalSurface.hh,v 1.11 2009/04/21 15:18:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4LogicalSurface.hh,v 1.12 2010/07/05 09:22:58 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
    … …  
    5353// Data members:
    5454//   G4String                       theName
    55 //   G4SurfaceProperty*              theSurfaceProperty
     55//   G4SurfaceProperty*             theSurfaceProperty
    5656//   G4TransitionRadiationSurface*  theTransRadSurface
    5757
    … …  
    8383 public:  // with description
    8484
    85    G4SurfaceProperty*  GetSurfaceProperty() const;
    86    void    SetSurfaceProperty(G4SurfaceProperty* ptrSurfaceProperty);
     85   inline G4SurfaceProperty*  GetSurfaceProperty() const;
     86   inline void SetSurfaceProperty(G4SurfaceProperty* ptrSurfaceProperty);
    8787
    88    const G4String& GetName() const;
    89    void    SetName(const G4String& name);
     88   inline const G4String& GetName() const;
     89   inline void SetName(const G4String& name);
    9090
    91    G4TransitionRadiationSurface* GetTransitionRadiationSurface() const;
    92    void SetTransitionRadiationSurface(G4TransitionRadiationSurface* tRadSurf);
     91   inline G4TransitionRadiationSurface* GetTransitionRadiationSurface() const;
     92   inline void SetTransitionRadiationSurface(G4TransitionRadiationSurface* tRadSurf);
    9393
    9494 public:  // without description
    … …  
    9696   virtual ~G4LogicalSurface();
    9797
    98    G4int operator==(const G4LogicalSurface &right) const;
    99    G4int operator!=(const G4LogicalSurface &right) const;
     98   inline G4int operator==(const G4LogicalSurface &right) const;
     99   inline G4int operator!=(const G4LogicalSurface &right) const;
    100100
    101101 protected:
    … …  
    103103   // There should be no instances of this class
    104104
    105    G4LogicalSurface(const G4String& name, G4SurfaceProperty* surfaceProperty);
     105   G4LogicalSurface(const G4String& name, G4SurfaceProperty* prop);
    106106     // Is the name more meaningful for the properties or the logical surface ?
    107107
    108  private:
     108 private:  // Copying restricted
    109109
    110    G4LogicalSurface(const G4LogicalSurface &right); // Copying restricted
    111    const G4LogicalSurface& operator=(const G4LogicalSurface& right);
     110   G4LogicalSurface(const G4LogicalSurface &right);
     111   inline const G4LogicalSurface& operator=(const G4LogicalSurface& right);
    112112
    113113 private:

  • trunk/source/geometry/management/include/G4LogicalSurface.icc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4LogicalSurface.icc,v 1.10 2009/04/21 15:18:15 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4LogicalSurface.icc,v 1.11 2010/07/05 09:22:58 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030////////////////////////////////////////////////////////////////////////
    … …  
    4242G4LogicalSurface::GetSurfaceProperty() const
    4343{
    44         return theSurfaceProperty;
     44  return theSurfaceProperty;
    4545}
    4646
    … …  
    4848G4LogicalSurface::SetSurfaceProperty(G4SurfaceProperty* ptrSurfaceProperty)
    4949{
    50         theSurfaceProperty = ptrSurfaceProperty;
     50  theSurfaceProperty = ptrSurfaceProperty;
    5151}
    5252
    … …  
    5454G4LogicalSurface::GetName() const
    5555{
    56         return theName;
     56  return theName;
    5757}
    5858
    … …  
    6060G4LogicalSurface::SetName(const G4String& name)
    6161{
    62         theName = name;
     62  theName = name;
    6363}
    6464
    … …  
    6666G4LogicalSurface::GetTransitionRadiationSurface() const
    6767{
    68         return theTransRadSurface;
     68  return theTransRadSurface;
    6969}
    7070
    … …  
    7373                                                transRadSurf )
    7474{
    75         theTransRadSurface= transRadSurf;
     75  theTransRadSurface= transRadSurf;
    7676}
    77 
    7877
    7978  //////////////
    … …  
    8483G4LogicalSurface::operator=(const G4LogicalSurface &right)
    8584{
    86         if (&right == this) return *this;
    87         theName = right.theName;
    88         theSurfaceProperty = right.theSurfaceProperty;
    89         theTransRadSurface = right.theTransRadSurface;
     85  if (&right == this)  { return *this; }
     86  theName = right.theName;
     87  theSurfaceProperty = right.theSurfaceProperty;
     88  theTransRadSurface = right.theTransRadSurface;
    9089       
    91         return *this;
     90  return *this;
    9291}
    9392
    … …  
    9594G4LogicalSurface::operator==(const G4LogicalSurface &right) const
    9695{
    97         return (this == (G4LogicalSurface *) &right);
     96  return (this == (G4LogicalSurface *) &right);
    9897}
    9998
    … …  
    101100G4LogicalSurface::operator!=(const G4LogicalSurface &right) const
    102101{
    103         return (this != (G4LogicalSurface *) &right);
     102  return (this != (G4LogicalSurface *) &right);
    104103}
    105 
    106   /////////////////
    107   // Constructors
    108   /////////////////
    109 
    110 inline G4LogicalSurface::G4LogicalSurface(const G4String&         name,
    111                                          G4SurfaceProperty* surfaceProperty)
    112   : theName(name), theSurfaceProperty(surfaceProperty)
    113 {
    114 }
    115 
    116 inline G4LogicalSurface::G4LogicalSurface(const G4LogicalSurface &right)
    117 {
    118         *this = right;
    119 }
    120 
    121 inline G4LogicalSurface::~G4LogicalSurface()
    122 {
    123 }

  • trunk/source/geometry/management/include/G4ReflectedSolid.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ReflectedSolid.hh,v 1.5 2006/06/29 18:31:03 gunter Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ReflectedSolid.hh,v 1.6 2010/10/19 15:20:18 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030//
    … …  
    9494    G4ThreeVector GetPointOnSurface() const;
    9595
     96    G4VSolid* Clone() const;
     97
    9698  public:  // with description
    9799
    … …  
    114116
    115117  public:  // without description
     118
     119    G4ReflectedSolid(const G4ReflectedSolid& rhs);
     120    G4ReflectedSolid& operator=(const G4ReflectedSolid& rhs);
     121      // Copy constructor and assignment operator.
    116122
    117123    void DescribeYourselfTo ( G4VGraphicsScene& scene ) const ;
    … …  
    149155    mutable G4Polyhedron* fpPolyhedron;  // Caches reflected G4Polyhedron.
    150156
    151   private:
    152 
    153     G4ReflectedSolid(const G4ReflectedSolid&);
    154     G4ReflectedSolid& operator=(const G4ReflectedSolid&);
    155       // Private copy constructor and assignment operator.
    156157} ;
    157158

  • trunk/source/geometry/management/include/G4VSolid.hh

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4VSolid.hh,v 1.29 2008/09/10 13:18:42 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4VSolid.hh,v 1.30 2010/10/19 15:19:37 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030//
    … …  
    193193      // Returns a random point located on the surface of the solid.
    194194
     195    virtual G4VSolid* Clone() const;
     196      // Returns a pointer of a dynamically allocated copy of the solid.
     197      // Returns NULL pointer with warning in case the concrete solid does not
     198      // implement this method. The caller has responsibility for ownership.
     199
    195200    virtual std::ostream& StreamInfo(std::ostream& os) const = 0;
    196201      // Dumps contents of the solid to a stream.

  • trunk/source/geometry/management/src/G4ErrorTarget.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ErrorTarget.cc,v 1.1 2007/05/16 12:50:52 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ErrorTarget.cc,v 1.2 2010/07/05 09:22:58 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030//
    … …  
    3535#include "G4ErrorTarget.hh"
    3636
    37 G4ErrorTarget::G4ErrorTarget(){}
    38 G4ErrorTarget::~G4ErrorTarget(){}
     37G4ErrorTarget::G4ErrorTarget()
     38 : theType(G4ErrorTarget_GeomVolume) {}
     39
     40G4ErrorTarget::~G4ErrorTarget() {}
    3941
    4042G4double G4ErrorTarget::GetDistanceFromPoint( const G4ThreeVector&,

  • trunk/source/geometry/management/src/G4GeometryManager.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4GeometryManager.cc,v 1.22 2008/05/16 13:46:48 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4GeometryManager.cc,v 1.24 2010/07/16 15:52:57 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030// class G4GeometryManager
    … …  
    261261#ifdef G4GEOMETRY_VOXELDEBUG
    262262     G4cout << "**** G4GeometryManager::BuildOptimisations" << G4endl
    263             << "     Skipping logical volume name = " << volume->GetName()
     263            << "     Skipping logical volume name = " << tVolume->GetName()
    264264            << G4endl;
    265265#endif
    … …  
    359359         << totalMemory/1024 << " kByte" << G4endl;
    360360  G4cout << "    Total CPU time elapsed for geometry optimisation: "
    361          << std::setprecision(2) << totalCpuTime << " seconds" << G4endl;
     361         << std::setprecision(2) << totalCpuTime << " seconds"
     362         << std::setprecision(6) << G4endl;
    362363 
    363364  //

  • trunk/source/geometry/management/src/G4LogicalVolume.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4LogicalVolume.cc,v 1.34 2009/09/24 13:22:57 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4LogicalVolume.cc,v 1.35 2010/07/05 09:22:58 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030//
    … …  
    211211    G4Exception("G4LogicalVolume::GetMass()", "InvalidSetup", FatalException,
    212212                "No material associated to the logical volume !");
     213    return 0;
    213214  }
    214215  if (!fSolid)
    … …  
    219220    G4Exception("G4LogicalVolume::GetMass()", "InvalidSetup", FatalException,
    220221                "No solid associated to the logical volume !");
     222    return 0;
    221223  }
    222224  G4double globalDensity = logMaterial->GetDensity();

  • trunk/source/geometry/management/src/G4ReflectedSolid.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4ReflectedSolid.cc,v 1.11 2006/11/08 09:56:33 gcosmo Exp $
    28 //
    29 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4ReflectedSolid.cc,v 1.13 2010/10/19 15:20:18 gcosmo Exp $
     28//
     29// GEANT4 tag $Name: geommng-V09-03-05 $
    3030//
    3131// Implementation for G4ReflectedSolid class for boolean
    … …  
    9393}
    9494
     95///////////////////////////////////////////////////////////////////
     96//
     97
     98G4ReflectedSolid::G4ReflectedSolid(const G4ReflectedSolid& rhs)
     99  : G4VSolid(rhs), fPtrSolid(rhs.fPtrSolid), fpPolyhedron(0)
     100{
     101  fPtrTransform      = new G4AffineTransform(*rhs.fPtrTransform);
     102  fDirectTransform   = new G4AffineTransform(*rhs.fDirectTransform);
     103  fPtrTransform3D    = new G4Transform3D(*rhs.fPtrTransform3D);
     104  fDirectTransform3D = new G4Transform3D(*rhs.fDirectTransform3D);
     105}
     106
     107///////////////////////////////////////////////////////////////////
     108//
     109
     110G4ReflectedSolid& G4ReflectedSolid::operator=(const G4ReflectedSolid& rhs)
     111{
     112  // Check assignment to self
     113  //
     114  if (this == &rhs)  { return *this; }
     115
     116  // Copy base class data
     117  //
     118  G4VSolid::operator=(rhs);
     119
     120  // Copy data
     121  //
     122  fPtrSolid= rhs.fPtrSolid; fpPolyhedron= 0;
     123  delete fPtrTransform;
     124  fPtrTransform= new G4AffineTransform(*rhs.fPtrTransform);
     125  delete fDirectTransform;
     126  fDirectTransform= new G4AffineTransform(*rhs.fDirectTransform);
     127  delete fPtrTransform3D;
     128  fPtrTransform3D= new G4Transform3D(*rhs.fPtrTransform3D);
     129  delete fDirectTransform3D;
     130  fDirectTransform3D= new G4Transform3D(*rhs.fDirectTransform3D);
     131
     132  return *this;
     133}
     134
     135///////////////////////////////////////////////////////////////////
     136//
     137
    95138G4GeometryType G4ReflectedSolid::GetEntityType() const
    96139{
    … …  
    257300  G4Point3D tmpPoint;
    258301
    259 // Calculate rotated vertex coordinates
     302  // Calculate rotated vertex coordinates
    260303
    261304  G4ThreeVectorList* vertices = new G4ThreeVectorList();
    262   vertices->reserve(8);
    263305
    264306  if (vertices)
    265307  {
     308    vertices->reserve(8);
     309
    266310    G4ThreeVector vertex0(x1,y1,z1) ;
    267311    tmpPoint    = transform3D*G4Point3D(vertex0);
    … …  
    494538{
    495539  DumpInfo();
    496   G4Exception("G4BooleanSolid::ComputeDimensions()",
     540  G4Exception("G4ReflectedSolid::ComputeDimensions()",
    497541               "NotApplicable", FatalException,
    498542               "Method not applicable in this context!");
    … …  
    511555  return G4ThreeVector(newPoint.x(),newPoint.y(),newPoint.z());
    512556}
     557
     558//////////////////////////////////////////////////////////////////////////
     559//
     560// Make a clone of this object
     561
     562G4VSolid* G4ReflectedSolid::Clone() const
     563{
     564  return new G4ReflectedSolid(*this);
     565}
     566
    513567
    514568//////////////////////////////////////////////////////////////////////////

  • trunk/source/geometry/management/src/G4SmartVoxelHeader.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4SmartVoxelHeader.cc,v 1.35 2010/02/09 16:50:25 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4SmartVoxelHeader.cc,v 1.39 2010/09/06 09:39:21 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030//
    … …  
    148148      }
    149149    }
    150       else
     150    else
    151151    {
    152152      dyingNode = fslices[node]->GetNode();
    … …  
    510510          pTestSlices->pop_back();
    511511          for (G4ProxyVector::iterator i=pTestSlices->begin();
    512                                        i!=pTestSlices->end(); i++)
     512                                       i!=pTestSlices->end(); )
    513513          {
    514514            if (*i==tmpProx)
    515515            {
    516               pTestSlices->erase(i); i--;
     516              i = pTestSlices->erase(i);
     517            }
     518            else
     519            {
     520              ++i;
    517521            }
    518522          }
    519523          if ( tmpProx ) { delete tmpProx; }
    520         } 
     524        }
    521525        delete pTestSlices;
    522526      }
    … …  
    534538                "InvalidSetup", FatalException,
    535539                "Cannot select more than 3 axis for optimisation.");
     540    return;
    536541  }
    537542
    … …  
    548553
    549554#ifdef G4GEOMETRY_VOXELDEBUG
    550   G4cout << G4endl << "     Selected axis = " << faxis << G4endl;
     555  G4cout << G4endl << "     Volume = " << pVolume->GetName()
     556         << G4endl << "     Selected axis = " << faxis << G4endl;
    551557  for (size_t islice=0; islice<fslices.size(); islice++)
    552558  {
    … …  
    808814      G4Exception("G4SmartVoxelHeader::BuildNodes()", "InvalidSetup",
    809815                  FatalException, "Missing parameterisation.");
     816      return 0;
    810817    }
    811818
    … …  
    864871    minExtents[nVol] = targetMinExtent;
    865872    maxExtents[nVol] = targetMaxExtent;
     873
     874#ifdef G4GEOMETRY_VOXELDEBUG
     875   G4cout << "---------------------------------------------------" << G4endl
     876          << "     Volume = " << pDaughter->GetName() << G4endl
     877          << " Min Extent = " << targetMinExtent << G4endl
     878          << " Max Extent = " << targetMaxExtent << G4endl
     879          << "---------------------------------------------------" << G4endl;
     880#endif
    866881
    867882    // Check not entirely outside mother when processing toplevel nodes
    … …  
    960975  G4double nodeWidth = (motherMaxExtent-motherMinExtent)/noNodes;
    961976
    962 // Create G4VoxelNodes. Will Add proxies before setting fslices
    963 //
     977  // Create G4VoxelNodes. Will Add proxies before setting fslices
     978  //
    964979  G4NodeVector* nodeList = new G4NodeVector();
    965   nodeList->reserve(noNodes);
    966980  if (!nodeList)
    967981  {
    … …  
    970984    G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
    971985                FatalException, "NodeList allocation error.");
    972   }
     986    return 0;
     987  }
     988  nodeList->reserve(noNodes);
     989
    973990  for (nNode=0; nNode<noNodes; nNode++)
    974991  {
    … …  
    981998      G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
    982999                  FatalException, "Node allocation error.");
     1000      return 0;
    9831001    }
    9841002    nodeList->push_back(pNode);
    … …  
    10261044  //
    10271045  G4ProxyVector* proxyList = new G4ProxyVector();
    1028   proxyList->reserve(noNodes);
    10291046  if (!proxyList)
    10301047  {
    … …  
    10331050    G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
    10341051                FatalException, "Proxy list allocation error.");
    1035   }
     1052    return 0;
     1053  }
     1054  proxyList->reserve(noNodes);
     1055
    10361056  //
    10371057  // Fill proxy List
    … …  
    10491069      G4Exception("G4SmartVoxelHeader::BuildNodes()", "FatalError",
    10501070                  FatalException, "Proxy node allocation failed.");
     1071      return 0;
    10511072    }
    10521073    proxyList->push_back(pProxyNode);
    … …  
    11961217        noContainedDaughters = targetNode->GetNoContained();
    11971218        targetList = new G4VolumeNosVector();
    1198         targetList->reserve(noContainedDaughters);
    11991219        if (!targetList)
    12001220        {
    … …  
    12051225                      "FatalError", FatalException,
    12061226                      "Target volume node list allocation error.");
    1207         }
     1227          return;
     1228        }
     1229        targetList->reserve(noContainedDaughters);
    12081230        for (i=0; i<noContainedDaughters; i++)
    12091231        {
    … …  
    12181240               << " - " << maxNo << " inclusive" << G4endl;
    12191241#endif
     1242        if (minNo > maxNo)    // Delete node and list to be replaced
     1243        {                     // and avoid further action ...
     1244          delete targetNode;
     1245          delete targetList;
     1246          return;
     1247        }
     1248
    12201249        // Delete node proxies at start of collected sets of nodes/headers
    12211250        //
    … …  
    12461275          G4Exception("G4SmartVoxelHeader::RefineNodes()", "FatalError",
    12471276                      FatalException, "Refined VoxelHeader allocation error.");
     1277          return;
    12481278        }
    12491279        replaceHeader->SetMinEquivalentSliceNo(minNo);
    12501280        replaceHeader->SetMaxEquivalentSliceNo(maxNo);
    12511281        replaceHeaderProxy = new G4SmartVoxelProxy(replaceHeader);
    1252         if (!replaceHeader)
     1282        if (!replaceHeaderProxy)
    12531283        {
    12541284          G4cerr << "ERROR - G4SmartVoxelHeader::RefineNodes()" << G4endl
    … …  
    12561286          G4Exception("G4SmartVoxelHeader::RefineNodes()", "FatalError",
    12571287                      FatalException, "Refined VoxelProxy allocation error.");
     1288          return;
    12581289        }
    12591290        for (replaceNo=minNo; replaceNo<=maxNo; replaceNo++)

  • trunk/source/geometry/management/src/G4VSolid.cc

    r1337 r1340  
    2525//
    2626//
    27 // $Id: G4VSolid.cc,v 1.39 2008/09/23 13:07:41 gcosmo Exp $
    28 // GEANT4 tag $Name: geant4-09-04-beta-01 $
     27// $Id: G4VSolid.cc,v 1.40 2010/10/19 15:19:37 gcosmo Exp $
     28// GEANT4 tag $Name: geommng-V09-03-05 $
    2929//
    3030// class G4VSolid
    … …  
    158158        JustWarning, "Not implemented for this solid ! Returning origin.");
    159159    return G4ThreeVector(0,0,0);
     160}
     161
     162//////////////////////////////////////////////////////////////////////////
     163//
     164// Dummy implementations ...
     165
     166const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
     167{ return 0; }
     168
     169G4VSolid* G4VSolid::GetConstituentSolid(G4int)
     170{ return 0; }
     171
     172const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
     173{ return 0; }
     174
     175G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
     176{ return 0; }
     177
     178////////////////////////////////////////////////////////////////
     179//
     180// Returns an estimation of the solid volume in internal units.
     181// The number of statistics and error accuracy is fixed.
     182// This method may be overloaded by derived classes to compute the
     183// exact geometrical quantity for solids where this is possible.
     184// or anyway to cache the computed value.
     185// This implementation does NOT cache the computed value.
     186
     187G4double G4VSolid::GetCubicVolume()
     188{
     189  G4int cubVolStatistics = 1000000;
     190  G4double cubVolEpsilon = 0.001;
     191  return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
     192}
     193
     194////////////////////////////////////////////////////////////////
     195//
     196// Calculate cubic volume based on Inside() method.
     197// Accuracy is limited by the second argument or the statistics
     198// expressed by the first argument.
     199// Implementation is courtesy of Vasiliki Despoina Mitsou,
     200// University of Athens.
     201
     202G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
     203{
     204  G4int iInside=0;
     205  G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
     206  G4bool yesno;
     207  G4ThreeVector p;
     208  EInside in;
     209
     210  // values needed for CalculateExtent signature
     211
     212  G4VoxelLimits limit;                // Unlimited
     213  G4AffineTransform origin;
     214
     215  // min max extents of pSolid along X,Y,Z
     216
     217  yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
     218  yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
     219  yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
     220
     221  // limits
     222
     223  if(nStat < 100)    nStat   = 100;
     224  if(epsilon > 0.01) epsilon = 0.01;
     225
     226  for(G4int i = 0; i < nStat; i++ )
     227  {
     228    px = minX+(maxX-minX)*G4UniformRand();
     229    py = minY+(maxY-minY)*G4UniformRand();
     230    pz = minZ+(maxZ-minZ)*G4UniformRand();
     231    p  = G4ThreeVector(px,py,pz);
     232    in = this->Inside(p);
     233    if(in != kOutside) iInside++;   
     234  }
     235  volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
     236  return volume;
     237}
     238
     239////////////////////////////////////////////////////////////////
     240//
     241// Returns an estimation of the solid surface area in internal units.
     242// The number of statistics and error accuracy is fixed.
     243// This method may be overloaded by derived classes to compute the
     244// exact geometrical quantity for solids where this is possible.
     245// or anyway to cache the computed value.
     246// This implementation does NOT cache the computed value.
     247
     248G4double G4VSolid::GetSurfaceArea()
     249{
     250  G4int stat = 1000000;
     251  G4double ell = -1.;
     252  return EstimateSurfaceArea(stat,ell);
     253}
     254
     255////////////////////////////////////////////////////////////////
     256//
     257// Estimate surface area based on Inside(), DistanceToIn(), and
     258// DistanceToOut() methods. Accuracy is limited by the statistics
     259// defined by the first argument. Implemented by Mikhail Kosov.
     260
     261G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
     262{
     263  G4int inside=0;
     264  G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
     265  G4bool yesno;
     266  G4ThreeVector p;
     267  EInside in;
     268
     269  // values needed for CalculateExtent signature
     270
     271  G4VoxelLimits limit;                // Unlimited
     272  G4AffineTransform origin;
     273
     274  // min max extents of pSolid along X,Y,Z
     275
     276  yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
     277  yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
     278  yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
     279
     280  // limits
     281
     282  if(nStat < 100) { nStat = 100; }
     283
     284  G4double dX=maxX-minX;
     285  G4double dY=maxY-minY;
     286  G4double dZ=maxZ-minZ;
     287  if(ell<=0.)          // Automatic definition of skin thickness
     288  {
     289    G4double minval=dX;
     290    if(dY<dX) { minval=dY; }
     291    if(dZ<minval) { minval=dZ; }
     292    ell=.01*minval;
     293  }
     294
     295  G4double dd=2*ell;
     296  minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
     297
     298  for(G4int i = 0; i < nStat; i++ )
     299  {
     300    px = minX+dX*G4UniformRand();
     301    py = minY+dY*G4UniformRand();
     302    pz = minZ+dZ*G4UniformRand();
     303    p  = G4ThreeVector(px,py,pz);
     304    in = this->Inside(p);
     305    if(in != kOutside)
     306    {
     307      if  (DistanceToOut(p)<ell) { inside++; }
     308    }
     309    else if(DistanceToIn(p)<ell) { inside++; }
     310  }
     311  // @@ The conformal correction can be upgraded
     312  surf = dX*dY*dZ*inside/dd/nStat;
     313  return surf;
     314}
     315
     316///////////////////////////////////////////////////////////////////////////
     317//
     318// Returns a pointer of a dynamically allocated copy of the solid.
     319// Returns NULL pointer with warning in case the concrete solid does not
     320// implement this method. The caller has responsibility for ownership.
     321//
     322
     323G4VSolid* G4VSolid::Clone() const
     324{
     325  G4String ErrMessage = "Clone() method not implemented for type: "
     326                      + GetEntityType() + "! Returning NULL pointer!";
     327  G4Exception("G4VSolid::Clone()", "NotImplemented",
     328              JustWarning, ErrMessage);
     329  return 0;
    160330}
    161331
    … …  
    447617}
    448618
    449 const G4VSolid* G4VSolid::GetConstituentSolid(G4int) const
    450 { return 0; }
    451 
    452 G4VSolid* G4VSolid::GetConstituentSolid(G4int)
    453 { return 0; }
    454 
    455 const G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr() const
    456 { return 0; }
    457 
    458 G4DisplacedSolid* G4VSolid::GetDisplacedSolidPtr()
    459 { return 0; }
    460 
    461619G4VisExtent G4VSolid::GetExtent () const
    462620{
    … …  
    491649  return 0;
    492650}
    493 
    494 ////////////////////////////////////////////////////////////////
    495 //
    496 // Returns an estimation of the solid volume in internal units.
    497 // The number of statistics and error accuracy is fixed.
    498 // This method may be overloaded by derived classes to compute the
    499 // exact geometrical quantity for solids where this is possible.
    500 // or anyway to cache the computed value.
    501 // This implementation does NOT cache the computed value.
    502 
    503 G4double G4VSolid::GetCubicVolume()
    504 {
    505   G4int cubVolStatistics = 1000000;
    506   G4double cubVolEpsilon = 0.001;
    507   return EstimateCubicVolume(cubVolStatistics, cubVolEpsilon);
    508 }
    509 
    510 ////////////////////////////////////////////////////////////////
    511 //
    512 // Calculate cubic volume based on Inside() method.
    513 // Accuracy is limited by the second argument or the statistics
    514 // expressed by the first argument.
    515 // Implementation is courtesy of Vasiliki Despoina Mitsou,
    516 // University of Athens.
    517 
    518 G4double G4VSolid::EstimateCubicVolume(G4int nStat, G4double epsilon) const
    519 {
    520   G4int iInside=0;
    521   G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,volume;
    522   G4bool yesno;
    523   G4ThreeVector p;
    524   EInside in;
    525 
    526   // values needed for CalculateExtent signature
    527 
    528   G4VoxelLimits limit;                // Unlimited
    529   G4AffineTransform origin;
    530 
    531   // min max extents of pSolid along X,Y,Z
    532 
    533   yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
    534   yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
    535   yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
    536 
    537   // limits
    538 
    539   if(nStat < 100)    nStat   = 100;
    540   if(epsilon > 0.01) epsilon = 0.01;
    541 
    542   for(G4int i = 0; i < nStat; i++ )
    543   {
    544     px = minX+(maxX-minX)*G4UniformRand();
    545     py = minY+(maxY-minY)*G4UniformRand();
    546     pz = minZ+(maxZ-minZ)*G4UniformRand();
    547     p  = G4ThreeVector(px,py,pz);
    548     in = this->Inside(p);
    549     if(in != kOutside) iInside++;   
    550   }
    551   volume = (maxX-minX)*(maxY-minY)*(maxZ-minZ)*iInside/nStat;
    552   return volume;
    553 }
    554 
    555 ////////////////////////////////////////////////////////////////
    556 //
    557 // Returns an estimation of the solid surface area in internal units.
    558 // The number of statistics and error accuracy is fixed.
    559 // This method may be overloaded by derived classes to compute the
    560 // exact geometrical quantity for solids where this is possible.
    561 // or anyway to cache the computed value.
    562 // This implementation does NOT cache the computed value.
    563 
    564 G4double G4VSolid::GetSurfaceArea()
    565 {
    566   G4int stat = 1000000;
    567   G4double ell = -1.;
    568   return EstimateSurfaceArea(stat,ell);
    569 }
    570 
    571 ////////////////////////////////////////////////////////////////
    572 //
    573 // Estimate surface area based on Inside(), DistanceToIn(), and
    574 // DistanceToOut() methods. Accuracy is limited by the statistics
    575 // defined by the first argument. Implemented by Mikhail Kosov.
    576 
    577 G4double G4VSolid::EstimateSurfaceArea(G4int nStat, G4double ell) const
    578 {
    579   G4int inside=0;
    580   G4double px,py,pz,minX,maxX,minY,maxY,minZ,maxZ,surf;
    581   G4bool yesno;
    582   G4ThreeVector p;
    583   EInside in;
    584 
    585   // values needed for CalculateExtent signature
    586 
    587   G4VoxelLimits limit;                // Unlimited
    588   G4AffineTransform origin;
    589 
    590   // min max extents of pSolid along X,Y,Z
    591 
    592   yesno = this->CalculateExtent(kXAxis,limit,origin,minX,maxX);
    593   yesno = this->CalculateExtent(kYAxis,limit,origin,minY,maxY);
    594   yesno = this->CalculateExtent(kZAxis,limit,origin,minZ,maxZ);
    595 
    596   // limits
    597 
    598   if(nStat < 100) { nStat = 100; }
    599 
    600   G4double dX=maxX-minX;
    601   G4double dY=maxY-minY;
    602   G4double dZ=maxZ-minZ;
    603   if(ell<=0.)          // Automatic definition of skin thickness
    604   {
    605     G4double minval=dX;
    606     if(dY<dX) { minval=dY; }
    607     if(dZ<minval) { minval=dZ; }
    608     ell=.01*minval;
    609   }
    610 
    611   G4double dd=2*ell;
    612   minX-=ell; minY-=ell; minZ-=ell; dX+=dd; dY+=dd; dZ+=dd;
    613 
    614   for(G4int i = 0; i < nStat; i++ )
    615   {
    616     px = minX+dX*G4UniformRand();
    617     py = minY+dY*G4UniformRand();
    618     pz = minZ+dZ*G4UniformRand();
    619     p  = G4ThreeVector(px,py,pz);
    620     in = this->Inside(p);
    621     if(in != kOutside)
    622     {
    623       if  (DistanceToOut(p)<ell) { inside++; }
    624     }
    625     else if(DistanceToIn(p)<ell) { inside++; }
    626   }
    627   // @@ The conformal correction can be upgraded
    628   surf = dX*dY*dZ*inside/dd/nStat;
    629   return surf;
    630 }
Note: See TracChangeset for help on using the changeset viewer.