00001
00002
00003
00004
00005 #include <cstring>
00006 #include <iostream>
00007
00008 #include "TrkExtAlg/Ext_track.h"
00009
00010 #include "G4ParticleTable.hh"
00011 #include "G4Navigator.hh"
00012 #include "G4VPhysicalVolume.hh"
00013 #include "G4VSolid.hh"
00014 #include "G4GeometryManager.hh"
00015 #include "G4RegionStore.hh"
00016 #include "G4ProductionCuts.hh"
00017 #include "G4Region.hh"
00018 #include "G4LogicalVolume.hh"
00019 #include "G4TransportationManager.hh"
00020 #include "G4PrimaryParticle.hh"
00021 #include "G4DynamicParticle.hh"
00022 #include "G4SDManager.hh"
00023 #include "G4SteppingManager.hh"
00024
00025
00026 #include "TrkExtAlg/ExtBesDetectorConstruction.h"
00027 #include "TrkExtAlg/ExtPhysicsList.h"
00028 #include "G4RunManagerKernel.hh"
00029
00030 using namespace std;
00031
00032
00033
00034
00035
00036 Ext_track::Ext_track() : myMsgFlag(true),myBFieldOn(true),myGeomOptimization(true),m_dir(0),m_detVer(1),myUseMucKal(true),myMucWindow(6)
00037 {
00038
00039 bes3DetectorConstruction = new ExtBesDetectorConstruction(myBFieldOn,m_detVer);
00040 bes3WorldVolume = bes3DetectorConstruction->Construct();
00041 extPhysicsList = new ExtPhysicsList;
00042 extTrack = new G4Track;
00043
00044
00045 extTrackingManager = new G4TrackingManager;
00046
00047 extRunManagerKernel = new G4RunManagerKernel;
00048 }
00049
00050 Ext_track::Ext_track(const bool msgFlag, const bool BFieldOn, const bool GeomOptimization, const int m_detVer,const bool aUseMucKal,const int aMucWindow) : myMsgFlag(msgFlag),myBFieldOn(BFieldOn),myGeomOptimization(GeomOptimization),m_dir(0),myUseMucKal(aUseMucKal),myMucWindow(aMucWindow)
00051 {
00052
00053 bes3DetectorConstruction = new ExtBesDetectorConstruction(myBFieldOn,m_detVer);
00054 bes3WorldVolume = bes3DetectorConstruction->Construct();
00055 extPhysicsList = new ExtPhysicsList;
00056 extTrack = new G4Track;
00057
00058
00059 extTrackingManager = new G4TrackingManager;
00060
00061 extRunManagerKernel = new G4RunManagerKernel;
00062 }
00063
00064
00065
00066 Ext_track::~Ext_track()
00067 {
00068 if(extRunManagerKernel) delete extRunManagerKernel;
00069 if(extTrackingManager) delete extTrackingManager;
00070
00071 if(extTrack) delete extTrack;
00072 if(bes3DetectorConstruction) delete bes3DetectorConstruction;
00073 if(extPhysicsList) delete extPhysicsList;
00074
00075
00076 G4GeometryManager::GetInstance()->OpenGeometry();
00077
00078
00079 G4SDManager* fSDM = G4SDManager::GetSDMpointerIfExist();
00080 if(fSDM)
00081 {
00082 delete fSDM;
00083 }
00084 }
00085
00086
00087
00088
00089
00090 void Ext_track::Initialization( const bool aMsgFlag,const bool Bfield,const bool GeomOptimization,const bool aUseMucKal,const int aMucWindow)
00091 {
00092 myMsgFlag=aMsgFlag;
00093 myGeomOptimization = GeomOptimization;
00094 myBFieldOn=Bfield,
00095 myUseMucKal=aUseMucKal;
00096 myMucWindow = aMucWindow;
00097
00098 G4ParticleTable::GetParticleTable()->SetReadiness();
00099 extPhysicsList->ConstructParticle();
00100
00101 if(myMsgFlag) cout << "Ext_track::Init will execute geant initialization." << endl;
00102 if(!GeometryInitialization()) cout << "Error in Ext_track::GeometryInitialization()" << endl;
00103 PhysicsInitialization();
00104
00105 extSteppingAction = new ExtSteppingAction;
00106 extSteppingAction->SetMsgFlag(aMsgFlag);
00107 extSteppingAction->SetMucKalFlag(aUseMucKal);
00108 extSteppingAction->SetMucWindow(aMucWindow);
00109
00110 extTrackingManager->SetUserAction(extSteppingAction);
00111
00112
00113
00114 }
00115
00116
00117
00119 bool Ext_track::GeometryInitialization()
00120 {
00121 if(myMsgFlag) cout << "Ext_track::GeometryInitialization()." << endl;
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131 G4Region *defaultRegion=new G4Region("DefaultRegionForBesWorld");
00132 defaultRegion->SetProductionCuts(G4ProductionCutsTable::GetProductionCutsTable()->GetDefaultProductionCuts());
00133
00134
00135 if(bes3WorldVolume->GetLogicalVolume()->GetRegion())
00136 {
00137 if(bes3WorldVolume->GetLogicalVolume()->GetRegion()!=defaultRegion)
00138 {
00139 cout << "The world volume has a user-defined region <"
00140 << bes3WorldVolume->GetLogicalVolume()->GetRegion()->GetName()
00141 << ">." << G4endl;
00142 return 0;
00143 }
00144 }
00145
00146
00147 if(defaultRegion->GetNumberOfRootVolumes())
00148 {
00149 if(defaultRegion->GetNumberOfRootVolumes()>size_t(1))
00150 {
00151 cout <<"DefaultRegionHasMoreThanOneVolume,Default world region should have a unique logical volume."<<endl;
00152 return 0;
00153 }
00154 std::vector<G4LogicalVolume*>::iterator lvItr = defaultRegion->GetRootLogicalVolumeIterator();
00155 defaultRegion->RemoveRootLogicalVolume(*lvItr);
00156 cout << (*lvItr)->GetName()
00157 << " is removed from the default region." << endl;
00158 }
00159
00160
00161 G4LogicalVolume* bes3WorldLog = bes3WorldVolume->GetLogicalVolume();
00162 bes3WorldLog->SetRegion(defaultRegion);
00163 defaultRegion->AddRootLogicalVolume(bes3WorldLog);
00164
00165
00166 G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->SetWorldVolume(bes3WorldVolume);
00167
00168 return 1;
00169 }
00170
00171
00172
00173 bool Ext_track::PhysicsInitialization()
00174 {
00175 if(myMsgFlag) cout<<"Ext_track::PhysicsInitialization()."<<endl;
00176
00177
00178
00179
00180 extPhysicsList->Construct();
00181 extPhysicsList->SetCuts();
00182
00183 CheckRegions();
00184
00185
00186 G4RegionStore::GetInstance()->UpdateMaterialList(bes3WorldVolume);
00187 G4ProductionCutsTable::GetProductionCutsTable()->UpdateCoupleTable(bes3WorldVolume);
00188
00189
00190 if(myMsgFlag) cout<<"Build PhysicsTables"<<endl;
00191 extPhysicsList->BuildPhysicsTable();
00192 if(myMsgFlag) cout<<"Build PhysicsTables end."<<endl;
00193 G4ProductionCutsTable::GetProductionCutsTable()->PhysicsTableUpdated();
00194 extPhysicsList->DumpCutValuesTableIfRequested();
00195
00196
00197 if(myGeomOptimization)
00198 {
00199 cout<<"Geometry Optimization,please wait for a few minutes."<<endl;
00200 G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
00201 geomManager->OpenGeometry();
00202 geomManager->CloseGeometry(true, false);
00203 }
00204
00205 return 1;
00206 }
00207
00208
00209
00210
00211 void Ext_track::CheckRegions()
00212 {
00213
00214 G4RegionStore::GetInstance()->SetWorldVolume();
00215
00216 for(size_t i=0;i<G4RegionStore::GetInstance()->size();i++)
00217 {
00218 G4Region* region = (*(G4RegionStore::GetInstance()))[i];
00219
00220 if(region->GetWorldPhysical()!=bes3WorldVolume) continue;
00221 G4ProductionCuts* cuts = region->GetProductionCuts();
00222 if(!cuts)
00223 {
00224 region->SetProductionCuts(G4ProductionCutsTable::GetProductionCutsTable()->GetDefaultProductionCuts());
00225 }
00226 }
00227 }
00228
00229
00230
00231
00232
00233
00234
00235 bool Ext_track::Set( const Hep3Vector &xv3, const Hep3Vector &pv3,const HepSymMatrix &err, const std::string &particleName, const double pathInMDC, const double tofInMdc)
00236 {
00237 if( err.num_row() != 6 ){
00238 std::cerr << "%ERROR at Ext_track::Set. Dimension of error matrix: "
00239 << err.num_row() << " should be 6" << std::endl;
00240 exit( 0 );
00241 }
00242
00243 m_vect[0] = xv3.x();
00244 m_vect[1] = xv3.y();
00245 m_vect[2] = xv3.z();
00246
00247
00248
00249
00250
00251
00252 if(!CheckVertexInsideWorld(xv3)) return 0;
00253
00254 float p( pv3.mag() );
00255 m_vect[3] = pv3.x()/p;
00256 m_vect[4] = pv3.y()/p;
00257 m_vect[5] = pv3.z()/p;
00258 m_vect[6] = p;
00259
00260
00261 if(particleName!="e+"&&particleName!="e-"&&
00262 particleName!="mu+"&&particleName!="mu-"&&
00263 particleName!="pi+"&&particleName!="pi-"&&
00264 particleName!="kaon+"&&particleName!="kaon-"&&
00265 particleName!="proton"&&particleName!="anti_proton"&&
00266 particleName!="gamma")
00267 {
00268 std::cerr <<"Unknown or unconstructed Particle."<< std::endl;
00269 return 0;
00270 }
00271
00272
00273 double mass;
00274 double Q;
00275
00276 G4ParticleDefinition* particleDefinition=G4ParticleTable::GetParticleTable()->FindParticle(particleName);
00277 Q = particleDefinition->GetPDGCharge();
00278 mass = particleDefinition->GetPDGMass();
00279
00280
00281 Hep3Vector xv( m_vect[0], m_vect[1], m_vect[2] );
00282 Hep3Vector pv(m_vect[3]*m_vect[6], m_vect[4]*m_vect[6], m_vect[5]*m_vect[6]);
00283
00284 m_xp_err.set_err( err, xv, pv, Q, mass );
00285
00286 extSteppingAction->SetXpErrPointer(&m_xp_err);
00287 extSteppingAction->SetInitialPath(pathInMDC);
00288 extSteppingAction->SetInitialTof(tofInMdc);
00289
00290 double betaInMDC = p/sqrt(mass*mass+p*p);
00291 extSteppingAction->SetBetaInMDC(betaInMDC);
00292
00293
00294
00295
00296 extSteppingAction->MucReset();
00297
00298
00299
00300
00301
00302
00303
00304
00305
00306
00307
00308
00309 G4DynamicParticle* DP = new G4DynamicParticle(particleDefinition,pv);
00310
00311 delete extTrack;
00312 extTrack = new G4Track(DP,0.0,xv3);
00313
00314
00315
00316 Hep3Vector center(0,0,0);
00317 G4Navigator* navigator= G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
00318 navigator->LocateGlobalPointAndSetup(center,0,false);
00319
00320 return 1;
00321 }
00322
00323
00324
00325 bool Ext_track::CheckVertexInsideWorld(const Hep3Vector& pos)
00326 {
00327 G4Navigator* navigator= G4TransportationManager::GetTransportationManager()-> GetNavigatorForTracking();
00328
00329 G4VPhysicalVolume* world= navigator-> GetWorldVolume();
00330 G4VSolid* solid= world-> GetLogicalVolume()-> GetSolid();
00331 EInside qinside= solid-> Inside(pos);
00332
00333 if( qinside != kInside) return false;
00334 else return true;
00335 }
00336
00337
00338
00339
00340
00341
00342
00343 void Ext_track::TrackExtrapotation()
00344 {
00345 extTrackingManager->ProcessOneTrack(extTrack);
00346 }
00347
00348