MaGeGeneratorPositionSampling.cc

Go to the documentation of this file.

//
// ********************************************************************
// * DISCLAIMER                                                       *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.                                                             *
// *                                                                  *
// ********************************************************************
//
// This class is supposed to sample the position of the primary
// particle randomly in a volume or on a surface. The actual
// algorithms are not implemented, but the interfaces are ok
//
// History:
// --------
// 28 Oct 2004   L.Pandola    First implementation (interfaces ok)
// 10 Nov 2004   Xiang        Implemented SampleUniformlyInVolume for Tubs 
// 11 Nov 2004   Luciano      Added check if the sampled point is inside (now it should work also if 
//                            the Tub has an asymmetric hole, as Ge detectors)
// 19 Nov 2004   Luciano      A lot of changes: implemented a modified 
//                            version of Dave Jordan algorithm
//
#include "generators/MaGeGeneratorPositionSampling.hh"
#include "generators/MJGeneratorUtil.hh"
#include "io/MJLogger.hh"
#include "Randomize.hh"
#include "G4LogicalVolume.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4SolidStore.hh"
#include "G4Tubs.hh"
#include "G4Box.hh"
#include "G4VisExtent.hh"
#include "G4Navigator.hh"
#include "G4TransportationManager.hh"

//Da mettere dentro:
// intera catena di inclusioni
// numero di copia

MaGeGeneratorPositionSampling::MaGeGeneratorPositionSampling() :
  myVolume(0)
{
  volumeName="NULL";
  myNavigator = G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
  generatorUtil = new MJGeneratorUtil();
  copyNumber = 0;
  fSolid_type = "NULL";
  radius = 0;

}

MaGeGeneratorPositionSampling::~MaGeGeneratorPositionSampling()
{
  delete generatorUtil;
}

G4ThreeVector MaGeGeneratorPositionSampling::SampleUniformlyInVolume(G4String volName,G4int copyN)
{
  MJLog(debugging) << "Starting SampleUniformlyInVolume" << endlog;
  G4ThreeVector vertex_world(0.,0.,0.); //vertex in world coordinates
  
  if (volName == "NULL") 
     {
       MJLog(warning) << 
        "The confinement is set but the volume was not specified" << endlog;
       MJLog(warning) << "Return default position" << endlog; 
       return vertex_world;
     }

  //Step 1: check if the volume changed since the last call
  G4bool ifound = true;
  if (volName != volumeName || copyN != copyNumber) //volume changed
    {  
      volumeName = volName;
      copyNumber = copyN;
      ifound = InitializeSamplingVolume();
      if (!ifound)
        {
          MJLog(warning) << 
            "The confinement is set but the volume was not specified" << endlog;
          MJLog(warning) << "Return default position" << endlog; 
          return vertex_world;
        }
      else 
        {
          MJLog(debugging) << "Volume initialized" << endlog;
        }
    }

  //Step 2: sample the point in the three cases
  G4ThreeVector vertex;       // vertex in local volume's coordinates
 
  //If the type is G4Box or G4Tubs it's not necessary 
  //to shoot and check many times
  if (fSolid_type == "G4Box") 
    {
      G4double x_lo = fSolid_par_arr[0];
      G4double x_hi = fSolid_par_arr[1];
      G4double y_lo = fSolid_par_arr[2];
      G4double y_hi = fSolid_par_arr[3];
      G4double z_lo = fSolid_par_arr[4];
      G4double z_hi = fSolid_par_arr[5];
      vertex = generatorUtil->pick_point_in_box(x_lo, x_hi, y_lo, y_hi, z_lo, z_hi);
      vertex_world = (GetGlobalRotation()*vertex) + GetGlobalTranslation();
    }
  if (fSolid_type == "G4Tubs")
    {
      G4double r1      = fSolid_par_arr[0];
      G4double r2      = fSolid_par_arr[1];
      G4double h       = fSolid_par_arr[2] * 2.0;
      G4double theta0  = fSolid_par_arr[3];
      G4double d_theta = fSolid_par_arr[4];
      vertex = generatorUtil->pick_point_in_annulus(r1, r2, h, theta0, d_theta);
      vertex_world = (GetGlobalRotation()*vertex) + GetGlobalTranslation();
    }

  //for other volumes and for non-located boxes/tubs one must instead shoot and try
  if (fSolid_type == "BoundingSphere")
    {
      G4bool point_in_volume = false;
      G4int iloop=0;
      while (!point_in_volume && iloop<10000) //prevents infinite loops
        {  
          // loop until point falls in physical vol.
          vertex = generatorUtil->pick_point_in_sphere(radius);
          //It is a sphere: no need to rotate!
          vertex_world = vertex + GetGlobalTranslation();
          G4VPhysicalVolume* pVol_sampled = 
            myNavigator->LocateGlobalPointAndSetup(vertex_world);
          if (pVol_sampled == myVolume) {
            point_in_volume = true;
          }
          iloop++;
        }
    }

  MJLog(debugging) << "Sampled vertex " << vertex_world.x() << " " << 
    vertex_world.y() << " " << vertex_world.z() << " mm" << endlog;
  MJLog(debugging) << "Volume: " <<  
    myNavigator->LocateGlobalPointAndSetup(vertex_world)->GetName() << endlog;

  if (myNavigator->LocateGlobalPointAndSetup(vertex_world) != myVolume) 
    MJLog(warning) << "The sampled point is not in the right volume!?!" << endlog;

  return vertex_world;
}

