gtestROOTGeometry.cxx File Reference

#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"

Go to the source code of this file.

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

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;

}

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();

}

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;
}

TVector3 kDefOptRayDirection	(	1	,
		0	,
		0
	)

TVector3 kDefOptRaySurf	(	0	,
		0	,
		0
	)

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;
}

Variable Documentation

string gOptGeomFile

Definition at line 93 of file gtestROOTGeometry.cxx.

int gOptNVtx

Definition at line 98 of file gtestROOTGeometry.cxx.

TVector3 gOptRayDirection

Definition at line 95 of file gtestROOTGeometry.cxx.

double gOptRayR

Definition at line 97 of file gtestROOTGeometry.cxx.

TVector3 gOptRaySurf

Definition at line 96 of file gtestROOTGeometry.cxx.

string gOptRootGeomTopVol

Definition at line 94 of file gtestROOTGeometry.cxx.

int gOptTgtPdg

Definition at line 99 of file gtestROOTGeometry.cxx.

double kDefOptRayR = 100

Definition at line 101 of file gtestROOTGeometry.cxx.