composite_calorimeter/src/CCalRotationMatrixFactory.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 ///////////////////////////////////////////////////////////////////////////////
0027 // File: CCalRotationMatrixFactory.cc
0028 // Description: CCalRotationFactory is a factory class to define all rotation
0029 //              matrices used in geometry building
0030 ///////////////////////////////////////////////////////////////////////////////
0031 #include <fstream>
0032 #include <stdlib.h>
0033
0034 #include "CCalRotationMatrixFactory.hh"
0035
0036 #include "CCalutils.hh"
0037
0038 #include "G4SystemOfUnits.hh"
0039
0040 //#define debug
0041 //#define ddebug
0042
0043 CCalRotationMatrixFactory * CCalRotationMatrixFactory::instance = 0;
0044 G4String CCalRotationMatrixFactory::file="";
0045
0046 CCalRotationMatrixFactory* CCalRotationMatrixFactory::getInstance(const G4String & rotfile){
0047   if (rotfile=="" || rotfile==file)
0048     return getInstance();
0049   else if (file=="") {
0050     file=rotfile;
0051     return getInstance();
0052   } else {
0053     G4cerr << "ERROR: Trying to get Rotation Matrices from " << rotfile
0054          << " when previously were retrieved from " << file <<"." << G4endl;
0055     return 0;
0056   }
0057 }
0058
0059
0060 CCalRotationMatrixFactory* CCalRotationMatrixFactory::getInstance(){
0061   if (file=="") {
0062     G4cerr << "ERROR: You haven't defined which file to use for materials in "
0063          << "CCalRotationMatrixFactory::getInstance(G4String)" << G4endl;
0064     return 0;
0065   }
0066
0067   if (instance==0) {
0068     instance = new CCalRotationMatrixFactory;
0069     return instance;
0070   }
0071   else
0072     return instance;
0073 }
0074
0075 void CCalRotationMatrixFactory::setFileName(const G4String& rotfile) {
0076   if (rotfile!=file && file!="") {
0077     G4cerr << "ERROR: Trying to change Rotation Matrices file name to "
0078          << rotfile << "." << G4endl;
0079     G4cerr << "       Using previous file: " << file << G4endl;
0080   }
0081   file=rotfile;
0082 }
0083
0084 CCalRotationMatrixFactory::~CCalRotationMatrixFactory(){
0085   G4RotationMatrixTableIterator i;
0086   for(i=theMatrices.begin(); i != theMatrices.end(); ++i) {
0087     delete (*i).second;
0088   };
0089   theMatrices.clear();
0090 }
0091
0092 G4RotationMatrix* CCalRotationMatrixFactory::findMatrix(const G4String & rot) {
0093   G4RotationMatrix* retrot=0;
0094   //Rotation :NULL is no rotation so a null pointer is returned
0095   if (rot != ":NULL") {
0096     //retrot untouched if rot not found!!!
0097     G4RotationMatrixTableIterator it = theMatrices.find(rot);
0098     if (it != theMatrices.end())
0099       retrot = (*it).second;
0100   }
0101
0102   return retrot; //!!!Maybe a treatment on not-found case needed.
0103 }
0104
0105 G4RotationMatrix* CCalRotationMatrixFactory::AddMatrix(const G4String& name,
0106                                                        G4double th1,
0107                                                        G4double phi1,
0108                                                        G4double th2,
0109                                                        G4double phi2,
0110                                                        G4double th3,
0111                                                        G4double phi3){
0112   G4double sinth1, sinth2,  sinth3, costh1, costh2, costh3;
0113   G4double sinph1, sinph2, sinph3, cosph1, cosph2, cosph3;
0114   G4double TH1 = th1/deg, TH2 = th2/deg, TH3 = th3/deg;
0115   G4double PH1 = phi1/deg, PH2 = phi2/deg, PH3 = phi3/deg;
0116
0117   if (TH1 == 0.0 || TH1 == 360) {
0118     sinth1 = 0.0; costh1 = 1.0;
0119   } else if (TH1 == 90.0 || TH1 == -270) {
0120     sinth1 = 1.0; costh1 = 0.0;
0121   } else if (TH1 == 180.0 || TH1 == -180.0) {
0122     sinth1 = 0.0; costh1 = -1.0;
0123   } else if (TH1 == 270.0 || TH1 == -90.0) {
0124     sinth1 = -1.0; costh1 = 0.0;
0125   } else {
0126     sinth1 = std::sin(th1); costh1 = std::cos(th1);
0127   }
0128
0129   if (TH2 == 0.0 || TH2 == 360) {
0130     sinth2 = 0.0; costh2 = 1.0;
0131   } else if (TH2 == 90.0 || TH2 == -270) {
0132     sinth2 = 1.0; costh2 = 0.0;
0133   } else if (TH2 == 180.0 || TH2 == -180.0) {
0134     sinth2 = 0.0; costh2 = -1.0;
0135   } else if (TH2 == 270.0 || TH2 == -90.0) {
0136     sinth2 = -1.0; costh2 = 0.0;
0137   } else {
0138     sinth2 = std::sin(th2); costh2 = std::cos(th2);
0139   }
0140
0141   if (TH3 == 0.0 || TH3 == 360) {
0142     sinth3 = 0.0; costh3 = 1.0;
0143   } else if (TH3 == 90.0 || TH2 == -270) {
0144     sinth3 = 1.0; costh3 = 0.0;
0145   } else if (TH3 == 180.0 || TH3 == -180.0) {
0146     sinth3 = 0.0; costh3 = -1.0;
0147   } else if (TH3 == 270.0 || TH3 == -90.0) {
0148     sinth3 = -1.0; costh3 = 0.0;
0149   } else {
0150     sinth3 = std::sin(th3); costh3 = std::cos(th3);
0151   }
0152
0153   if (PH1 == 0.0 || PH1 == 360) {
0154     sinph1 = 0.0; cosph1 = 1.0;
0155   } else if (PH1 == 90.0 || PH1 == -270) {
0156     sinph1 = 1.0; cosph1 = 0.0;
0157   } else if (PH1 == 180.0 || PH1 == -180.0) {
0158     sinph1 = 0.0; cosph1 = -1.0;
0159   } else if (PH1 == 270.0 || PH1 == -90.0) {
0160     sinph1 = -1.0; cosph1 = 0.0;
0161   } else {
0162     sinph1 = std::sin(phi1); cosph1 = std::cos(phi1);
0163   }
0164
0165   if (PH2 == 0.0 || PH2 == 360) {
0166     sinph2 = 0.0; cosph2 = 1.0;
0167   } else if (PH2 == 90.0 || PH2 == -270) {
0168     sinph2 = 1.0; cosph2 = 0.0;
0169   } else if (PH2 == 180.0 || PH2 == -180.0) {
0170     sinph2 = 0.0; cosph2 = -1.0;
0171   } else if (PH2 == 270.0 || PH2 == -90.0) {
0172     sinph2 = -1.0; cosph2 = 0.0;
0173   } else {
0174     sinph2 = std::sin(phi2); cosph2 = std::cos(phi2);
0175   }
0176
0177   if (PH3 == 0.0 || PH3 == 360) {
0178     sinph3 = 0.0; cosph3 = 1.0;
0179   } else if (PH3 == 90.0 || PH3 == -270) {
0180     sinph3 = 1.0; cosph3 = 0.0;
0181   } else if (PH3 == 180.0 || PH3 == -180.0) {
0182     sinph3 = 0.0; cosph3 = -1.0;
0183   } else if (PH3 == 270.0 || PH3 == -90.0) {
0184     sinph3 = -1.0; cosph3 = 0.0;
0185   } else {
0186     sinph3 = std::sin(phi3); cosph3 = std::cos(phi3);
0187   }
0188
0189   //xprime axis coordinates
0190   CLHEP::Hep3Vector xprime(sinth1*cosph1,sinth1*sinph1,costh1);
0191   //yprime axis coordinates
0192   CLHEP::Hep3Vector yprime(sinth2*cosph2,sinth2*sinph2,costh2);
0193   //zprime axis coordinates
0194   CLHEP::Hep3Vector zprime(sinth3*cosph3,sinth3*sinph3,costh3);
0195
0196 #ifdef ddebug
0197   G4cout << xprime << '\t';    G4cout << yprime << '\t';    G4cout << zprime << G4endl;
0198 #endif
0199   G4RotationMatrix *rotMat = new G4RotationMatrix();
0200   rotMat->rotateAxes(xprime, yprime, zprime);
0201   if (*rotMat == G4RotationMatrix()) {
0202     // G4cerr << "WARNING: Matrix " << name << " will not be created as a rotation matrix."
0203     // G4cerr << "WARNING: Matrix " << name << " is = identity matrix. It will not be created as a rotation matrix." << G4endl;
0204     delete rotMat;
0205     rotMat=0;
0206   } else {
0207     rotMat->invert();
0208     theMatrices[name]=rotMat;
0209 #ifdef ddebug
0210     G4cout << *rotMat << G4endl;
0211 #endif
0212   }
0213
0214   return rotMat;
0215 }
0216
0217 CCalRotationMatrixFactory::CCalRotationMatrixFactory():theMatrices(){
0218
0219   G4String path = "NULL";
0220   if (std::getenv("CCAL_GLOBALPATH"))
0221     path = std::getenv("CCAL_GLOBALPATH");
0222
0223   G4cout << " ==> Opening file " << file << "..." << G4endl;
0224   std::ifstream is;
0225   G4bool ok = openGeomFile(is, path, file);
0226   if (!ok) {
0227     G4ExceptionDescription ed;
0228     ed << "Could not open file " << file << " ... Exiting!" << G4endl;
0229     G4Exception("CCalRotationMatrixFactory::CCalRotationMatrixFactory()",
0230                 "ccal002",
0231                 FatalException,ed);
0232   }
0233
0234   //////////////////////////////////////////////////
0235   // Find *DO ROTM
0236   findDO(is, G4String("ROTM"));
0237
0238   char rubish[256];
0239   G4String name;
0240
0241 #ifdef debug
0242   G4cout << "     ==> Reading Rotation Matrices... " << G4endl;
0243   G4cout << "       Name\tTheta1\tPhi1\tTheta2\tPhi2\tTheta3\tPhi3"<< G4endl;
0244 #endif
0245
0246   is >> name;
0247   while(name!="*ENDDO") {
0248     if (name.find("#.")==0) { //It is a comment.Skip line.
0249       is.getline(rubish,256,'\n');
0250     } else {
0251 #ifdef debug
0252       G4cout << "       " << name <<'\t';
0253 #endif
0254       G4double th1, phi1, th2, phi2, th3, phi3;
0255       //Get xprime axis angles
0256       is >> th1 >> phi1;
0257 #ifdef debug
0258       G4cout << th1 << '\t' << phi1 << '\t';
0259 #endif
0260       //Get yprime axis angles
0261       is >> th2 >> phi2;
0262 #ifdef debug
0263       G4cout << th2 << '\t' << phi2 << '\t';
0264 #endif
0265       //Get zprime axis angles
0266       is >> th3 >> phi3;
0267 #ifdef debug
0268       G4cout << th3 << '\t' << phi3 << '\t';
0269 #endif
0270
0271       is.getline(rubish,256,'\n');
0272 #ifdef debug
0273       G4cout << rubish << G4endl;
0274 #endif
0275
0276       AddMatrix(name, th1*deg, phi1*deg, th2*deg, phi2*deg, th3*deg, phi3*deg);
0277     }
0278
0279     is >> name;
0280   };
0281
0282   is.close();
0283   G4cout << "       "  << theMatrices.size() << " rotation matrices read in." << G4endl;
0284 }
0285
0286
0287 // 29-Jan-2004 A.R. : commented to avoid clashes with CLHEP.
0288 //                    Streaming operators for rotation matrices are
0289 //                    already defined in CLHEP::HepRotation.
0290 // std::ostream& operator<<(std::ostream& os , const G4RotationMatrix & rot){
0291 //   //  os << "( " << rot.xx() << tab << rot.xy() << tab << rot.xz() << " )" << G4endl;
0292 //   //  os << "( " << rot.yx() << tab << rot.yy() << tab << rot.yz() << " )" << G4endl;
0293 //   //  os << "( " << rot.zx() << tab << rot.zy() << tab << rot.zz() << " )" << G4endl;
0294 //
0295 //   os << "["
0296 //      << rot.thetaX()/deg << tab << rot.phiX()/deg << tab
0297 //      << rot.thetaY()/deg << tab << rot.phiY()/deg << tab
0298 //      << rot.thetaZ()/deg << tab << rot.phiZ()/deg << "]"
0299 //      << G4endl;
0300 //
0301 //   return os;
0302 // }