G4bool MaGeGeneratorPositionSampling::InitializeSamplingVolume()
{
  MJLog(debugging) << "Starting InitializeSamplingVolume" << endlog;  
  MJLog(debugging) << "Volume to be found: " << volumeName << endlog; 

  //For now it works if there are no copies
  G4PhysicalVolumeStore* volumeStore = G4PhysicalVolumeStore::GetInstance();
  G4int nVolumes = (G4int) volumeStore->size();
  G4bool ifound = false;
  G4int nfound=0;
  G4int i;
  for (i=0; i<nVolumes; i++) {
    if ((*volumeStore)[i]->GetName() == volumeName) {
      ifound = true;
      nfound++;
      myVolume = (*volumeStore)[i];   
    }
  }

  if (nfound > 1) {
    ifound = false;
    MJLog(warning) << "The volume name is not unique: retrieve copy number" << endlog;
    for (i=0; i<nVolumes; i++) {
      if ((*volumeStore)[i]->GetName() == volumeName && 
          (*volumeStore)[i]->GetCopyNo() == copyNumber && !ifound){
        ifound = true;
        myVolume = (*volumeStore)[i];
      }
    }
  }

  if (!ifound) return ifound;

  MJLog(debugging) << "Found volume: " << myVolume->GetName() << ", copy " << myVolume->GetCopyNo() << 
    endlog;
 
  //Store the dimensions of solids for Tubs and Boxes
  G4VSolid* pSolid = myVolume->GetLogicalVolume()->GetSolid();
  fSolid_type = pSolid->GetEntityType();

  MJLog(debugging) << "pSolid->GetEntityType()= " << fSolid_type << endlog;

  if (fSolid_type == "G4Box") 
    {
      MJLog(debugging) << "Identified source solid type as G4Box." << endlog;
      MJLog(debugging) << "Geometry parameters for box shape follow:" << endlog;
      G4double x_half_length = ((G4Box*)pSolid)->GetXHalfLength();
      G4double y_half_length = ((G4Box*)pSolid)->GetYHalfLength();
      G4double z_half_length = ((G4Box*)pSolid)->GetZHalfLength();
      
      MJLog(debugging) << "X half-length (mm): " << x_half_length << endlog;
      MJLog(debugging) << "Y half-length (mm): " << y_half_length << endlog;
      MJLog(debugging) << "Z half-length (mm): " << z_half_length << endlog;
      
      fSolid_par_arr[0] = -x_half_length;
      fSolid_par_arr[1] =  x_half_length;
      fSolid_par_arr[2] = -y_half_length;
      fSolid_par_arr[3] =  y_half_length;
      fSolid_par_arr[4] = -z_half_length;
      fSolid_par_arr[5] =  z_half_length;

    }
    else if (fSolid_type == "G4Tubs") 
      {
        MJLog(debugging) << "Identified source solid type as G4Tubs." << endlog;
        MJLog(debugging) << "Geometry parameters for annulus shape follow:" << endlog;
        fSolid_par_arr[0] = ((G4Tubs*)pSolid)->GetInnerRadius();
        fSolid_par_arr[1] = ((G4Tubs*)pSolid)->GetOuterRadius();
        fSolid_par_arr[2] = ((G4Tubs*)pSolid)->GetZHalfLength();
        fSolid_par_arr[3] = ((G4Tubs*)pSolid)->GetStartPhiAngle();
        fSolid_par_arr[4] = ((G4Tubs*)pSolid)->GetDeltaPhiAngle();
        MJLog(debugging) << "Inner radius (mm):        " << fSolid_par_arr[0] << endlog;
        MJLog(debugging) << "Outer radius (mm):        " << fSolid_par_arr[1] << endlog;
        MJLog(debugging) << "Z half-length (mm):       " << fSolid_par_arr[2] << endlog;
        MJLog(debugging) << "Starting phi angle (rad): " << fSolid_par_arr[3] << endlog;
        MJLog(debugging) << "Delta phi angle (rad):    " << fSolid_par_arr[4] << endlog;
      }
  else {
    MJLog(debugging) << "The source volume solid is neither G4Box nor G4Tubs." << endlog;
    fSolid_type = "BoundingSphere";
  }

  //MJLog(debugging) << "Retrieving bounding sphere dimensions...." << endlog;
  radius = pSolid->GetExtent().GetExtentRadius();
  MJLog(debugging) << "Sphere for solid. " << "Radius " << radius/mm << " mm" << endlog;

  G4VPhysicalVolume* worldV = myNavigator->GetWorldVolume();
  if (!worldV) MJLog(error) << "World volume not defined" << endlog;
  
  G4VPhysicalVolume* testVolume = myVolume;
  G4ThreeVector globTrasl(0,0,0);
  G4RotationMatrix  rotM;   // Initialised to identity
  G4int iloop = 0;
  while (testVolume != worldV && iloop < nVolumes) 
    {
      MJLog(debugging) << "Going up: volume " << testVolume->GetName() << endlog;
      globTrasl = globTrasl + testVolume->GetObjectTranslation();
      rotM = (*(testVolume->GetObjectRotation())) * rotM;
      testVolume = FindTheDirectMother(testVolume);
      iloop++; 
    }
  SetGlobalRotation(rotM);
  SetGlobalTranslation(globTrasl);
  
  //Just dump infos
  MJLog(debugging) << "Global volume translation: " << 
    globTrasl.x() << " " << globTrasl.y() << " " << globTrasl.z() << " mm" << endlog;
  MJLog(debugging) << "Global volume rotation: " << endlog;
  for (i=0;i<3;i++) {
    MJLog(debugging) << rotM(i,0) << " " << rotM(i,1) << " " << rotM(i,2) << 
      endlog;
  }
  return ifound;
}         

