g4sim: src/MyFieldSvc.cc Source File

00001 //---------------------------------------------------------------------------//
00002 //Description: Setup MagField
00003 //Author: Wu Chen(wuchen@mail.ihep.ac.cn)
00004 //Created: 17 Oct, 2012
00005 //Comment:
00006 //---------------------------------------------------------------------------//
00007 
00008 #include "MyFieldSvc.hh"
00009 
00010 #include "G4UnitsTable.hh"
00011 #include "myglobals.hh"
00012 
00013 #include "G4UniformElectricField.hh"
00014 #include "G4UniformMagField.hh"
00015 #include "G4MagneticField.hh"
00016 #include "G4FieldManager.hh"
00017 #include "G4TransportationManager.hh"
00018 #include "G4EquationOfMotion.hh"
00019 #include "G4EqMagElectricField.hh"
00020 #include "G4Mag_UsualEqRhs.hh"
00021 #include "G4MagIntegratorStepper.hh"
00022 #include "G4MagIntegratorDriver.hh"
00023 #include "G4ChordFinder.hh"
00024 
00025 #include "G4ExplicitEuler.hh"
00026 #include "G4ImplicitEuler.hh"
00027 #include "G4SimpleRunge.hh"
00028 #include "G4SimpleHeum.hh"
00029 #include "G4ClassicalRK4.hh"
00030 #include "G4HelixExplicitEuler.hh"
00031 #include "G4HelixImplicitEuler.hh"
00032 #include "G4HelixSimpleRunge.hh"
00033 #include "G4CashKarpRKF45.hh"
00034 #include "G4RKG3_Stepper.hh"
00035 
00036 #include <sstream>
00037 #include <fstream>
00038 #include <string.h>
00039 #include <iostream>
00040 
00041 #include "MyGlobalField.hh"
00042 #include "MyFieldMap.hh"
00043 
00044 #include "MyFieldSvcMessenger.hh"
00045 
00046 MyFieldSvc* MyFieldSvc::fMyFieldSvc = 0;
00047 
00048         MyFieldSvc::MyFieldSvc()
00049 : fChordFinder(0), fStepper(0), fIntgrDriver(0), fMagField(0), fEleField(0)
00050 {
00051         if (fMyFieldSvc){
00052                 G4Exception("MyFieldSvc::MyFieldSvc()","Run0031",
00053                                 FatalException, "MyFieldSvc constructed twice.");
00054         }
00055         fMyFieldSvc = this;
00056         fFieldManager = GetGlobalFieldManager();
00057 
00058         //Magnetic Field
00059     (void)MyGlobalField::getObject();
00060 
00061         //Electric Field
00062         fEquation = new G4EqMagElectricField(fEleField); 
00063 
00064         fMyFieldSvcMessenger= new MyFieldSvcMessenger(this);
00065 }
00066 
00067 MyFieldSvc::~MyFieldSvc()
00068 {
00069         if(fMagField) delete fMagField;
00070         delete fMyFieldSvcMessenger;
00071 }
00072 
00073 MyFieldSvc* MyFieldSvc::GetMyFieldSvc(){
00074         if ( !fMyFieldSvc ){
00075                 fMyFieldSvc = new MyFieldSvc;
00076         }
00077         return fMyFieldSvc;
00078 }
00079 
00080 void MyFieldSvc::SetField(G4LogicalVolume* fLogicWorld){
00081         Dump();
00082         //Setup the field
00083         G4String opt = "";
00084         if(fMagField) delete fMagField;
00085         if(fEleField) delete fEleField;
00086         if(fFieldMaps.size()>0){
00087                 for (int ifm = 0; ifm < fFieldMaps.size(); ifm++){
00088                         if (fFieldMaps[ifm]) delete fFieldMaps[ifm];
00089                 }
00090                 fFieldMaps.clear();
00091         }
00092 
00093         if ( fType == "Uniform" ){
00094                 G4ThreeVector MagField_vec(1,0,0);
00095                 MagField_vec.setTheta(UniF_Theta);
00096                 MagField_vec.setPhi(UniF_Phi);
00097                 MagField_vec = MagField_vec.unit() * UniF_Intensity;
00098                 if(MagField_vec!= G4ThreeVector(0.,0.,0.))
00099                 { 
00100                         fMagField = new G4UniformMagField(MagField_vec);
00101                         opt = "Mag";
00102                 }
00103                 else{
00104                         fMagField = 0;
00105                         opt = "Mag_0";
00106                 }
00107         }
00108         else if ( fType == "Electric_Uniform" ){
00109                 G4ThreeVector EleField_vec(1,0,0);
00110                 EleField_vec.setTheta(UniEF_Theta);
00111                 EleField_vec.setPhi(UniEF_Phi);
00112                 EleField_vec = EleField_vec.unit() * UniEF_Intensity;
00113                 if(EleField_vec!= G4ThreeVector(0.,0.,0.))
00114                 { 
00115                         fEleField = new  G4UniformElectricField(EleField_vec);
00116                         fEquation->SetFieldObj(fEleField);  // must now point to the new field
00117                         opt = "Ele";
00118                 }
00119                 else 
00120                 {
00121                         // If the new field's value is Zero, then it is best to
00122                         //  insure that it is not used for propagation.
00123                         fEleField = 0;
00124                         fEquation->SetFieldObj(fEleField);   // As a double check ...
00125                         fEleField = 0;
00126                         opt = "Ele_0";
00127                 }
00128         }
00129         else if ( fType == "fieldMap" ){
00130                 for(G4int ifieldmap = 0; ifieldmap < fFieldMapFilenames.size(); ifieldmap++){
00131                         fFieldMaps.push_back(new MyFieldMap(fFieldMapFilenames[ifieldmap],
00132                                   fFieldMapScalings[ifieldmap], // "current" in G4Beamline field parlance
00133                                   1.0, // gradient
00134                                   0.0, // time offset
00135                                   fLogicWorld));
00136                                 G4cout << "Adding MyFieldMap: "<< fFieldMapFilenames[ifieldmap]<<", "<<
00137                                 fFieldMapScalings[ifieldmap]<<G4endl;
00138                 }
00139                 opt = "fieldMap";
00140         }
00141         else if ( fType != "None" && fType != "none" && fType != "NONE"){
00142                 std::cout<<"Field Type "<<fType<<" is not not supported yet!!!"<<std::endl;
00143                 G4Exception("MyFieldSvc::SetField()","Run0031",
00144                                 FatalException, "Field Type is not not supported.");
00145         }
00146         UpdateField(opt);
00147 }
00148 
00150 //
00151 // Register this field to 'global' Field Manager and 
00152 // Create Stepper and Chord Finder with predefined type, minstep (resp.)
00153 //
00154 
00155 void MyFieldSvc::UpdateField(G4String opt)
00156 {
00157         if ( opt == "Mag_0" || opt == "Mag" ){
00158                 fFieldManager->SetDetectorField(fMagField);
00159                 if (opt == "Mag"){
00160                         fFieldManager->CreateChordFinder(fMagField);
00161                 }
00162         } 
00163         else if ( opt == "Ele" || opt == "Ele_0" ){
00164                 SetStepper();
00165                 fFieldManager->SetDetectorField(fEleField );
00166                 if ( opt == "Ele" ){
00167                         if(fChordFinder) delete fChordFinder;
00168                         // fChordFinder = new G4ChordFinder( fEleField, UniEF_StepL, fStepper);
00169                         fIntgrDriver = new G4MagInt_Driver(UniEF_StepL,fStepper,fStepper->GetNumberOfVariables() );
00170                         fChordFinder = new G4ChordFinder(fIntgrDriver);
00171                         fFieldManager->SetChordFinder( fChordFinder );
00172                 }
00173         }
00174         else if ( opt == "fieldMap" ){
00175                 // Nothing else should be done. this part is handled by RunAction instead
00176         }
00177 }
00178 
00180 //
00181 // Set stepper according to the stepper type
00182 //
00183 
00184 void MyFieldSvc::SetStepper()
00185 {
00186         G4int nvar = 8;
00187 
00188         if(fStepper) delete fStepper;
00189 
00190         switch ( UniEF_StepT ) 
00191         {
00192                 case 0:  
00193                         fStepper = new G4ExplicitEuler( fEquation, nvar ); 
00194                         G4cout<<"G4ExplicitEuler is calledS"<<G4endl;     
00195                         break;
00196                 case 1:  
00197                         fStepper = new G4ImplicitEuler( fEquation, nvar );      
00198                         G4cout<<"G4ImplicitEuler is called"<<G4endl;     
00199                         break;
00200                 case 2:  
00201                         fStepper = new G4SimpleRunge( fEquation, nvar );        
00202                         G4cout<<"G4SimpleRunge is called"<<G4endl;     
00203                         break;
00204                 case 3:  
00205                         fStepper = new G4SimpleHeum( fEquation, nvar );         
00206                         G4cout<<"G4SimpleHeum is called"<<G4endl;     
00207                         break;
00208                 case 4:  
00209                         fStepper = new G4ClassicalRK4( fEquation, nvar );       
00210                         G4cout<<"G4ClassicalRK4 (default) is called"<<G4endl;     
00211                         break;
00212                 case 5:  
00213                         fStepper = new G4CashKarpRKF45( fEquation, nvar );      
00214                         G4cout<<"G4CashKarpRKF45 is called"<<G4endl;     
00215                         break;
00216                 case 6:  
00217                         fStepper = 0; // new G4RKG3_Stepper( fEquation, nvar );       
00218                         G4cout<<"G4RKG3_Stepper is not currently working for Electric Field"<<G4endl;     
00219                         break;
00220                 case 7:  
00221                         fStepper = 0; // new G4HelixExplicitEuler( fEquation ); 
00222                         G4cout<<"G4HelixExplicitEuler is not valid for Electric Field"<<G4endl;     
00223                         break;
00224                 case 8:  
00225                         fStepper = 0; // new G4HelixImplicitEuler( fEquation ); 
00226                         G4cout<<"G4HelixImplicitEuler is not valid for Electric Field"<<G4endl;     
00227                         break;
00228                 case 9:  
00229                         fStepper = 0; // new G4HelixSimpleRunge( fEquation );   
00230                         G4cout<<"G4HelixSimpleRunge is not valid for Electric Field"<<G4endl;     
00231                         break;
00232                 default: fStepper = 0;
00233         }
00234 }
00235 
00237 //
00238 //  Utility method
00239 
00240 G4FieldManager*  MyFieldSvc::GetGlobalFieldManager()
00241 {
00242         return G4TransportationManager::GetTransportationManager()
00243                 ->GetFieldManager();
00244 }
00245 
00247 //
00248 // Interface
00249 //
00250 
00251 void MyFieldSvc::ReadCard( G4String file_name ){
00252         if(file_name[0] != '/'){ // Relative Dir
00253                 G4String dir_name = getenv("CONFIGUREROOT");
00254                 if (dir_name[dir_name.size()-1] != '/') dir_name.append("/");
00255                 file_name = dir_name + file_name;
00256         }
00257         std::ifstream fin_card(file_name);
00258         if ( !fin_card ){
00259                 std::cout<<"MagField card"<<file_name<<" is not available!!!"<<std::endl;
00260                 G4Exception("MyFieldSvc::SetField()","Run0031",
00261                                 FatalException, "MagField card is not available.");
00262         }
00263         std::stringstream buf_card;
00264         std::string s_card;
00265         while(getline(fin_card,s_card)){
00266                 buf_card.str("");
00267                 buf_card.clear();
00268                 buf_card<<s_card;
00269                 const char* c = s_card.c_str();
00270                 int length = strlen(c);
00271                 int offset = 0;
00272                 for ( ; offset < length; offset++ ){
00273                         if ( c[offset] != ' ' ) break;
00274                 }
00275                 if ( c[offset] == '#' ) continue;
00276                 else if ( c[offset] == '/' && c[offset+1] == '/' ) continue;
00277                 else if ( length - offset == 0 ) continue;
00278                 std::string keyword;
00279                 buf_card>>keyword;
00280                 if ( keyword == "TYPE:" ){
00281                         buf_card>>fType;
00282                         continue;
00283                 }
00284                 else if ( keyword == "UniF:" ){
00285                         buf_card>>UniF_Intensity>>UniF_Theta>>UniF_Phi;
00286                         UniF_Intensity *= tesla;
00287                         UniF_Theta *= deg;
00288                         UniF_Phi *= deg;
00289                 }
00290                 else if ( keyword == "UniEF:" ){
00291                         buf_card>>UniEF_Intensity>>UniEF_Theta>>UniEF_Phi>>UniEF_StepL>>UniEF_StepT;
00292                         UniEF_Intensity *= kilovolt/cm;
00293                         UniEF_Theta *= deg;
00294                         UniEF_Phi *= deg;
00295                         UniEF_StepL *= mm;
00296                 }
00297                 else if ( keyword == "fMap:" ){
00298                         G4String name;
00299                         G4double scale;
00300                         buf_card>>name>>scale;
00301                         G4String dir_name = getenv("FIELDMAPSROOT");
00302                         if (dir_name[dir_name.size()-1] != '/') dir_name.append("/");
00303                         name = dir_name + name;
00304                         fFieldMapFilenames.push_back(name);
00305                         fFieldMapScalings.push_back(scale);
00306                 }
00307                 else{
00308                         std::cout<<"In MyFieldSvc::ReadCard, unknown name: "<<keyword<<" in file "<<file_name<<std::endl;
00309                         std::cout<<"Will ignore this line!"<<std::endl;
00310                 }
00311         }
00312         fin_card.close();
00313         buf_card.str("");
00314         buf_card.clear();
00315 }
00316 
00317 void MyFieldSvc::Dump(){
00318         std::cout<<"**********************MagField Settings***************************"<<std::endl;
00319         std::cout<<"Type: "<<fType<<std::endl;
00320         if ( fType == "Uniform" ){
00321                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10) <<"Intensity"
00322                         <<std::setiosflags(std::ios::left)<<std::setw(10) <<"Direction(Theta, Phi)"
00323                         <<std::endl;
00324                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10)<<"T"
00325                         <<std::setiosflags(std::ios::left)<<std::setw(5) <<"deg"
00326                         <<std::setiosflags(std::ios::left)<<std::setw(5) <<"deg"
00327                         <<std::endl;
00328                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10)<<UniF_Intensity/tesla
00329                         <<std::setiosflags(std::ios::left)<<std::setw(5) <<UniF_Theta/deg
00330                         <<std::setiosflags(std::ios::left)<<std::setw(5) <<UniF_Phi/deg
00331                         <<std::endl;
00332         }
00333         else if ( fType == "Electric_Uniform" ){
00334                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10) <<"Intensity"
00335                         <<std::setiosflags(std::ios::left)<<std::setw(20) <<"Dir(Theta,Phi)"
00336                         <<std::setiosflags(std::ios::left)<<std::setw(6)  <<"StepL"
00337                         <<std::setiosflags(std::ios::left)<<std::setw(6)  <<"StepT"
00338                         <<std::endl;
00339                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10)<<"kV/cm"
00340                         <<std::setiosflags(std::ios::left)<<std::setw(10)<<"deg"
00341                         <<std::setiosflags(std::ios::left)<<std::setw(10)<<"deg"
00342                         <<std::setiosflags(std::ios::left)<<std::setw(6) <<"mm"
00343                         <<std::endl;
00344                 std::cout<<std::setiosflags(std::ios::left)<<std::setw(10)<<UniEF_Intensity/kilovolt*cm
00345                         <<std::setiosflags(std::ios::left)<<std::setw(10)<<UniEF_Theta/deg
00346                         <<std::setiosflags(std::ios::left)<<std::setw(10)<<UniEF_Phi/deg
00347                         <<std::setiosflags(std::ios::left)<<std::setw(6) <<UniEF_StepL/mm
00348                         <<std::setiosflags(std::ios::left)<<std::setw(6) <<UniEF_StepT
00349                         <<std::endl;
00350         }
00351         else if ( fType == "fieldMap" ){
00352                 for(G4int ifieldmap = 0; ifieldmap < fFieldMapFilenames.size(); ifieldmap++){
00353                         std::cout<<std::setiosflags(std::ios::left)<<std::setw(30) <<fFieldMapFilenames[ifieldmap]<<" "
00354                                 <<std::setiosflags(std::ios::left)<<std::setw(6) <<fFieldMapScalings[ifieldmap]
00355                                 <<std::endl;
00356                 }
00357         }
00358         else if ( fType == "none" || fType == "NONE" || fType == "None" ){
00359                 std::cout<<"No Field!"<<std::endl;
00360         }
00361         std::cout<<"******************************************************************"<<std::endl;
00362 }