#include <string>
#include <vector>
#include <cassert>
#include <TFile.h>
#include <TNtupleD.h>
#include <TSystem.h>
#include <TLorentzVector.h>
#include <TVector3.h>
#include <TApplication.h>
#include <TPolyMarker3D.h>
#include "Framework/Conventions/Constants.h"
#include "Tools/Geometry/ROOTGeomAnalyzer.h"
#include "Framework/EventGen/PathLengthList.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/Utils/StringUtils.h"
#include "Framework/Utils/CmdLnArgParser.h"

Include dependency graph for gtestROOTGeometry.cxx:

Functions
void	GetCommandLineArgs (int argc, char **argv)
void	GetRandomRay (TLorentzVector &x, TLorentzVector &p)
int	GetTargetMaterial (const PathLengthList &pl)
TVector3	kDefOptRayDirection (1, 0, 0)
TVector3	kDefOptRaySurf (0, 0, 0)
int	main (int argc, char **argv)

Variables
string	gOptGeomFile
string	gOptRootGeomTopVol
TVector3	gOptRayDirection
TVector3	gOptRaySurf
double	gOptRayR
int	gOptNVtx
int	gOptTgtPdg
double	kDefOptRayR = 100

Function Documentation

◆ GetCommandLineArgs()

void GetCommandLineArgs	(	int	argc,
		char **	argv )

Definition at line 268 of file gtestROOTGeometry.cxx.

{
  CmdLnArgParser parser(argc,argv);
 
  // geometry file
  if( parser.OptionExists('f') ) {
    LOG("test", pINFO) << "Getting ROOT geometry filename";
    gOptGeomFile = parser.ArgAsString('f');
  } else {
    string base_dir = string( gSystem->Getenv("GENIE") );
    string filename = base_dir + 
          string("/data/geo/samples/BoxWithLArPbLayers.root");
    gOptGeomFile = filename;
  }
 
  // check whether an event generation volume name has been 
  // specified -- default is the 'top volume'
  if( parser.OptionExists('v') ) {
    LOG("test", pDEBUG) << "Checking for input volume name";
    gOptRootGeomTopVol = parser.ArgAsString('v');
  } else {
    LOG("test", pDEBUG) << "Using the <master volume>";
  } 
 
  // direction
  if( parser.OptionExists('d') ) {
    LOG("test", pINFO) << "Reading ray direction";
    // split the comma separated list
    vector<double> dirv = parser.ArgAsDoubleTokens('d', ",");
    assert(dirv.size() == 3);
    gOptRayDirection.SetXYZ(dirv[0],dirv[1],dirv[2]);
  } else {
    LOG("test", pINFO) << "No input ray direction - Using default";
    gOptRayDirection = kDefOptRayDirection;
  }
 
  // ray surface:
  if( parser.OptionExists('s') ) {
    LOG("test", pINFO) << "Reading ray generation surface";
    // split the comma separated list
    vector<double> rsv = parser.ArgAsDoubleTokens('s', ",");
    assert(rsv.size() == 3);
    gOptRaySurf.SetXYZ(rsv[0],rsv[1],rsv[2]);
  } else {
    LOG("test", pINFO) << "No input ray generation surface - Using default";
    gOptRaySurf = kDefOptRaySurf;
  }
 
  // ray generation area radius:
  if( parser.OptionExists('r') ) {
    LOG("test", pINFO) << "Reading radius of ray generation area";
    gOptRayR = parser.ArgAsDouble('r');
  } else {
    LOG("test", pINFO) << "No input radius of ray generation area - Using default";
    gOptRayR = kDefOptRayR;
  }
  gOptRayR = TMath::Abs(gOptRayR); // must be positive
 
  // number of vertices to generate:
  if( parser.OptionExists('n') ) {
    LOG("test", pINFO) << "Getting number of vertices to generate";
    gOptNVtx = parser.ArgAsInt('n');
  } else {
    gOptNVtx = 0;
  }
 
  // 'forced' target pdg:
  if( parser.OptionExists('p') ) {
    LOG("test", pINFO) << "Getting 'forced' target pdg";
    gOptTgtPdg = parser.ArgAsInt('p');
  } else {
    gOptTgtPdg = -1;
  }
 
  LOG("test", pNOTICE) 
    << "\n Options: "
    << "\n ROOT geometry file: " << gOptGeomFile
    << "\n ROOT geometry top volume: " << gOptRootGeomTopVol
    << "\n Ray direction: (" 
           << gOptRayDirection.X() << ", "
           << gOptRayDirection.Y() << ", "
           << gOptRayDirection.Z() << ") "
    << "\n Ray generation surface : (" 
           << gOptRaySurf.X() << ", "
           << gOptRaySurf.Y() << ", "
           << gOptRaySurf.Z() << ") "
    << "\n Ray generation area radius : " << gOptRayR 
    << "\n Number of vertices : " << gOptNVtx
    << "\n Forced targer PDG : "  << gOptTgtPdg;
 
}