G4ThreeVector MaGeGeneratorPositionSampling::SampleOnSurface(G4String volName,G4int copyN)
{
  
  MJLog(debugging) << "Starting SampleOnSurface" << endlog;
  G4ThreeVector vertex_world(0.,0.,0.); //vertex in world coordinates
  
  if (volName == "NULL") 
     {
       MJLog(warning) << 
        "The confinement is set but the volume was not specified" << endlog;
       MJLog(warning) << "Return default position" << endlog; 
       return vertex_world;
     }

  //Step 1: check if the volume changed since the last call
  G4bool ifound = true;
  if (volName != volumeName || copyN != copyNumber) //volume changed
    {  
      volumeName = volName;
      copyNumber = copyN;
      ifound = InitializeSamplingVolume();
      if (!ifound)
        {
          MJLog(warning) << 
            "The confinement is set but the volume was not specified" << endlog;
          MJLog(warning) << "Return default position" << endlog; 
          return vertex_world;
        }
      else 
        {
          MJLog(debugging) << "Volume initialized" << endlog;
        }
    }

  //Step 2: sample the point in the three cases
  G4ThreeVector vertex;       // vertex in local volume's coordinates
 
  return vertex_world;
}

G4VPhysicalVolume* MaGeGeneratorPositionSampling::FindTheDirectMother(G4VPhysicalVolume* vol)
{
  //Find the direct mother of vol
  G4PhysicalVolumeStore* volumeStore = G4PhysicalVolumeStore::GetInstance();
  G4int nVolumes = (G4int) volumeStore->size();
  std::vector<G4VPhysicalVolume*> auxVect; 

  G4int i,j;
  //1st step: look for all the ancestors
  for (i=0; i<nVolumes; i++) {
    if ((*volumeStore)[i]->GetLogicalVolume()->IsAncestor(vol)) 
      {
        //MJLog(debugging) << "Ancestor found: " << (*volumeStore)[i]->GetName() << endlog;
        auxVect.push_back((*volumeStore)[i]);
      }
  }
  G4int nAncestors = (G4int) auxVect.size();
  std::vector<G4VPhysicalVolume*>::iterator iter = auxVect.begin();
  //MJLog(debugging) << "Initial number of ancestors: " << nAncestors << endlog;
  
  //2nd step: remove the non-direct ancestors
  G4bool erased = false;
  while (nAncestors > 1) {
    erased = false;
    //MJLog(debugging) << "Number of ancestors: " << nAncestors << endlog;
    for (i=0;i<nAncestors;i++) 
      {
        for (j=0;j<nAncestors;j++) 
          {
            iter = auxVect.begin();
            G4VPhysicalVolume* vol1 = auxVect[i];
            G4VPhysicalVolume* vol2 = auxVect[j];
            if (vol1->GetLogicalVolume()->IsAncestor(vol2)) 
              //vol1 is ancestor of vol2
              //so vol1 is not a direct ancestor of the target volume!
              {
                iter += i; //moves to the i-th element of the vector
                auxVect.erase(iter); //removes vol1 from the list
                erased = true;
              }
            if (erased) break;
          }
        if (erased) break;
      }
    nAncestors = (G4int) auxVect.size(); //updates the dimension of the list
    //Here ends the loop
  }
  return auxVect[0];
}

---