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