/search.css" rel="stylesheet" type="text/css"/> /search.js">
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

In This Package:

GtRockGammaTool.cc
Go to the documentation of this file.
00001 #include "GtRockGammaTool.h"
00002 
00003 #include "GaudiKernel/IRndmGenSvc.h"
00004 #include "GaudiKernel/IParticlePropertySvc.h"
00005 #include "GaudiKernel/ParticleProperty.h"
00006 #include "GaudiKernel/Vector3DTypes.h"
00007 #include "GaudiKernel/PhysicalConstants.h"
00008 #include "HepMC/GenEvent.h"
00009 #include "HepMC/GenParticle.h"
00010 #include "HepMC/GenVertex.h"
00011 
00012 #include "CLHEP/Units/SystemOfUnits.h"
00013 #include "CLHEP/Units/PhysicalConstants.h"
00014 
00015 #include "DetHelpers/IPositionerTool.h"
00016 
00017 using namespace CLHEP;
00018 
00019 GtSourceRect::GtSourceRect(const CLHEP::Hep3Vector& origin,
00020                            const CLHEP::Hep3Vector& edgeA,
00021                            const CLHEP::Hep3Vector& edgeB)
00022   : m_origin(origin),
00023     m_edgeA(edgeA),
00024     m_edgeB(edgeB)
00025 {
00026   // Pre-calculate area and normal vector
00027   m_area = edgeA.mag()*edgeB.mag();
00028   m_normal = edgeA.cross(edgeB).unit();
00029 }
00030 
00031 GtSourceRect::~GtSourceRect()
00032 {
00033   ;
00034 }
00035 
00036 double GtSourceRect::area()
00037 {
00038   return m_area;
00039 }
00040 
00041 CLHEP::Hep3Vector GtSourceRect::origin()
00042 {
00043   return m_origin;
00044 }
00045 
00046 CLHEP::Hep3Vector GtSourceRect::normal()
00047 {
00048   return m_normal;
00049 }
00050 
00051 CLHEP::Hep3Vector GtSourceRect::randomPoint(Rndm::Numbers& random)
00052 {
00053   // Choose a random point on the rectangular surface
00054   double x=m_origin.x();
00055   double y=m_origin.y();
00056   double z=m_origin.z();
00057   double randA = random()-0.5;
00058   double randB = random()-0.5;
00059   x+=randA*m_edgeA.x();
00060   y+=randA*m_edgeA.y();
00061   z+=randA*m_edgeA.z();
00062   x+=randB*m_edgeB.x();
00063   y+=randB*m_edgeB.y();
00064   z+=randB*m_edgeB.z();
00065   return CLHEP::Hep3Vector(x,y,z);
00066 }
00067 
00068 GtSourceOct::GtSourceOct(const CLHEP::Hep3Vector& origin,
00069                          const CLHEP::Hep3Vector& edgeA,
00070                          const CLHEP::Hep3Vector& edgeB,
00071                          const double& bevelLength)
00072   : GtSourceRect(origin,edgeA,edgeB),
00073     m_bevelLength(bevelLength)
00074 {
00075   // Pre-calculate area
00076   m_area = edgeA.mag()*edgeB.mag() - 2*bevelLength*bevelLength;
00077 }
00078 
00079 GtSourceOct::~GtSourceOct()
00080 {
00081   ;
00082 }
00083 
00084 CLHEP::Hep3Vector GtSourceOct::randomPoint(Rndm::Numbers& random)
00085 {
00086   // Choose a random point on the octagonal surface
00087   static const double sqrt2 = sqrt(2);
00088   double cutDist = ((m_edgeA.mag()+m_edgeB.mag())/(2*sqrt2)) 
00089     - m_bevelLength/sqrt2;
00090   CLHEP::Hep3Vector point;
00091   while(true){
00092     point = GtSourceRect::randomPoint(random);
00093     // Check if point is in beveled corners
00094     double pointDist = (fabs(point.x()) + fabs(point.y())) / sqrt2;
00095     if( pointDist < cutDist ) break;
00096   }
00097   return point;
00098 }
00099 
00100 
00101 GtRockGammaTool::GtRockGammaTool(const std::string& type,
00102                                  const std::string& name, 
00103                                  const IInterface* parent) 
00104   :GaudiTool(type,name,parent),
00105    m_totalSurfaceArea(0),
00106    m_walls()
00107 {
00108   // Initialization
00109   m_pdfGammas.clear();   
00110   m_pdfGammas.push_back(1);
00111   m_edgesGammas.clear(); 
00112   m_edgesGammas.push_back(0.5*MeV); 
00113   m_edgesGammas.push_back(10*MeV);
00114 
00115   declareInterface<IHepMCEventMutator>(this);
00116   declareProperty("GammaSpectrum",m_pdfGammas,"User-Defined PDF for Gamma Spectrum");
00117   declareProperty("GammaBinEdges",m_edgesGammas,"Bin Edges of Gamma Spectrum");
00118 
00119 }
00120 
00121 GtRockGammaTool::~GtRockGammaTool() 
00122 {
00123   for(unsigned int wallIdx=0;wallIdx<m_walls.size();wallIdx++){
00124     delete m_walls[wallIdx];
00125     m_walls[wallIdx] = 0;
00126   }
00127   m_walls.clear();
00128 }
00129 
00130 StatusCode GtRockGammaTool::initialize() 
00131 {
00132   // Initialize the rock gamma tool
00133   // Random number service
00134   IRndmGenSvc *rgs = 0;
00135   if (service("RndmGenSvc",rgs,true).isFailure()) {
00136     fatal() << "Failed to get random service" << endreq;
00137     return StatusCode::FAILURE;        
00138   }
00139   
00140   StatusCode sc;
00141   sc = m_uni.initialize(rgs, Rndm::Flat(0,1));
00142   if (sc.isFailure()) {
00143     fatal() << "Failed to initialize uniform random numbers" << endreq;
00144     return StatusCode::FAILURE;
00145   }
00146 
00147   // check user input for user-defined PDFs
00148   if (m_pdfGammas.size()+1 != m_edgesGammas.size() ) {
00149     fatal() << "There should be 1 more number of bin edges than bins when defining pdf" << endreq;
00150     return StatusCode::FAILURE; 
00151   }
00152 
00153   // Fill pre-calculated Cumulative Distribution Functions
00154   m_cdfGammas.push_back( 0 );
00155   for (unsigned int bin=0; bin<m_pdfGammas.size(); bin++){
00156     m_cdfGammas.push_back( m_cdfGammas[bin] + m_pdfGammas[bin] );
00157   }
00158   // Normalize
00159   for (unsigned int bin=0; bin<m_cdfGammas.size(); bin++){
00160     m_cdfGammas[bin] /= m_cdfGammas.back();
00161   }
00162 
00163   // Pre-calculate some useful quantities 
00164 
00165   // FIXME: Current surfaces are defined as an octagon consistent with
00166   // the water pool dead volume, but the top surface of the octagon
00167   // has been shifted vertically to encompass most of the RPC.  It
00168   // would be better to define two volumes (Pool+RPC) to completely
00169   // encompass the RPC.
00170 
00171   // This defines 10 walls for the octagonal space:
00172   //        +Y
00173   // -X+Y /----\ +X+Y
00174   //  -X |      | +X    and (+Z, -Z)
00175   // -X-Y \----/ +X-Y
00176   //        -Y
00177 
00178   // FIXME: should extract these numbers from DetDesc Geometry
00179   double poolDeadSizeZ=10.0*m;
00180   double nearPoolDeadSizeX=16.0*m;
00181   double nearPoolDeadSizeY=10.0*m;
00182   //double farPoolDeadSizeX=16.0*m;
00183   //double farPoolDeadSizeY=16.0*m;
00184   // Bevel size is the distance that the octagonal corner walls remove
00185   // from the initial rectangular pool shape.  It is twice the
00186   // distance from the bevel wall to the corner that would be formed
00187   // by the X,Y walls.
00188   double poolDeadBevelSize=4.249*m;
00189   // RPC space above detector
00190   // Note: I have pre-rotated the RPC XY coords by 90deg to match the pool
00191   double rpcRoofSpace=0.50*m;
00192   //double rpcNearRoofSizeX=18.5*m;
00193   //double rpcNearRoofSizeY=12.5*m;
00194   double rpcNearRoofSizeZ=0.30*m;
00195   double rpcNearRoofHeightZ= rpcRoofSpace + rpcNearRoofSizeZ;
00196   //double rpcFarRoofSizeX=18.5*m;
00197   //double rpcFarRoofSizeY=18.5*m;
00198   //double rpcFarRoofSizeZ=0.30*m;
00199   //double rpcFarRoofHeightZ= rpcRoofSpace + rpcFarRoofSizeZ;
00200 
00201   // Calculate some helper parameters
00202   // FIXME: Hardcoded near side walls.  Must change for Far site.
00203   double poolAndRpcZ = poolDeadSizeZ + rpcNearRoofHeightZ;
00204   double sqrt2 = sqrt(2);
00205   // Extra RPC space shifts wall center higher by half the RPC thickness
00206   double zOrigin = rpcNearRoofHeightZ / 2.;
00207   // The bevelXYOffset is the distance along the X or Y axis removed by the bevel.
00208   double bevelXYOffset = poolDeadBevelSize/sqrt2;
00209 
00210   m_walls.clear();
00211   {
00212     // +X wall
00213     CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.,0,zOrigin);
00214     CLHEP::Hep3Vector edgeA(0,
00215                             0,
00216                             poolAndRpcZ);
00217     CLHEP::Hep3Vector edgeB(0,
00218                             (nearPoolDeadSizeY-2*bevelXYOffset),
00219                             0);
00220     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00221   }
00222   {
00223     // +X+Y wall
00224     CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.-bevelXYOffset/2.,
00225                              nearPoolDeadSizeY/2.-bevelXYOffset/2.,
00226                              zOrigin);
00227     CLHEP::Hep3Vector edgeA(0,
00228                             0,
00229                             poolAndRpcZ);
00230     CLHEP::Hep3Vector edgeB(-bevelXYOffset,
00231                             bevelXYOffset,
00232                             0);
00233     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00234   }
00235   {
00236     // +Y wall
00237     CLHEP::Hep3Vector origin(0,nearPoolDeadSizeY/2.,zOrigin);
00238     CLHEP::Hep3Vector edgeA(0,
00239                             0,
00240                             poolAndRpcZ);
00241     CLHEP::Hep3Vector edgeB(-(nearPoolDeadSizeX-2*bevelXYOffset),
00242                             0,
00243                             0);
00244     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00245   }
00246   {
00247     // -X+Y wall
00248     CLHEP::Hep3Vector origin(-(nearPoolDeadSizeX/2.-bevelXYOffset/2.),
00249                              nearPoolDeadSizeY/2.-bevelXYOffset/2.,
00250                              zOrigin);
00251     CLHEP::Hep3Vector edgeA(0,
00252                             0,
00253                             poolAndRpcZ);
00254     CLHEP::Hep3Vector edgeB(-bevelXYOffset,
00255                             -bevelXYOffset,
00256                             0);
00257     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00258   }
00259   {
00260     // -X wall
00261     CLHEP::Hep3Vector origin(-nearPoolDeadSizeX/2.,0,zOrigin);
00262     CLHEP::Hep3Vector edgeA(0,
00263                             0,
00264                             poolAndRpcZ);
00265     CLHEP::Hep3Vector edgeB(0,
00266                             -(nearPoolDeadSizeY-2*bevelXYOffset),
00267                             0);
00268     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00269   }
00270   {
00271     // -X-Y wall
00272     CLHEP::Hep3Vector origin(-(nearPoolDeadSizeX/2.-bevelXYOffset/2.),
00273                              -(nearPoolDeadSizeY/2.-bevelXYOffset/2.),
00274                              zOrigin);
00275     CLHEP::Hep3Vector edgeA(0,
00276                             0,
00277                             poolAndRpcZ);
00278     CLHEP::Hep3Vector edgeB(bevelXYOffset,
00279                             -bevelXYOffset,
00280                             0);
00281     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00282   }
00283   {
00284     // -Y wall
00285     CLHEP::Hep3Vector origin(0,-nearPoolDeadSizeY/2.,zOrigin);
00286     CLHEP::Hep3Vector edgeA(0,
00287                             0,
00288                             poolAndRpcZ);
00289     CLHEP::Hep3Vector edgeB(nearPoolDeadSizeX-2*bevelXYOffset,
00290                             0,
00291                             0);
00292     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00293   }
00294   {
00295     // +X-Y wall
00296     CLHEP::Hep3Vector origin(nearPoolDeadSizeX/2.-bevelXYOffset/2.,
00297                              -(nearPoolDeadSizeY/2.-bevelXYOffset/2.),
00298                              zOrigin);
00299     CLHEP::Hep3Vector edgeA(0,
00300                             0,
00301                             poolAndRpcZ);
00302     CLHEP::Hep3Vector edgeB(bevelXYOffset,
00303                             bevelXYOffset,
00304                             0);
00305     m_walls.push_back( new GtSourceRect(origin,edgeA,edgeB) );
00306   }
00307   {
00308     // +Z wall
00309     CLHEP::Hep3Vector origin(0,0,poolDeadSizeZ/2. + rpcNearRoofHeightZ);
00310     CLHEP::Hep3Vector edgeA(0,
00311                             nearPoolDeadSizeY,
00312                             0);
00313     CLHEP::Hep3Vector edgeB(nearPoolDeadSizeX,
00314                             0,
00315                             0);
00316     m_walls.push_back( new GtSourceOct(origin,edgeA,edgeB,bevelXYOffset) );
00317   }
00318   {
00319     // -Z wall
00320     CLHEP::Hep3Vector origin(0,0,-poolDeadSizeZ/2.);
00321     CLHEP::Hep3Vector edgeA(nearPoolDeadSizeX,
00322                             0,
00323                             0);
00324     CLHEP::Hep3Vector edgeB(0,
00325                             nearPoolDeadSizeY,
00326                             0);
00327     m_walls.push_back( new GtSourceOct(origin,edgeA,edgeB,bevelXYOffset) );
00328   }
00329 
00330   m_totalSurfaceArea = 0;
00331   for(unsigned int wallIdx=0;wallIdx<m_walls.size();wallIdx++){
00332     debug() << "Adding wall with area: " << m_walls[wallIdx]->area() << endreq;
00333     m_totalSurfaceArea += m_walls[wallIdx]->area();
00334   }
00335   debug() << "Total surface area: " << m_totalSurfaceArea << endreq;
00336 
00337   return this->GaudiTool::initialize();
00338 }
00339 
00340 StatusCode GtRockGammaTool::finalize()
00341 {
00342   return this->GaudiTool::finalize();
00343 }
00344 
00345 StatusCode GtRockGammaTool::mutate(HepMC::GenEvent& event)
00346 { 
00347   if (this->oneVertex(event).isFailure()) 
00348     return StatusCode::FAILURE;
00349   return StatusCode::SUCCESS;
00350 }
00351 
00352 StatusCode GtRockGammaTool::oneVertex(HepMC::GenEvent& event)
00353 {
00354   // generate vertex time (within single event)
00355   double vertex_t = 0;
00356 
00357   // generate vertex position in local pool coordinates
00358   // (random point around pool edges, top, or bottom)
00359   unsigned int wallIdx = GetRandomWall();
00360   // Choose a random point on the wall
00361   CLHEP::Hep3Vector point = m_walls[wallIdx]->randomPoint(m_uni);
00362   HepMC::GenVertex* vertex = new HepMC::GenVertex(HepMC::FourVector(point.x(),
00363                                                                     point.y(),
00364                                                                     point.z(),
00365                                                                     vertex_t));
00366   if( !event.add_vertex(vertex) ){
00367     error() << "Failed to add vertex to event" << endreq;
00368     return StatusCode::FAILURE;
00369   }
00370 
00371   // generate momentum direction
00372   CLHEP::Hep3Vector dir;
00373   double dir_costh = m_uni(); // cos(theta) range 0 to 1
00374   double dir_phi = m_uni()*twopi; // range from 0 to 2*pi radians
00375   
00376   dir.setRThetaPhi(1.0,acos(dir_costh),dir_phi);
00377 
00378   // rotate along Y-axis by vertex direction's theta
00379   dir.rotateY(m_walls[wallIdx]->normal().theta());
00380   // rotate along Z-axis by vertex direction's phi
00381   dir.rotateZ(m_walls[wallIdx]->normal().phi());
00382   // Normalize direction vector
00383   dir = dir.unit();
00384 
00385   // generate random energy
00386   double energy = ConvertCdfRand(m_uni(),m_cdfGammas,m_edgesGammas); 
00387   HepMC::FourVector fourmom(energy*dir.x(),
00388                             energy*dir.y(),
00389                             energy*dir.z(),
00390                             energy);
00391   
00392   HepMC::GenParticle* particle = new HepMC::GenParticle(fourmom,22,
00393                                                         1/*=status*/);
00394   // Set photon polarization
00395   double pol_phi = m_uni()*CLHEP::twopi;
00396   particle->set_polarization(HepMC::Polarization(0.5*pi,pol_phi));
00397   // Add gamma to vertex
00398   vertex->add_particle_out(particle);
00399 
00400   return StatusCode::SUCCESS;
00401 }
00402 
00403 double GtRockGammaTool::ConvertCdfRand(const double& rand, 
00404                                        const std::vector<double>& cdf, 
00405                                        const std::vector<double>& edges) {
00406   // Convert a uniform random number in [0,1] to a random sample from
00407   // the given cumulative distribution function.
00408   int left = 0;
00409   int right = cdf.size();
00410   const double* cdfStart = &cdf[0];
00411   const double* edgesStart = &edges[0];
00412 
00413   while( right-left > 1 ){
00414     int mid = (right+left)/2;
00415     if( rand < *(cdfStart+mid) ) right = mid;
00416     else left = mid;
00417   }
00418   double value = *(edgesStart + left);
00419   double slope = (*(edgesStart+right) - *(edgesStart+left))
00420                 /(*(cdfStart+right) - *(cdfStart+left));
00421   value += slope * (rand - *(cdfStart+left));
00422   return value;
00423 }
00424 
00425 unsigned int GtRockGammaTool::GetRandomWall() {
00426   // Choose a random wall, based on surface area
00427   double rand = m_uni();
00428   double randArea = rand*m_totalSurfaceArea;
00429   unsigned int wallIdx = 0;
00430   while(wallIdx<m_walls.size()){
00431     if( randArea < m_walls[wallIdx]->area() ){
00432       // Found the correct wall
00433       break;
00434     }
00435     randArea -= m_walls[wallIdx]->area();
00436     wallIdx++;
00437   }
00438   return wallIdx;
00439 }
| Classes | Job Modules | Data Objects | Services | Algorithms | Tools | Packages | Directories | Tracs |

Generated on Fri May 16 2014 10:18:50 for GenTools by doxygen 1.7.4