References genie::CmdLnArgParser::ArgAsDouble(), genie::CmdLnArgParser::ArgAsDoubleTokens(), genie::CmdLnArgParser::ArgAsInt(), genie::CmdLnArgParser::ArgAsString(), gOptGeomFile, gOptNVtx, gOptRayDirection, gOptRayR, gOptRaySurf, gOptRootGeomTopVol, gOptTgtPdg, kDefOptRayDirection(), kDefOptRayR, kDefOptRaySurf(), LOG, genie::CmdLnArgParser::OptionExists(), pDEBUG, pINFO, and pNOTICE.

Referenced by main().

◆ GetRandomRay()

void GetRandomRay	(	TLorentzVector &	x,
		TLorentzVector &	p )

Definition at line 200 of file gtestROOTGeometry.cxx.

{
// generate a random ray (~flux neutrino)
//
  RandomGen * rnd = RandomGen::Instance();
 
  TVector3 vec0(gOptRayDirection);
  TVector3 vec = vec0.Orthogonal();
 
  double phi = 2.*kPi * rnd->RndFlux().Rndm();
  double Rt  = -1;
  bool accept = false;
  while(!accept) {
    double r = gOptRayR * rnd->RndFlux().Rndm();
    double y = gOptRayR * rnd->RndFlux().Rndm();
    if(y<r) {
      accept = true;
      Rt = r;
    }
  }
 
  vec.Rotate(phi,vec0);
  vec.SetMag(Rt);
 
  vec = vec + gOptRaySurf;
 
  TLorentzVector xx(vec, 0.);
  TLorentzVector pp(gOptRayDirection, gOptRayDirection.Mag());
 
  x = xx;
  p = pp;
 
  LOG("test", pNOTICE) 
   << "** Curr ray:";
  LOG("test", pNOTICE) 
   << "    x = " << x.X() << ",  y = " << x.Y() << ",  z = " << x.Z();
  LOG("test", pNOTICE) 
   << "    px = " << p.X() << ", py = " << p.Y() << ", pz = " << p.Z();
 
}

References gOptRayDirection, gOptRayR, gOptRaySurf, genie::RandomGen::Instance(), genie::constants::kPi, LOG, pNOTICE, and genie::RandomGen::RndFlux().

Referenced by main().

◆ GetTargetMaterial()

int GetTargetMaterial ( const PathLengthList & pl )

Definition at line 241 of file gtestROOTGeometry.cxx.

{
  if(pl.AreAllZero()) return -1;
 
  if(gOptTgtPdg > 0) {
    if(pl.PathLength(gOptTgtPdg) > 0) return gOptTgtPdg;
  } 
  else {
    RandomGen * rnd = RandomGen::Instance();
 
    PathLengthList::const_iterator pliter;
    double sum = 0;
    for(pliter = pl.begin(); pliter != pl.end(); ++pliter) {
        sum += pliter->second;
    }
    double cpl = sum * rnd->RndFlux().Rndm();
    sum = 0;
    for(pliter = pl.begin(); pliter != pl.end(); ++pliter) {
       sum += pliter->second;
       if(cpl < sum) {
          return pliter->first;
       }
    }
  }
  return -1;
}

