| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 ////////////////////////////////////////////// 26 ///////////////////////////////////////////////////////////////////////////////
27 // File: CCalRotationMatrixFactory.cc 27 // File: CCalRotationMatrixFactory.cc
28 // Description: CCalRotationFactory is a facto 28 // Description: CCalRotationFactory is a factory class to define all rotation
29 // matrices used in geometry buil 29 // matrices used in geometry building
30 ////////////////////////////////////////////// 30 ///////////////////////////////////////////////////////////////////////////////
31 #include <fstream> <<
32 #include <stdlib.h> <<
33 <<
34 #include "CCalRotationMatrixFactory.hh" 31 #include "CCalRotationMatrixFactory.hh"
35 32
36 #include "CCalutils.hh" 33 #include "CCalutils.hh"
37 34
38 #include "G4SystemOfUnits.hh" << 35 #include <fstream>
>> 36 #include <stdlib.h>
39 37
40 //#define debug 38 //#define debug
41 //#define ddebug 39 //#define ddebug
42 40
43 CCalRotationMatrixFactory * CCalRotationMatrix 41 CCalRotationMatrixFactory * CCalRotationMatrixFactory::instance = 0;
44 G4String CCalRotationMatrixFactory::file=""; 42 G4String CCalRotationMatrixFactory::file="";
45 43
46 CCalRotationMatrixFactory* CCalRotationMatrixF 44 CCalRotationMatrixFactory* CCalRotationMatrixFactory::getInstance(const G4String & rotfile){
47 if (rotfile=="" || rotfile==file) 45 if (rotfile=="" || rotfile==file)
48 return getInstance(); 46 return getInstance();
49 else if (file=="") { << 47 else if (file="") {
50 file=rotfile; 48 file=rotfile;
51 return getInstance(); 49 return getInstance();
52 } else { 50 } else {
53 G4cerr << "ERROR: Trying to get Rotation M 51 G4cerr << "ERROR: Trying to get Rotation Matrices from " << rotfile
54 << " when previously were retrieved f << 52 << " when previously were retrieved from " << file <<"." << G4endl;
55 return 0; 53 return 0;
56 } 54 }
57 } 55 }
58 56
59 57
60 CCalRotationMatrixFactory* CCalRotationMatrixF 58 CCalRotationMatrixFactory* CCalRotationMatrixFactory::getInstance(){
61 if (file=="") { 59 if (file=="") {
62 G4cerr << "ERROR: You haven't defined whic 60 G4cerr << "ERROR: You haven't defined which file to use for materials in "
63 << "CCalRotationMatrixFactory::getIns << 61 << "CCalRotationMatrixFactory::getInstance(G4String)" << G4endl;
64 return 0; 62 return 0;
65 } 63 }
66 64
67 if (instance==0) { 65 if (instance==0) {
68 instance = new CCalRotationMatrixFactory; 66 instance = new CCalRotationMatrixFactory;
69 return instance; 67 return instance;
70 } 68 }
71 else 69 else
72 return instance; 70 return instance;
73 } 71 }
74 72
75 void CCalRotationMatrixFactory::setFileName(co 73 void CCalRotationMatrixFactory::setFileName(const G4String& rotfile) {
76 if (rotfile!=file && file!="") { 74 if (rotfile!=file && file!="") {
77 G4cerr << "ERROR: Trying to change Rotatio 75 G4cerr << "ERROR: Trying to change Rotation Matrices file name to "
78 << rotfile << "." << G4endl; << 76 << rotfile << "." << G4endl;
79 G4cerr << " Using previous file: " < 77 G4cerr << " Using previous file: " << file << G4endl;
80 } 78 }
81 file=rotfile; 79 file=rotfile;
82 } 80 }
83 81
84 CCalRotationMatrixFactory::~CCalRotationMatrix 82 CCalRotationMatrixFactory::~CCalRotationMatrixFactory(){
85 G4RotationMatrixTableIterator i; 83 G4RotationMatrixTableIterator i;
86 for(i=theMatrices.begin(); i != theMatrices. 84 for(i=theMatrices.begin(); i != theMatrices.end(); ++i) {
87 delete (*i).second; 85 delete (*i).second;
88 }; 86 };
89 theMatrices.clear(); 87 theMatrices.clear();
90 } 88 }
91 89
92 G4RotationMatrix* CCalRotationMatrixFactory::f 90 G4RotationMatrix* CCalRotationMatrixFactory::findMatrix(const G4String & rot) {
93 G4RotationMatrix* retrot=0; 91 G4RotationMatrix* retrot=0;
94 //Rotation :NULL is no rotation so a null po 92 //Rotation :NULL is no rotation so a null pointer is returned
95 if (rot != ":NULL") { 93 if (rot != ":NULL") {
96 //retrot untouched if rot not found!!! 94 //retrot untouched if rot not found!!!
97 G4RotationMatrixTableIterator it = theMatr 95 G4RotationMatrixTableIterator it = theMatrices.find(rot);
98 if (it != theMatrices.end()) 96 if (it != theMatrices.end())
99 retrot = (*it).second; 97 retrot = (*it).second;
100 } 98 }
101 99
102 return retrot; //!!!Maybe a treatment on not 100 return retrot; //!!!Maybe a treatment on not-found case needed.
103 } 101 }
104 102
105 G4RotationMatrix* CCalRotationMatrixFactory::A 103 G4RotationMatrix* CCalRotationMatrixFactory::AddMatrix(const G4String& name,
106 << 104 G4double th1,
107 << 105 G4double phi1,
108 << 106 G4double th2,
109 << 107 G4double phi2,
110 << 108 G4double th3,
111 << 109 G4double phi3){
112 G4double sinth1, sinth2, sinth3, costh1, co 110 G4double sinth1, sinth2, sinth3, costh1, costh2, costh3;
113 G4double sinph1, sinph2, sinph3, cosph1, cos 111 G4double sinph1, sinph2, sinph3, cosph1, cosph2, cosph3;
114 G4double TH1 = th1/deg, TH2 = th2/deg, TH3 = 112 G4double TH1 = th1/deg, TH2 = th2/deg, TH3 = th3/deg;
115 G4double PH1 = phi1/deg, PH2 = phi2/deg, PH3 113 G4double PH1 = phi1/deg, PH2 = phi2/deg, PH3 = phi3/deg;
116 << 114
117 if (TH1 == 0.0 || TH1 == 360) { 115 if (TH1 == 0.0 || TH1 == 360) {
118 sinth1 = 0.0; costh1 = 1.0; 116 sinth1 = 0.0; costh1 = 1.0;
119 } else if (TH1 == 90.0 || TH1 == -270) { 117 } else if (TH1 == 90.0 || TH1 == -270) {
120 sinth1 = 1.0; costh1 = 0.0; 118 sinth1 = 1.0; costh1 = 0.0;
121 } else if (TH1 == 180.0 || TH1 == -180.0) { 119 } else if (TH1 == 180.0 || TH1 == -180.0) {
122 sinth1 = 0.0; costh1 = -1.0; 120 sinth1 = 0.0; costh1 = -1.0;
123 } else if (TH1 == 270.0 || TH1 == -90.0) { 121 } else if (TH1 == 270.0 || TH1 == -90.0) {
124 sinth1 = -1.0; costh1 = 0.0; 122 sinth1 = -1.0; costh1 = 0.0;
125 } else { 123 } else {
126 sinth1 = std::sin(th1); costh1 = std::cos( 124 sinth1 = std::sin(th1); costh1 = std::cos(th1);
127 } 125 }
128 126
129 if (TH2 == 0.0 || TH2 == 360) { 127 if (TH2 == 0.0 || TH2 == 360) {
130 sinth2 = 0.0; costh2 = 1.0; 128 sinth2 = 0.0; costh2 = 1.0;
131 } else if (TH2 == 90.0 || TH2 == -270) { 129 } else if (TH2 == 90.0 || TH2 == -270) {
132 sinth2 = 1.0; costh2 = 0.0; 130 sinth2 = 1.0; costh2 = 0.0;
133 } else if (TH2 == 180.0 || TH2 == -180.0) { 131 } else if (TH2 == 180.0 || TH2 == -180.0) {
134 sinth2 = 0.0; costh2 = -1.0; 132 sinth2 = 0.0; costh2 = -1.0;
135 } else if (TH2 == 270.0 || TH2 == -90.0) { 133 } else if (TH2 == 270.0 || TH2 == -90.0) {
136 sinth2 = -1.0; costh2 = 0.0; 134 sinth2 = -1.0; costh2 = 0.0;
137 } else { 135 } else {
138 sinth2 = std::sin(th2); costh2 = std::cos( 136 sinth2 = std::sin(th2); costh2 = std::cos(th2);
139 } 137 }
140 << 138
141 if (TH3 == 0.0 || TH3 == 360) { 139 if (TH3 == 0.0 || TH3 == 360) {
142 sinth3 = 0.0; costh3 = 1.0; 140 sinth3 = 0.0; costh3 = 1.0;
143 } else if (TH3 == 90.0 || TH2 == -270) { 141 } else if (TH3 == 90.0 || TH2 == -270) {
144 sinth3 = 1.0; costh3 = 0.0; 142 sinth3 = 1.0; costh3 = 0.0;
145 } else if (TH3 == 180.0 || TH3 == -180.0) { 143 } else if (TH3 == 180.0 || TH3 == -180.0) {
146 sinth3 = 0.0; costh3 = -1.0; 144 sinth3 = 0.0; costh3 = -1.0;
147 } else if (TH3 == 270.0 || TH3 == -90.0) { 145 } else if (TH3 == 270.0 || TH3 == -90.0) {
148 sinth3 = -1.0; costh3 = 0.0; 146 sinth3 = -1.0; costh3 = 0.0;
149 } else { 147 } else {
150 sinth3 = std::sin(th3); costh3 = std::cos( 148 sinth3 = std::sin(th3); costh3 = std::cos(th3);
151 } 149 }
152 150
153 if (PH1 == 0.0 || PH1 == 360) { 151 if (PH1 == 0.0 || PH1 == 360) {
154 sinph1 = 0.0; cosph1 = 1.0; 152 sinph1 = 0.0; cosph1 = 1.0;
155 } else if (PH1 == 90.0 || PH1 == -270) { 153 } else if (PH1 == 90.0 || PH1 == -270) {
156 sinph1 = 1.0; cosph1 = 0.0; 154 sinph1 = 1.0; cosph1 = 0.0;
157 } else if (PH1 == 180.0 || PH1 == -180.0) { 155 } else if (PH1 == 180.0 || PH1 == -180.0) {
158 sinph1 = 0.0; cosph1 = -1.0; 156 sinph1 = 0.0; cosph1 = -1.0;
159 } else if (PH1 == 270.0 || PH1 == -90.0) { 157 } else if (PH1 == 270.0 || PH1 == -90.0) {
160 sinph1 = -1.0; cosph1 = 0.0; 158 sinph1 = -1.0; cosph1 = 0.0;
161 } else { 159 } else {
162 sinph1 = std::sin(phi1); cosph1 = std::cos 160 sinph1 = std::sin(phi1); cosph1 = std::cos(phi1);
163 } 161 }
164 162
165 if (PH2 == 0.0 || PH2 == 360) { 163 if (PH2 == 0.0 || PH2 == 360) {
166 sinph2 = 0.0; cosph2 = 1.0; 164 sinph2 = 0.0; cosph2 = 1.0;
167 } else if (PH2 == 90.0 || PH2 == -270) { 165 } else if (PH2 == 90.0 || PH2 == -270) {
168 sinph2 = 1.0; cosph2 = 0.0; 166 sinph2 = 1.0; cosph2 = 0.0;
169 } else if (PH2 == 180.0 || PH2 == -180.0) { 167 } else if (PH2 == 180.0 || PH2 == -180.0) {
170 sinph2 = 0.0; cosph2 = -1.0; 168 sinph2 = 0.0; cosph2 = -1.0;
171 } else if (PH2 == 270.0 || PH2 == -90.0) { 169 } else if (PH2 == 270.0 || PH2 == -90.0) {
172 sinph2 = -1.0; cosph2 = 0.0; 170 sinph2 = -1.0; cosph2 = 0.0;
173 } else { 171 } else {
174 sinph2 = std::sin(phi2); cosph2 = std::cos 172 sinph2 = std::sin(phi2); cosph2 = std::cos(phi2);
175 } 173 }
176 << 174
177 if (PH3 == 0.0 || PH3 == 360) { 175 if (PH3 == 0.0 || PH3 == 360) {
178 sinph3 = 0.0; cosph3 = 1.0; 176 sinph3 = 0.0; cosph3 = 1.0;
179 } else if (PH3 == 90.0 || PH3 == -270) { 177 } else if (PH3 == 90.0 || PH3 == -270) {
180 sinph3 = 1.0; cosph3 = 0.0; 178 sinph3 = 1.0; cosph3 = 0.0;
181 } else if (PH3 == 180.0 || PH3 == -180.0) { 179 } else if (PH3 == 180.0 || PH3 == -180.0) {
182 sinph3 = 0.0; cosph3 = -1.0; 180 sinph3 = 0.0; cosph3 = -1.0;
183 } else if (PH3 == 270.0 || PH3 == -90.0) { 181 } else if (PH3 == 270.0 || PH3 == -90.0) {
184 sinph3 = -1.0; cosph3 = 0.0; 182 sinph3 = -1.0; cosph3 = 0.0;
185 } else { 183 } else {
186 sinph3 = std::sin(phi3); cosph3 = std::cos 184 sinph3 = std::sin(phi3); cosph3 = std::cos(phi3);
187 } 185 }
188 << 186
189 //xprime axis coordinates 187 //xprime axis coordinates
190 CLHEP::Hep3Vector xprime(sinth1*cosph1,sinth 188 CLHEP::Hep3Vector xprime(sinth1*cosph1,sinth1*sinph1,costh1);
191 //yprime axis coordinates 189 //yprime axis coordinates
192 CLHEP::Hep3Vector yprime(sinth2*cosph2,sinth 190 CLHEP::Hep3Vector yprime(sinth2*cosph2,sinth2*sinph2,costh2);
193 //zprime axis coordinates 191 //zprime axis coordinates
194 CLHEP::Hep3Vector zprime(sinth3*cosph3,sinth 192 CLHEP::Hep3Vector zprime(sinth3*cosph3,sinth3*sinph3,costh3);
195 193
196 #ifdef ddebug 194 #ifdef ddebug
197 G4cout << xprime << '\t'; G4cout << yprim 195 G4cout << xprime << '\t'; G4cout << yprime << '\t'; G4cout << zprime << G4endl;
198 #endif 196 #endif
199 G4RotationMatrix *rotMat = new G4RotationMat 197 G4RotationMatrix *rotMat = new G4RotationMatrix();
200 rotMat->rotateAxes(xprime, yprime, zprime); 198 rotMat->rotateAxes(xprime, yprime, zprime);
201 if (*rotMat == G4RotationMatrix()) { 199 if (*rotMat == G4RotationMatrix()) {
202 // G4cerr << "WARNING: Matrix " << name << 200 // G4cerr << "WARNING: Matrix " << name << " will not be created as a rotation matrix."
203 // G4cerr << "WARNING: Matrix " << name << 201 // G4cerr << "WARNING: Matrix " << name << " is = identity matrix. It will not be created as a rotation matrix." << G4endl;
204 delete rotMat; 202 delete rotMat;
205 rotMat=0; 203 rotMat=0;
206 } else { 204 } else {
207 rotMat->invert(); 205 rotMat->invert();
208 theMatrices[name]=rotMat; 206 theMatrices[name]=rotMat;
209 #ifdef ddebug 207 #ifdef ddebug
210 G4cout << *rotMat << G4endl; 208 G4cout << *rotMat << G4endl;
211 #endif 209 #endif
212 } 210 }
213 211
214 return rotMat; 212 return rotMat;
215 } 213 }
216 214
217 CCalRotationMatrixFactory::CCalRotationMatrixF 215 CCalRotationMatrixFactory::CCalRotationMatrixFactory():theMatrices(){
218 << 216
219 G4String path = "NULL"; << 217 G4String path = getenv("CCAL_GLOBALPATH");
220 if (std::getenv("CCAL_GLOBALPATH")) <<
221 path = std::getenv("CCAL_GLOBALPATH"); <<
222 <<
223 G4cout << " ==> Opening file " << file << ". 218 G4cout << " ==> Opening file " << file << "..." << G4endl;
224 std::ifstream is; 219 std::ifstream is;
225 G4bool ok = openGeomFile(is, path, file); << 220 bool ok = openGeomFile(is, path, file);
226 if (!ok) { 221 if (!ok) {
227 G4ExceptionDescription ed; << 222 G4cerr << "ERROR: Could not open file " << file << " ... Exiting!" << G4endl;
228 ed << "Could not open file " << file << " << 223 exit(-1);
229 G4Exception("CCalRotationMatrixFactory::CC <<
230 "ccal002", <<
231 FatalException,ed); <<
232 } 224 }
233 225
234 //////////////////////////////////////////// 226 //////////////////////////////////////////////////
235 // Find *DO ROTM 227 // Find *DO ROTM
236 findDO(is, G4String("ROTM")); 228 findDO(is, G4String("ROTM"));
237 229
238 char rubish[256]; 230 char rubish[256];
239 G4String name; 231 G4String name;
240 232
241 #ifdef debug 233 #ifdef debug
242 G4cout << " ==> Reading Rotation Matrice 234 G4cout << " ==> Reading Rotation Matrices... " << G4endl;
243 G4cout << " Name\tTheta1\tPhi1\tTheta2 235 G4cout << " Name\tTheta1\tPhi1\tTheta2\tPhi2\tTheta3\tPhi3"<< G4endl;
244 #endif 236 #endif
245 237
246 is >> name; 238 is >> name;
247 while(name!="*ENDDO") { 239 while(name!="*ENDDO") {
248 if (name.find("#.")==0) { //It is a commen << 240 if (name.index("#.")==0) { //It is a comment.Skip line.
249 is.getline(rubish,256,'\n'); 241 is.getline(rubish,256,'\n');
250 } else { 242 } else {
251 #ifdef debug 243 #ifdef debug
252 G4cout << " " << name <<'\t'; 244 G4cout << " " << name <<'\t';
253 #endif 245 #endif
254 G4double th1, phi1, th2, phi2, th3, phi3 246 G4double th1, phi1, th2, phi2, th3, phi3;
255 //Get xprime axis angles 247 //Get xprime axis angles
256 is >> th1 >> phi1; 248 is >> th1 >> phi1;
257 #ifdef debug 249 #ifdef debug
258 G4cout << th1 << '\t' << phi1 << '\t'; 250 G4cout << th1 << '\t' << phi1 << '\t';
259 #endif 251 #endif
260 //Get yprime axis angles 252 //Get yprime axis angles
261 is >> th2 >> phi2; 253 is >> th2 >> phi2;
262 #ifdef debug 254 #ifdef debug
263 G4cout << th2 << '\t' << phi2 << '\t'; 255 G4cout << th2 << '\t' << phi2 << '\t';
264 #endif 256 #endif
265 //Get zprime axis angles 257 //Get zprime axis angles
266 is >> th3 >> phi3; 258 is >> th3 >> phi3;
267 #ifdef debug 259 #ifdef debug
268 G4cout << th3 << '\t' << phi3 << '\t'; 260 G4cout << th3 << '\t' << phi3 << '\t';
269 #endif 261 #endif
270 262
271 is.getline(rubish,256,'\n'); 263 is.getline(rubish,256,'\n');
272 #ifdef debug 264 #ifdef debug
273 G4cout << rubish << G4endl; 265 G4cout << rubish << G4endl;
274 #endif 266 #endif
275 267
276 AddMatrix(name, th1*deg, phi1*deg, th2*d 268 AddMatrix(name, th1*deg, phi1*deg, th2*deg, phi2*deg, th3*deg, phi3*deg);
277 } 269 }
278 270
279 is >> name; 271 is >> name;
280 }; 272 };
281 273
282 is.close(); 274 is.close();
283 G4cout << " " << theMatrices.size() < 275 G4cout << " " << theMatrices.size() << " rotation matrices read in." << G4endl;
284 } 276 }
285 277
286 278
287 // 29-Jan-2004 A.R. : commented to avoid clash 279 // 29-Jan-2004 A.R. : commented to avoid clashes with CLHEP.
288 // Streaming operators for 280 // Streaming operators for rotation matrices are
289 // already defined in CLHEP 281 // already defined in CLHEP::HepRotation.
290 // std::ostream& operator<<(std::ostream& os , 282 // std::ostream& operator<<(std::ostream& os , const G4RotationMatrix & rot){
291 // // os << "( " << rot.xx() << tab << rot. 283 // // os << "( " << rot.xx() << tab << rot.xy() << tab << rot.xz() << " )" << G4endl;
292 // // os << "( " << rot.yx() << tab << rot. 284 // // os << "( " << rot.yx() << tab << rot.yy() << tab << rot.yz() << " )" << G4endl;
293 // // os << "( " << rot.zx() << tab << rot. 285 // // os << "( " << rot.zx() << tab << rot.zy() << tab << rot.zz() << " )" << G4endl;
294 // 286 //
295 // os << "[" 287 // os << "["
296 // << rot.thetaX()/deg << tab << rot.phiX 288 // << rot.thetaX()/deg << tab << rot.phiX()/deg << tab
297 // << rot.thetaY()/deg << tab << rot.phiY 289 // << rot.thetaY()/deg << tab << rot.phiY()/deg << tab
298 // << rot.thetaZ()/deg << tab << rot.phiZ 290 // << rot.thetaZ()/deg << tab << rot.phiZ()/deg << "]"
299 // << G4endl; 291 // << G4endl;
300 // 292 //
301 // return os; 293 // return os;
302 // } 294 // }
303 295