| 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
27 #include "G4ChannelingECHARM.hh"
28 #include "G4PhysicsLinearVector.hh"
29 #include "G4Physics2DVector.hh"
30 #include "G4SystemOfUnits.hh"
31
32 G4ChannelingECHARM::G4ChannelingECHARM(const G4String& fileName,G4double vConversion):
33 fVectorEC(0),
34 fDistances{0.,0.,0.},
35 fPoints{0,0,0},
36 fMaximum(-DBL_MAX),
37 fMinimum(DBL_MAX){
38 fDistances[0] = 0;
39 fDistances[1] = 0;
40 fDistances[2] = 0;
41 fPoints[0] = 0;
42 fPoints[1] = 0;
43 fPoints[2] = 0;
44 fVectorEC2D = 0;
45 ReadFromECHARM(fileName,vConversion);
46 }
47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48
49 G4ChannelingECHARM::~G4ChannelingECHARM(){
50 delete(fVectorEC);
51 delete(fVectorEC2D);
52 }
53
54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
55
56 G4double G4ChannelingECHARM::GetEC(G4ThreeVector& vPosition){
57 G4double vX = vPosition.x();
58 if (vX < 0.0) {
59 vX += ((int( - vX / fDistances[0]) + 1.0 ) * fDistances[0]);
60 }
61 else if( vX > fDistances[0] ){
62 vX -= ( int( vX / fDistances[0]) * fDistances[0] );
63 }
64 if(fPoints[1]==1){
65 return fVectorEC->Value(vX);
66 }
67 else{
68 G4double vY = vPosition.y();
69 if (vY < 0.0) {
70 vY += ((int( - vY / fDistances[1]) + 1.0 ) * fDistances[1]);
71 }
72 else if( vY > fDistances[1] ){
73 vY -= ( int( vY / fDistances[1]) * fDistances[1] );
74 }
75 return fVectorEC2D->Value((vX),(vY));
76 }
77 }
78
79 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81 void G4ChannelingECHARM::ReadFromECHARM(const G4String& filename,
82 G4double vConversion){
83 std::ifstream vFileIn;
84 vFileIn.open(filename);
85
86 vFileIn >> fPoints[0] >> fPoints[1] >> fPoints[2];
87 vFileIn >> fDistances[0] >> fDistances[1] >> fDistances[2];
88
89 fDistances[0] *= CLHEP::meter;
90 fDistances[1] *= CLHEP::meter;
91 fDistances[2] *= CLHEP::meter;
92 fMaximum = -DBL_MAX;
93 fMinimum = +DBL_MAX;
94
95 if(fPoints[1]<1){
96 G4ExceptionDescription ed;
97 ed << "No Points not found !" << G4endl;
98 G4Exception("G4ChannelingECHARM::ReadFromECHARM(...)",
99 "G4ChannelingECHARM",
100 FatalException,
101 ed);
102 return;
103 }
104 else if(fPoints[1]==1){
105 fVectorEC = new G4PhysicsLinearVector(0,fDistances[0],fPoints[0]);
106 }
107 else{
108 fVectorEC2D = new G4Physics2DVector(fPoints[0],fPoints[1]);
109 }
110 G4double stepX = fDistances[0]/fPoints[0];
111 G4double stepY = fDistances[1]/fPoints[1];
112 for(G4int i1=0;i1<fPoints[1]; i1++){
113 if(fPoints[1]>1){
114 fVectorEC2D->PutY(i1,i1*stepY);
115 }
116 for(G4int i0=0;i0<fPoints[0]; i0++){
117 double vTempX;
118 vFileIn >> vTempX;
119
120 vTempX *= vConversion;
121 if(vTempX > fMaximum) {fMaximum = vTempX;}
122 if(vTempX < fMinimum) {fMinimum = vTempX;}
123 if(fPoints[1]==1){
124 fVectorEC->PutValue(i0,vTempX);
125 }
126 else{
127 fVectorEC2D->PutValue(i0,i1,vTempX);
128 fVectorEC2D->PutX(i0,i0*stepX);
129 }
130 }
131 }
132 G4cout << "G4ChannelingECHARM::ReadFromECHARM() - " << vConversion << " " << fPoints[0] << " " << fDistances[0] << " " << fPoints[1] << " " << fDistances[1] << " " << fMinimum << " " << fMaximum << G4endl;
133
134 vFileIn.close();
135
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139