// PosGen_FillShellRopes.cc // Contact person: Juan-Pablo Yanez // See PosGen_FillShellRopes.hh for more details //———————————————————————// #include using namespace CLHEP; #include #include #include #include #include #include #include #include #include #include #include namespace RAT { PosGen_FillShellRopes::PosGen_FillShellRopes(const char *arg_dbname) : GLG4PosGen(arg_dbname), pos(0.,0.,0.), max_iterations(10000000) { } // PosGen_FillShellRopes::SetState parses a string of generator parameters // to set the appropriate local variables. void PosGen_FillShellRopes::SetState(G4String newValues) { std::vector params = split(newValues, " "); // it is a GLG4 convention that SetState with a null string argument // should print usage information if (newValues.length() == 0) { G4cout << "Current state of this GLG4PosGen_PointPaintFill:\n" << " \"" << GetState() << "\"\n" << G4endl; G4cout << "Usage: x_mm y_mm z_mm r_in r_out volname" << G4endl; return; } // Setting the values periodicity = (360./10.)*pi/180.; // 36deg in rad mincosphi = 0.6; maxcostheta = 0.72; double x, y, z; G4ThreeVector vpos; x = to_double(params[0]); y = to_double(params[1]); z = to_double(params[2]); ri = to_double(params[3]); ro = to_double(params[4]); volname = params[5]; // hold_up or hold_down pos = G4ThreeVector(x, y, z); if(ri > ro) { info << "PosGen_FillShellRopes::SetState: Inner radius is greater than outer radius, flipping them.\n"; double radius_temp = ro; ro = ri; ri = radius_temp; } } // PosGen_FillShellRopes::GetState returns a G4String containing the values // of all generator parameters, formatted like the input string G4String PosGen_FillShellRopes::GetState() const { std::string rv = dformat("%d %d %d %d %d %s", pos.x(), pos.y(), pos.z(), ri, ro,volname.c_str()); return G4String(rv); } // PosGen_FillShellRopes::GeneratePosition sets argResult to a point uniformly // sampled from a spherical shell centered at the given origin (pos) and // within the specified physical volume (volname). void PosGen_FillShellRopes::GeneratePosition(G4ThreeVector &argResult) { G4Navigator* gNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking(); G4ThreeVector rpos; G4bool holdup_system; int iterations = 0; while(1) { // a point is found by sampling uniformly in r^3, cos(theta), and phi double r = pow(G4UniformRand() * (pow(ro, 3.0) - pow(ri, 3.0)) + pow(ri, 3.0), 1.0/3); double phi = pi * G4UniformRand(); double theta = acos(2.0 * G4UniformRand() - 1); if(G4UniformRand() > 0.5) theta = -theta; double x = r * sin(theta) * cos(phi); double y = r * sin(theta) * sin(phi); double z = r * cos(theta); double phi2 = atan2(y,x); double costheta = z/r; rpos = G4ThreeVector(x, y, z); rpos += pos; iterations++; Log::Assert(iterations mincosphi) && (costhetaLocateGlobalPointAndSetup(rpos,0,true)->GetName(); //G4cout<<"Volume generated "<