Geant4_10
CCalRotationMatrixFactory.cc
Go to the documentation of this file.
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 //
27 // File: CCalRotationMatrixFactory.cc
28 // Description: CCalRotationFactory is a factory class to define all rotation
29 // matrices used in geometry building
31 #include <fstream>
32 #include <stdlib.h>
33 
35 
36 #include "CCalutils.hh"
37 
38 #include "G4SystemOfUnits.hh"
39 
40 //#define debug
41 //#define ddebug
42 
43 CCalRotationMatrixFactory * CCalRotationMatrixFactory::instance = 0;
44 G4String CCalRotationMatrixFactory::file="";
45 
47  if (rotfile=="" || rotfile==file)
48  return getInstance();
49  else if (file="") {
50  file=rotfile;
51  return getInstance();
52  } else {
53  G4cerr << "ERROR: Trying to get Rotation Matrices from " << rotfile
54  << " when previously were retrieved from " << file <<"." << G4endl;
55  return 0;
56  }
57 }
58 
59 
61  if (file=="") {
62  G4cerr << "ERROR: You haven't defined which file to use for materials in "
63  << "CCalRotationMatrixFactory::getInstance(G4String)" << G4endl;
64  return 0;
65  }
66 
67  if (instance==0) {
68  instance = new CCalRotationMatrixFactory;
69  return instance;
70  }
71  else
72  return instance;
73 }
74 
76  if (rotfile!=file && file!="") {
77  G4cerr << "ERROR: Trying to change Rotation Matrices file name to "
78  << rotfile << "." << G4endl;
79  G4cerr << " Using previous file: " << file << G4endl;
80  }
81  file=rotfile;
82 }
83 
86  for(i=theMatrices.begin(); i != theMatrices.end(); ++i) {
87  delete (*i).second;
88  };
89  theMatrices.clear();
90 }
91 
93  G4RotationMatrix* retrot=0;
94  //Rotation :NULL is no rotation so a null pointer is returned
95  if (rot != ":NULL") {
96  //retrot untouched if rot not found!!!
97  G4RotationMatrixTableIterator it = theMatrices.find(rot);
98  if (it != theMatrices.end())
99  retrot = (*it).second;
100  }
101 
102  return retrot;
103 }
104 
106  G4double th1,
107  G4double phi1,
108  G4double th2,
109  G4double phi2,
110  G4double th3,
111  G4double phi3){
112  G4double sinth1, sinth2, sinth3, costh1, costh2, costh3;
113  G4double sinph1, sinph2, sinph3, cosph1, cosph2, cosph3;
114  G4double TH1 = th1/deg, TH2 = th2/deg, TH3 = th3/deg;
115  G4double PH1 = phi1/deg, PH2 = phi2/deg, PH3 = phi3/deg;
116 
117  if (TH1 == 0.0 || TH1 == 360) {
118  sinth1 = 0.0; costh1 = 1.0;
119  } else if (TH1 == 90.0 || TH1 == -270) {
120  sinth1 = 1.0; costh1 = 0.0;
121  } else if (TH1 == 180.0 || TH1 == -180.0) {
122  sinth1 = 0.0; costh1 = -1.0;
123  } else if (TH1 == 270.0 || TH1 == -90.0) {
124  sinth1 = -1.0; costh1 = 0.0;
125  } else {
126  sinth1 = std::sin(th1); costh1 = std::cos(th1);
127  }
128 
129  if (TH2 == 0.0 || TH2 == 360) {
130  sinth2 = 0.0; costh2 = 1.0;
131  } else if (TH2 == 90.0 || TH2 == -270) {
132  sinth2 = 1.0; costh2 = 0.0;
133  } else if (TH2 == 180.0 || TH2 == -180.0) {
134  sinth2 = 0.0; costh2 = -1.0;
135  } else if (TH2 == 270.0 || TH2 == -90.0) {
136  sinth2 = -1.0; costh2 = 0.0;
137  } else {
138  sinth2 = std::sin(th2); costh2 = std::cos(th2);
139  }
140 
141  if (TH3 == 0.0 || TH3 == 360) {
142  sinth3 = 0.0; costh3 = 1.0;
143  } else if (TH3 == 90.0 || TH2 == -270) {
144  sinth3 = 1.0; costh3 = 0.0;
145  } else if (TH3 == 180.0 || TH3 == -180.0) {
146  sinth3 = 0.0; costh3 = -1.0;
147  } else if (TH3 == 270.0 || TH3 == -90.0) {
148  sinth3 = -1.0; costh3 = 0.0;
149  } else {
150  sinth3 = std::sin(th3); costh3 = std::cos(th3);
151  }
152 
153  if (PH1 == 0.0 || PH1 == 360) {
154  sinph1 = 0.0; cosph1 = 1.0;
155  } else if (PH1 == 90.0 || PH1 == -270) {
156  sinph1 = 1.0; cosph1 = 0.0;
157  } else if (PH1 == 180.0 || PH1 == -180.0) {
158  sinph1 = 0.0; cosph1 = -1.0;
159  } else if (PH1 == 270.0 || PH1 == -90.0) {
160  sinph1 = -1.0; cosph1 = 0.0;
161  } else {
162  sinph1 = std::sin(phi1); cosph1 = std::cos(phi1);
163  }
164 
165  if (PH2 == 0.0 || PH2 == 360) {
166  sinph2 = 0.0; cosph2 = 1.0;
167  } else if (PH2 == 90.0 || PH2 == -270) {
168  sinph2 = 1.0; cosph2 = 0.0;
169  } else if (PH2 == 180.0 || PH2 == -180.0) {
170  sinph2 = 0.0; cosph2 = -1.0;
171  } else if (PH2 == 270.0 || PH2 == -90.0) {
172  sinph2 = -1.0; cosph2 = 0.0;
173  } else {
174  sinph2 = std::sin(phi2); cosph2 = std::cos(phi2);
175  }
176 
177  if (PH3 == 0.0 || PH3 == 360) {
178  sinph3 = 0.0; cosph3 = 1.0;
179  } else if (PH3 == 90.0 || PH3 == -270) {
180  sinph3 = 1.0; cosph3 = 0.0;
181  } else if (PH3 == 180.0 || PH3 == -180.0) {
182  sinph3 = 0.0; cosph3 = -1.0;
183  } else if (PH3 == 270.0 || PH3 == -90.0) {
184  sinph3 = -1.0; cosph3 = 0.0;
185  } else {
186  sinph3 = std::sin(phi3); cosph3 = std::cos(phi3);
187  }
188 
189  //xprime axis coordinates
190  CLHEP::Hep3Vector xprime(sinth1*cosph1,sinth1*sinph1,costh1);
191  //yprime axis coordinates
192  CLHEP::Hep3Vector yprime(sinth2*cosph2,sinth2*sinph2,costh2);
193  //zprime axis coordinates
194  CLHEP::Hep3Vector zprime(sinth3*cosph3,sinth3*sinph3,costh3);
195 
196 #ifdef ddebug
197  G4cout << xprime << '\t'; G4cout << yprime << '\t'; G4cout << zprime << G4endl;
198 #endif
199  G4RotationMatrix *rotMat = new G4RotationMatrix();
200  rotMat->rotateAxes(xprime, yprime, zprime);
201  if (*rotMat == G4RotationMatrix()) {
202  // G4cerr << "WARNING: Matrix " << name << " will not be created as a rotation matrix."
203  // G4cerr << "WARNING: Matrix " << name << " is = identity matrix. It will not be created as a rotation matrix." << G4endl;
204  delete rotMat;
205  rotMat=0;
206  } else {
207  rotMat->invert();
208  theMatrices[name]=rotMat;
209 #ifdef ddebug
210  G4cout << *rotMat << G4endl;
211 #endif
212  }
213 
214  return rotMat;
215 }
216 
217 CCalRotationMatrixFactory::CCalRotationMatrixFactory():theMatrices(){
218 
219  G4String path = getenv("CCAL_GLOBALPATH");
220  G4cout << " ==> Opening file " << file << "..." << G4endl;
221  std::ifstream is;
222  bool ok = openGeomFile(is, path, file);
223  if (!ok) {
224  G4cerr << "ERROR: Could not open file " << file << " ... Exiting!" << G4endl;
225  exit(-1);
226  }
227 
229  // Find *DO ROTM
230  findDO(is, G4String("ROTM"));
231 
232  char rubish[256];
233  G4String name;
234 
235 #ifdef debug
236  G4cout << " ==> Reading Rotation Matrices... " << G4endl;
237  G4cout << " Name\tTheta1\tPhi1\tTheta2\tPhi2\tTheta3\tPhi3"<< G4endl;
238 #endif
239 
240  is >> name;
241  while(name!="*ENDDO") {
242  if (name.index("#.")==0) { //It is a comment.Skip line.
243  is.getline(rubish,256,'\n');
244  } else {
245 #ifdef debug
246  G4cout << " " << name <<'\t';
247 #endif
248  G4double th1, phi1, th2, phi2, th3, phi3;
249  //Get xprime axis angles
250  is >> th1 >> phi1;
251 #ifdef debug
252  G4cout << th1 << '\t' << phi1 << '\t';
253 #endif
254  //Get yprime axis angles
255  is >> th2 >> phi2;
256 #ifdef debug
257  G4cout << th2 << '\t' << phi2 << '\t';
258 #endif
259  //Get zprime axis angles
260  is >> th3 >> phi3;
261 #ifdef debug
262  G4cout << th3 << '\t' << phi3 << '\t';
263 #endif
264 
265  is.getline(rubish,256,'\n');
266 #ifdef debug
267  G4cout << rubish << G4endl;
268 #endif
269 
270  AddMatrix(name, th1*deg, phi1*deg, th2*deg, phi2*deg, th3*deg, phi3*deg);
271  }
272 
273  is >> name;
274  };
275 
276  is.close();
277  G4cout << " " << theMatrices.size() << " rotation matrices read in." << G4endl;
278 }
279 
280 
281 // 29-Jan-2004 A.R. : commented to avoid clashes with CLHEP.
282 // Streaming operators for rotation matrices are
283 // already defined in CLHEP::HepRotation.
284 // std::ostream& operator<<(std::ostream& os , const G4RotationMatrix & rot){
285 // // os << "( " << rot.xx() << tab << rot.xy() << tab << rot.xz() << " )" << G4endl;
286 // // os << "( " << rot.yx() << tab << rot.yy() << tab << rot.yz() << " )" << G4endl;
287 // // os << "( " << rot.zx() << tab << rot.zy() << tab << rot.zz() << " )" << G4endl;
288 //
289 // os << "["
290 // << rot.thetaX()/deg << tab << rot.phiX()/deg << tab
291 // << rot.thetaY()/deg << tab << rot.phiY()/deg << tab
292 // << rot.thetaZ()/deg << tab << rot.phiZ()/deg << "]"
293 // << G4endl;
294 //
295 // return os;
296 // }
bool openGeomFile(std::ifstream &is, const G4String &pathname, const G4String &filename)
Definition: CCalutils.cc:116
CLHEP::HepRotation G4RotationMatrix
HepRotation & rotateAxes(const Hep3Vector &newX, const Hep3Vector &newY, const Hep3Vector &newZ)
Definition: Rotation.cc:105
const XML_Char * name
Definition: expat.h:151
HepRotation & invert()
std::ifstream & findDO(std::ifstream &, const G4String &)
Definition: CCalutils.cc:72
G4RotationMatrix * AddMatrix(const G4String &name, G4double th1, G4double phi1, G4double th2, G4double phi2, G4double th3, G4double phi3)
std::map< G4String, G4RotationMatrixPtr, std::less< G4String > >::iterator G4RotationMatrixTableIterator
G4GLOB_DLL std::ostream G4cout
static CCalRotationMatrixFactory * getInstance()
G4RotationMatrix * findMatrix(const G4String &)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static void setFileName(const G4String &rotfile)
G4GLOB_DLL std::ostream G4cerr