GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gtestROOTGeometry.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*!
3 
4 \program gtestROOTGeometry
5 
6 \brief Tests the GENIE ROOT geometry driver by generating random rays,
7  following them through the detector, computing density-weighted
8  path-lengths for each material & generating vertices
9 
10 \syntax gtestROOTGeometry [-f geom] [-n nvtx] [-d dx,dy,dz] [-s x,y,z]
11  [-r size] [-p pdg]
12 
13  Options:
14 
15  -f A ROOT file containing a ROOT/GEANT geometry description
16  [default: $GENIE/data/geo/samples/BoxWithLArPbLayers.root]
17  -n Number of vertices to generate
18  -d Direction of generated rays (dirx,diry,dirz)
19  [default: 1,0,0 - along x]
20  -s Specifies the ray generation surface.
21  That surface is perpendicular to the beam direction and
22  contains the point specified here.
23  Use the same units as your input geometry.
24  [default: 0,0,0]
25  -r Specifies the radius of the area over which rays will be
26  generated. That event generation area is the area contained
27  within the circle centered at the point specified with the
28  -s option and a radius specified here.
29  Use the same units as your input geometry.
30  [default: 100 in your geometry unit]
31  -p Can be used to specify a target pdg code. If that option is set
32  then the vertex is generated only at that material. If not set
33  then vertices will be generated in all detector materials (by
34  weight).
35  -v Can specify specific volumes of the input geometry. If not set
36  will use the master volume.
37 
38  Examples:
39 
40  gtestROOTGeometry -n 20000 -d 1,0,0 -s -10,0,0 -r 10 -p 1000180400
41 
42  will take a test geom ($GENIE/data/geo/samples/BoxWithLArPbLayers.root)
43  and generate 20k vertices in Ar40 (pdg=1000180400) only, using rays
44  (flux) generated uniformly (in area) within a circle of radius = 10
45  (in units of that input geometry, m in this case), centered at
46  (-10,0,0) (see prev comment on units) and perpendicular to the ray
47  direction (1,0,0).
48 
49 \Author Anselmo Meregaglia <anselmo.meregaglia@cern.ch>
50  ETH Zurich
51 
52  Costas Andreopoulos <c.andreopoulos \at cern.ch>
53  STFC, Rutherford Appleton Lab
54 
55 \created August 11, 2005
56 
57 \cpright Copyright (c) 2003-2024, The GENIE Collaboration
58  For the full text of the license visit http://copyright.genie-mc.org
59 
60 */
61 //____________________________________________________________________________
62 
63 #include <string>
64 #include <vector>
65 #include <cassert>
66 
67 #include <TFile.h>
68 #include <TNtupleD.h>
69 #include <TSystem.h>
70 #include <TLorentzVector.h>
71 #include <TVector3.h>
72 #include <TApplication.h>
73 #include <TPolyMarker3D.h>
74 
82 
83 using std::string;
84 using std::vector;
85 
86 using namespace genie;
87 using namespace genie::constants;
88 
89 void GetCommandLineArgs (int argc, char ** argv);
90 void GetRandomRay (TLorentzVector & x, TLorentzVector & p);
91 int GetTargetMaterial (const PathLengthList & pl);
92 
93 string gOptGeomFile; // input ROOT geom file
94 string gOptRootGeomTopVol; // top volume / can be used to override the master volume
95 TVector3 gOptRayDirection; // ray direction
96 TVector3 gOptRaySurf; // ray generation surface
97 double gOptRayR; // ray generation area radius
98 int gOptNVtx; // number of vertices to generate
99 int gOptTgtPdg; //
100 
101 double kDefOptRayR = 100;
102 TVector3 kDefOptRayDirection (1,0,0);
103 TVector3 kDefOptRaySurf (0,0,0);
104 
105 //____________________________________________________________________________
106 int main(int argc, char ** argv)
107 {
108  GetCommandLineArgs(argc, argv);
109 
110 #ifdef __GENIE_GEOM_DRIVERS_ENABLED__
111 
112  TApplication theApp("App", &argc, argv);
113 
114  // Create & configure the geometry driver
115  //
116  LOG("test", pINFO)
117  << "Creating a geometry driver for ROOT geometry at: "
118  << gOptGeomFile;
119  geometry::ROOTGeomAnalyzer * geom_driver =
121  geom_driver ->SetTopVolName(gOptRootGeomTopVol);
122 
123  // Draw the geometry
124  // & define TPolyMarker3D for drawing vertices later on
125  //
126  LOG("test", pINFO)
127  << "Drawing the ROOT geometry";
128  geom_driver->GetGeometry()->GetTopVolume()->Draw();
129 
130  TPolyMarker3D * vtxp = new TPolyMarker3D();
131  vtxp->SetMarkerColor(kRed);
132  vtxp->SetMarkerStyle(8);
133  vtxp->SetMarkerSize(0.4);
134 
135  // Compute & Printout the (density weighted) max path lengths
136  //
137  LOG("test", pINFO)
138  << "Computing max {density-weighted path lengths}";
139  const PathLengthList & maxpl = geom_driver->ComputeMaxPathLengths();
140 
141  LOG("test", pINFO) << "Maximum math lengths: " << maxpl;
142 
143  TFile f("geomtest.root","recreate");
144  TNtupleD vtxnt("vtxnt","","x:y:z:A:Z");
145 
146  TLorentzVector x(0,0,0,0);
147  TLorentzVector p(0,0,0,0);
148 
149  int n = 0;
150  while (n < gOptNVtx) {
151 
152  // generate a random ray
153  GetRandomRay(x,p);
154 
155  // compute density-weighted path lengths for each geometry
156  // material for the current ray
157  const PathLengthList & pl = geom_driver->ComputePathLengths(x,p);
158  LOG("test",pINFO)
159  << "Current path lengths: " << pl;
160 
161  // select detector material (amongst all materials defined in the
162  // detector geometry -- do so based on density-weighted path lengths)
163  // or force it to the user-selected material
164  int tpdg = GetTargetMaterial(pl);
165  if (tpdg == -1) continue;
166  LOG("test",pINFO) << "Selected target material: " << tpdg;
167 
168  // generate an 'interaction vertex' in the selected material
169  const TVector3 & vtx = geom_driver->GenerateVertex(x,p,tpdg);
170  LOG("test",pINFO)
171  << "Generated vtx: (x = " << vtx.X()
172  << ", y = " << vtx.Y() << ", z = " <<vtx.Z() << ")";
173 
174  // add it at the ntuple & at the vtx marker
175  vtxnt.Fill(vtx.X(),vtx.Y(),vtx.Z(),
177  vtxp->SetNextPoint(vtx.X(),vtx.Y(),vtx.Z());
178 
179  n++;
180  LOG("test", pNOTICE)
181  << " *** Vertices generated so far: " << n;
182  }
183 
184  // draw vertices
185  vtxp->Draw("same");
186 
187  vtxnt.Write();
188  f.Close();
189 
190  theApp.Run(kTRUE);
191 
192 #else
193  LOG("test", pERROR)
194  << "*** You should have enabled the geometry drivers first!";
195 #endif
196 
197  return 0;
198 }
199 //____________________________________________________________________________
200 void GetRandomRay(TLorentzVector & x, TLorentzVector & p)
201 {
202 // generate a random ray (~flux neutrino)
203 //
204  RandomGen * rnd = RandomGen::Instance();
205 
206  TVector3 vec0(gOptRayDirection);
207  TVector3 vec = vec0.Orthogonal();
208 
209  double phi = 2.*kPi * rnd->RndFlux().Rndm();
210  double Rt = -1;
211  bool accept = false;
212  while(!accept) {
213  double r = gOptRayR * rnd->RndFlux().Rndm();
214  double y = gOptRayR * rnd->RndFlux().Rndm();
215  if(y<r) {
216  accept = true;
217  Rt = r;
218  }
219  }
220 
221  vec.Rotate(phi,vec0);
222  vec.SetMag(Rt);
223 
224  vec = vec + gOptRaySurf;
225 
226  TLorentzVector xx(vec, 0.);
227  TLorentzVector pp(gOptRayDirection, gOptRayDirection.Mag());
228 
229  x = xx;
230  p = pp;
231 
232  LOG("test", pNOTICE)
233  << "** Curr ray:";
234  LOG("test", pNOTICE)
235  << " x = " << x.X() << ", y = " << x.Y() << ", z = " << x.Z();
236  LOG("test", pNOTICE)
237  << " px = " << p.X() << ", py = " << p.Y() << ", pz = " << p.Z();
238 
239 }
240 //____________________________________________________________________________
242 {
243  if(pl.AreAllZero()) return -1;
244 
245  if(gOptTgtPdg > 0) {
246  if(pl.PathLength(gOptTgtPdg) > 0) return gOptTgtPdg;
247  }
248  else {
249  RandomGen * rnd = RandomGen::Instance();
250 
251  PathLengthList::const_iterator pliter;
252  double sum = 0;
253  for(pliter = pl.begin(); pliter != pl.end(); ++pliter) {
254  sum += pliter->second;
255  }
256  double cpl = sum * rnd->RndFlux().Rndm();
257  sum = 0;
258  for(pliter = pl.begin(); pliter != pl.end(); ++pliter) {
259  sum += pliter->second;
260  if(cpl < sum) {
261  return pliter->first;
262  }
263  }
264  }
265  return -1;
266 }
267 //____________________________________________________________________________
268 void GetCommandLineArgs(int argc, char ** argv)
269 {
270  CmdLnArgParser parser(argc,argv);
271 
272  // geometry file
273  if( parser.OptionExists('f') ) {
274  LOG("test", pINFO) << "Getting ROOT geometry filename";
275  gOptGeomFile = parser.ArgAsString('f');
276  } else {
277  string base_dir = string( gSystem->Getenv("GENIE") );
278  string filename = base_dir +
279  string("/data/geo/samples/BoxWithLArPbLayers.root");
280  gOptGeomFile = filename;
281  }
282 
283  // check whether an event generation volume name has been
284  // specified -- default is the 'top volume'
285  if( parser.OptionExists('v') ) {
286  LOG("test", pDEBUG) << "Checking for input volume name";
287  gOptRootGeomTopVol = parser.ArgAsString('v');
288  } else {
289  LOG("test", pDEBUG) << "Using the <master volume>";
290  }
291 
292  // direction
293  if( parser.OptionExists('d') ) {
294  LOG("test", pINFO) << "Reading ray direction";
295  // split the comma separated list
296  vector<double> dirv = parser.ArgAsDoubleTokens('d', ",");
297  assert(dirv.size() == 3);
298  gOptRayDirection.SetXYZ(dirv[0],dirv[1],dirv[2]);
299  } else {
300  LOG("test", pINFO) << "No input ray direction - Using default";
302  }
303 
304  // ray surface:
305  if( parser.OptionExists('s') ) {
306  LOG("test", pINFO) << "Reading ray generation surface";
307  // split the comma separated list
308  vector<double> rsv = parser.ArgAsDoubleTokens('s', ",");
309  assert(rsv.size() == 3);
310  gOptRaySurf.SetXYZ(rsv[0],rsv[1],rsv[2]);
311  } else {
312  LOG("test", pINFO) << "No input ray generation surface - Using default";
314  }
315 
316  // ray generation area radius:
317  if( parser.OptionExists('r') ) {
318  LOG("test", pINFO) << "Reading radius of ray generation area";
319  gOptRayR = parser.ArgAsDouble('r');
320  } else {
321  LOG("test", pINFO) << "No input radius of ray generation area - Using default";
323  }
324  gOptRayR = TMath::Abs(gOptRayR); // must be positive
325 
326  // number of vertices to generate:
327  if( parser.OptionExists('n') ) {
328  LOG("test", pINFO) << "Getting number of vertices to generate";
329  gOptNVtx = parser.ArgAsInt('n');
330  } else {
331  gOptNVtx = 0;
332  }
333 
334  // 'forced' target pdg:
335  if( parser.OptionExists('p') ) {
336  LOG("test", pINFO) << "Getting 'forced' target pdg";
337  gOptTgtPdg = parser.ArgAsInt('p');
338  } else {
339  gOptTgtPdg = -1;
340  }
341 
342  LOG("test", pNOTICE)
343  << "\n Options: "
344  << "\n ROOT geometry file: " << gOptGeomFile
345  << "\n ROOT geometry top volume: " << gOptRootGeomTopVol
346  << "\n Ray direction: ("
347  << gOptRayDirection.X() << ", "
348  << gOptRayDirection.Y() << ", "
349  << gOptRayDirection.Z() << ") "
350  << "\n Ray generation surface : ("
351  << gOptRaySurf.X() << ", "
352  << gOptRaySurf.Y() << ", "
353  << gOptRaySurf.Z() << ") "
354  << "\n Ray generation area radius : " << gOptRayR
355  << "\n Number of vertices : " << gOptNVtx
356  << "\n Forced targer PDG : " << gOptTgtPdg;
357 
358 }
359 //____________________________________________________________________________
double PathLength(int pdgc) const
bool AreAllZero(void) const
virtual const PathLengthList & ComputeMaxPathLengths(void)
#define pERROR
Definition: Messenger.h:59
string gOptGeomFile
static RandomGen * Instance()
Access instance.
Definition: RandomGen.cxx:71
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
A singleton holding random number generator classes. All random number generation in GENIE should tak...
Definition: RandomGen.h:29
int main(int argc, char **argv)
Definition: gAtmoEvGen.cxx:327
double gOptRayR
virtual const PathLengthList & ComputePathLengths(const TLorentzVector &x, const TLorentzVector &p)
TVector3 gOptRayDirection
Object to be filled with the neutrino path-length, for all detector geometry materials, when starting from a position x and travelling along the direction of the neutrino 4-momentum.
TVector3 gOptRaySurf
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
string gOptRootGeomTopVol
Definition: gAtmoEvGen.cxx:301
double kDefOptRayR
A ROOT/GEANT4 geometry driver.
TVector3 kDefOptRaySurf(0, 0, 0)
void GetRandomRay(TLorentzVector &x, TLorentzVector &p)
int gOptNVtx
#define pINFO
Definition: Messenger.h:62
virtual const TVector3 & GenerateVertex(const TLorentzVector &x, const TLorentzVector &p, int tgtpdg)
virtual void SetTopVolName(string nm)
int GetTargetMaterial(const PathLengthList &pl)
TVector3 kDefOptRayDirection(1, 0, 0)
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
Command line argument parser.
#define pNOTICE
Definition: Messenger.h:61
void GetCommandLineArgs(int argc, char **argv)
Definition: gAtmoEvGen.cxx:563
TRandom3 & RndFlux(void) const
rnd number generator used by flux drivers
Definition: RandomGen.h:71
int gOptTgtPdg
virtual TGeoManager * GetGeometry(void) const
#define pDEBUG
Definition: Messenger.h:63