Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // G4DriverReporter 27 // 28 // Implementation 29 // 30 // 31 // Authors: J.Apostolakis - January/March 2020 32 // ------------------------------------------------------------------- 33 34 #include "G4DriverReporter.hh" 35 36 // --------------------------------------------------------------------------- 37 38 void G4DriverReporter::PrintStatus( const G4double* StartArr, 39 G4double xstart, 40 const G4double* CurrentArr, 41 G4double xcurrent, 42 G4double requestStep, 43 unsigned int subStepNo, 44 unsigned int noIntegrationVariables 45 ) 46 // Potentially add as arguments: 47 // <dydx> - as Initial Force 48 // stepTaken(hdid) - last step taken 49 // nextStep (hnext) - proposal for size 50 { 51 G4FieldTrack StartFT(G4ThreeVector(0,0,0), 52 G4ThreeVector(0,0,0), 0., 0., 0., 0. ); 53 G4FieldTrack CurrentFT (StartFT); 54 55 StartFT.LoadFromArray( StartArr, noIntegrationVariables); 56 StartFT.SetCurveLength( xstart); 57 CurrentFT.LoadFromArray( CurrentArr, noIntegrationVariables); 58 CurrentFT.SetCurveLength( xcurrent ); 59 60 PrintStatus(StartFT, CurrentFT, requestStep, subStepNo ); 61 } 62 63 // --------------------------------------------------------------------------- 64 const G4int noPrecision = 8; 65 const G4int prec7= noPrecision+2; 66 const G4int prec8= noPrecision+3; 67 const G4int prec9= noPrecision+4; 68 69 void G4DriverReporter::PrintStatus(const G4FieldTrack& StartFT, 70 const G4FieldTrack& CurrentFT, 71 G4double requestStep, 72 unsigned int subStepNo) 73 { 74 G4int verboseLevel= 2; // fVerboseLevel; 75 G4long oldPrec= G4cout.precision(noPrecision); 76 // G4cout.setf(ios_base::fixed,ios_base::floatfield); 77 78 const G4ThreeVector StartPosition= StartFT.GetPosition(); 79 const G4ThreeVector StartUnitVelocity= StartFT.GetMomentumDir(); 80 const G4ThreeVector CurrentPosition= CurrentFT.GetPosition(); 81 const G4ThreeVector CurrentUnitVelocity= CurrentFT.GetMomentumDir(); 82 83 G4double DotStartCurrentVeloc= StartUnitVelocity.dot(CurrentUnitVelocity); 84 85 G4double step_len= CurrentFT.GetCurveLength() - StartFT.GetCurveLength(); 86 G4double subStepSize = step_len; 87 88 if( (subStepNo <= 1) || (verboseLevel > 3) ) 89 { 90 subStepNo = - subStepNo; // To allow printing banner 91 92 G4cout << "------------------------------------------------------------------" 93 << G4endl; 94 G4cout << std::setw( 6) << " " << std::setw( 25) 95 << " G4DriverReporter: Current Position and Direction" << " " 96 << G4endl; 97 G4cout << std::setw( 5) << "Step#" << " " 98 << std::setw( prec7) << "s-curve" << " " 99 << std::setw( prec9) << "X(mm)" << " " 100 << std::setw( prec9) << "Y(mm)" << " " 101 << std::setw( prec9) << "Z(mm)" << " " 102 << std::setw( prec8) << " N_x " << " " 103 << std::setw( prec8) << " N_y " << " " 104 << std::setw( prec8) << " N_z " << " " 105 << std::setw( 6) << " N^2-1 " << " " 106 << std::setw(10) << " N(0).N " << " " 107 << std::setw( 7) << "KinEner " << " " 108 << std::setw(12) << "Track-l" << " " // Add the Sub-step ?? 109 << std::setw(12) << "Step-len" << " " 110 << std::setw(12) << "Step-len" << " " 111 << std::setw( 9) << "ReqStep" << " " 112 << G4endl; 113 } 114 115 G4cout.precision(noPrecision); 116 117 if( (subStepNo <= 0) ) 118 { 119 PrintStat_Aux( StartFT, requestStep, 0., 120 0, 0.0, 1.0); 121 } 122 123 // if( verboseLevel <= 3 ) 124 { 125 G4cout.precision(noPrecision); 126 PrintStat_Aux( CurrentFT, requestStep, step_len, 127 subStepNo, subStepSize, DotStartCurrentVeloc ); 128 } 129 G4cout << "------------------------------------------------------------------" 130 << G4endl; 131 G4cout.precision(oldPrec); 132 } 133 134 // --------------------------------------------------------------------------- 135 136 void G4DriverReporter::PrintStat_Aux(const G4FieldTrack& aFieldTrack, 137 G4double requestStep, 138 G4double step_len, 139 G4int subStepNo, 140 G4double subStepSize, 141 G4double dotVeloc_StartCurr) 142 { 143 const G4ThreeVector Position = aFieldTrack.GetPosition(); 144 const G4ThreeVector UnitVelocity = aFieldTrack.GetMomentumDir(); 145 146 G4long oldprec= G4cout.precision(noPrecision); 147 148 if( subStepNo >= 0) 149 { 150 G4cout << std::setw( 5) << subStepNo << " "; 151 } 152 else 153 { 154 G4cout << std::setw( 5) << "Start" << " "; 155 } 156 G4double curveLen= aFieldTrack.GetCurveLength(); 157 G4cout << std::setw( 7) << curveLen; 158 // G4cout.precision(noPrecision); 159 G4cout << std::setw( prec9) << Position.x() << " " 160 << std::setw( prec9) << Position.y() << " " 161 << std::setw( prec9) << Position.z() << " " 162 << std::setw( prec8) << UnitVelocity.x() << " " 163 << std::setw( prec8) << UnitVelocity.y() << " " 164 << std::setw( prec8) << UnitVelocity.z() << " "; 165 G4cout.precision(3); 166 G4double unitMagDif = UnitVelocity.mag2()-1.0; 167 if( std::fabs( unitMagDif ) < 1.0e-15 ) { unitMagDif = 0.0; } 168 G4cout << std::setw( 8) << unitMagDif << " "; 169 G4cout.precision(6); 170 G4cout << std::setw(10) << dotVeloc_StartCurr << " "; 171 G4cout.precision(oldprec); 172 G4cout << std::setw( prec7) << aFieldTrack.GetKineticEnergy(); 173 G4cout << std::setw(12) << step_len << " "; 174 175 static G4ThreadLocal G4double oldCurveLength = 0.0; 176 static G4ThreadLocal G4double oldSubStepLength = 0.0; 177 static G4ThreadLocal G4int oldSubStepNo = -1; 178 179 G4double subStep_len = 0.0; 180 if( curveLen > oldCurveLength ) 181 { 182 subStep_len= curveLen - oldCurveLength; 183 } 184 else if (subStepNo == oldSubStepNo) 185 { 186 subStep_len= oldSubStepLength; 187 } 188 oldCurveLength= curveLen; 189 oldSubStepLength= subStep_len; 190 191 G4cout << std::setw(12) << subStep_len << " "; 192 G4cout << std::setw(12) << subStepSize << " "; 193 if( requestStep != -1.0 ) 194 { 195 G4cout << std::setw( prec9) << requestStep << " "; 196 } 197 else 198 { 199 G4cout << std::setw( prec9) << " InitialStep " << " "; 200 } 201 G4cout << G4endl; 202 } 203