References genie::PathLengthList::AreAllZero(), gOptTgtPdg, genie::RandomGen::Instance(), genie::PathLengthList::PathLength(), and genie::RandomGen::RndFlux().

Referenced by main().

◆ kDefOptRayDirection()

TVector3 kDefOptRayDirection	(	1	,
		0	,
		0	)

Referenced by GetCommandLineArgs().

◆ kDefOptRaySurf()

TVector3 kDefOptRaySurf	(	,
		,
		)

Referenced by GetCommandLineArgs().

◆ main()

int main	(	int	argc,
		char **	argv )

Definition at line 106 of file gtestROOTGeometry.cxx.

{
  GetCommandLineArgs(argc, argv);
 
#ifdef __GENIE_GEOM_DRIVERS_ENABLED__
 
  TApplication theApp("App", &argc, argv);
 
  // Create & configure the geometry driver
  //
  LOG("test", pINFO) 
     << "Creating a geometry driver for ROOT geometry at: " 
     << gOptGeomFile;
  geometry::ROOTGeomAnalyzer * geom_driver = 
               new geometry::ROOTGeomAnalyzer(gOptGeomFile);
  geom_driver ->SetTopVolName(gOptRootGeomTopVol);
 
  // Draw the geometry
  // & define TPolyMarker3D for drawing vertices later on
  //
  LOG("test", pINFO) 
      << "Drawing the ROOT geometry";
  geom_driver->GetGeometry()->GetTopVolume()->Draw();
 
  TPolyMarker3D * vtxp = new TPolyMarker3D();
  vtxp->SetMarkerColor(kRed);
  vtxp->SetMarkerStyle(8);
  vtxp->SetMarkerSize(0.4);
 
  // Compute & Printout the (density weighted) max path lengths
  //
  LOG("test", pINFO) 
       << "Computing max {density-weighted path lengths}";
  const PathLengthList & maxpl = geom_driver->ComputeMaxPathLengths();
 
  LOG("test", pINFO) << "Maximum math lengths: " << maxpl;
 
  TFile f("geomtest.root","recreate");
  TNtupleD vtxnt("vtxnt","","x:y:z:A:Z");
 
  TLorentzVector x(0,0,0,0);
  TLorentzVector p(0,0,0,0);
 
  int n = 0;
  while (n < gOptNVtx) {
 
    // generate a random ray 
    GetRandomRay(x,p);
 
    // compute density-weighted path lengths for each geometry
    // material for the current ray
    const PathLengthList & pl = geom_driver->ComputePathLengths(x,p);
    LOG("test",pINFO)        
       << "Current path lengths: " << pl;
 
    // select detector material (amongst all materials defined in the 
    // detector geometry -- do so based on density-weighted path lengths) 
    // or force it to the user-selected material
    int tpdg = GetTargetMaterial(pl);
    if (tpdg == -1) continue;
    LOG("test",pINFO) << "Selected target material: " << tpdg;
 
    // generate an 'interaction vertex' in the selected material
    const TVector3 & vtx = geom_driver->GenerateVertex(x,p,tpdg);
    LOG("test",pINFO) 
      << "Generated vtx: (x = " << vtx.X() 
      << ", y = " << vtx.Y() << ", z = " <<vtx.Z() << ")";
 
    // add it at the ntuple & at the vtx marker
    vtxnt.Fill(vtx.X(),vtx.Y(),vtx.Z(), 
               pdg::IonPdgCodeToA(tpdg), pdg::IonPdgCodeToZ(tpdg));
    vtxp->SetNextPoint(vtx.X(),vtx.Y(),vtx.Z()); 
 
    n++;
    LOG("test", pNOTICE) 
      << " *** Vertices generated so far: " << n;
  }
 
  // draw vertices
  vtxp->Draw("same");
  
  vtxnt.Write();
  f.Close();
  
  theApp.Run(kTRUE);
 
#else
    LOG("test", pERROR) 
       << "*** You should have enabled the geometry drivers first!";
#endif
 
  return 0;
}