test_fill_dk2nu.C
Go to the documentation of this file.
1 /// \file test_fill_dk2nu.C
2 /// \brief Test creating and filling a TTree based on:
3 /// dk2nu.h (dk2nu.C) - decays of beam particles to a neutrino
4 /// dkmeta.h (dkmeta.C) - metadata for the file
5 ///
6 /// also show the use of reading location information, generating
7 /// weights at those locations, and how to extend the tree (for one-off
8 /// tests) without modifying the standard class.
9 ///
10 /// This script can be run using:
11 /// root -b -q test_fill_dk2nu.C+
12 ///
13 /// \author Robert Hatcher <rhatcher \at fnal.gov>
14 /// Fermi National Accelerator Laboratory
15 ///
16 /// \created 2012-04-03
17 /// \modified 2012-10-03
18 /// \version $Id: test_fill_dk2nu.C,v 1.3 2012-11-15 21:56:53 rhatcher Exp $
19 ///==========================================================================
20 
21 #include <iostream>
22 #include <iomanip>
23 
24 #include "TFile.h"
25 #include "TTree.h"
26 #include "TRandom3.h"
27 
28 #ifndef __CINT__
29 // hide header stuff from CINT, assume load_dk2nu.C run first
30 
31 #include "tree/dk2nu.h"
32 #include "tree/dkmeta.h"
33 
34 /// include standardized code for reading location text file
35 #include "tree/readWeightLocations.h"
36 
37 /// include standardized code for getting energy/weight vectors for locations
38 #include "tree/calcLocationWeights.h"
39 
40 #endif // ifndef __CINT__
41 
42 // if running bare CINT don't try adding the NonStd addition it won't work
43 #define ADD_NONSTD
44 
45 #ifdef ADD_NONSTD
46 /// example class for extending the tree with non-standard extras
47 /// that doesn't require modifying the basic "dk2nu" class
48 class NonStd {
49  public:
50  NonStd() { }
51  virtual ~NonStd() { }
52  void Clear() { }
53  double foo; ///< data for my one-off test
54  double bar; ///< more data
55  ClassDef(NonStd,1)
56 };
57 
58 /// make a dictionary for classes used in the tree
59 /// again do this because we have no external linkages to libraries
60 #pragma link C++ class NonStd+;
61 #endif
62 
63 // flugg 500K POT lowth files seem to have 510000 as an upper limit on
64 // # of entries. So to test for estimate of file size one needs to have
65 // that many entries _and_ semi-sensible values for all branches (so
66 // compression isn't better than it would be in real life).
67 void test_fill_dk2nu(unsigned int nentries=1000)
68 {
69 
70  // stuff...
71  TRandom3* rndm = new TRandom3();
72 
73  ///-----------------------------------------------------------------------
74  ///
75  /// equivalent to NumiAnalysis::NumiAnalysis() in g4numi
76  ///
77  ///-----------------------------------------------------------------------
78 
79  // create objects
82 #ifdef ADD_NONSTD
83  NonStd* nonstd = new NonStd;
84 #endif
85 
86  // read the text file for locations, fill the dkmeta object
87  std::string locfilename = "$DK2NU/etc/locations.txt";
88  bsim::readWeightLocations(locfilename,dkmeta);
89 
90  // print out what we have for locations
91  size_t nloc = dkmeta->location.size();
92  std::cout << "Read " << nloc << " locations read from \""
93  << locfilename << "\"" << std::endl;
94  for (size_t iloc = 0; iloc < nloc; ++iloc ) {
95  std::cout << "[ " << std::setw(2) << iloc << "] "
96  << dkmeta->location[iloc] << std::endl;
97  }
98 
99  ///-----------------------------------------------------------------------
100  ///
101  /// equivalent to NumiAnalysis::book() in g4numi
102  ///
103  ///-----------------------------------------------------------------------
104 
105  // create file, book tree, set branch address to created object
106  TFile* treeFile = new TFile("test_dk2nu.root","RECREATE");
107 
108  TTree* dk2nuTree = new TTree("dk2nuTree","neutrino ntuple");
109  dk2nuTree->Branch("dk2nu","bsim::Dk2Nu",&dk2nu,32000,99);
110 #ifdef ADD_NONSTD
111  // extend the tree with additional branches without modifying std class
112  dk2nuTree->Branch("nonstdb","NonStd",&nonstd,32000,99);
113 #endif
114 
115  TTree* dkmetaTree = new TTree("dkmetaTree","neutrino ntuple metadata");
116  dkmetaTree->Branch("dkmeta","bsim::DkMeta",&dkmeta,32000,99);
117 
118  int myjob = 42; // unique identifying job # for this series
119 
120  ///-----------------------------------------------------------------------
121  ///
122  /// equivalent to NumiAnalysis::? in g4numi
123  /// this is the main loop, making entries as the come about
124  ///
125  ///-----------------------------------------------------------------------
126  // fill a few element of a few entries
127  for (unsigned int ipot=1; ipot <= nentries; ++ipot) {
128 
129  ///
130  /// equivalent to NumiAnalysis::FillNeutrinoNtuple() in g4numi
131  /// (only the part within the loop over ipot)
132  ///
133 
134  // clear the object in preparation for filling an entry
135  dk2nu->clear();
136 
137  // fill with info ... only a few elements, just for test purposes
138  dk2nu->job = myjob;
139  dk2nu->potnum = ipot;
140 
141  // pick a bogus particle type to decay, and a neutrino flavor
142  int ptype = 211; // pi+
143  if ( ipot % 5 == 0 ) ptype = 321; // k+
144  if ( ipot % 50 == 0 ) ptype = 13; // mu-
145  int ntype = ( ( ptype == 321 ) ? 12 : 14 );
146  TVector3 p3nu(1,2,3); // bogus random neutrino decay vector
147 
148  // calcLocationWeights needs these filled if it isn't going assert()
149  // really need to fill the other bits at this point as well:
150  // ntype, ptype, vx, vy, vz, pdpx, pdpy, pdpz, necm,
151  // ppenergy, ppdxdz, ppdydz, pppz,
152  // muparpx, muparpy, muparpz, mupare
153  dk2nu->decay.ptype = ptype;
154  dk2nu->decay.ntype = ntype;
155 
156  // fill nupx, nupy, nupz, nuenergy, nuwgt(=1) for random decay
157  // should be the 0-th entry
158  if ( dkmeta->location[0].name == "random decay" ) {
159  bsim::NuRay nurndm(p3nu.x(),p3nu.y(),p3nu.z(),p3nu.Mag(),1.0);
160  dk2nu->nuray.push_back(nurndm);
161  }
162  // fill location specific p3, energy and weights; locations in metadata
163  bsim::calcLocationWeights(dkmeta,dk2nu);
164 
165  // test the filling of vector where entries vary in length
166  // ... really need to fill whole dk2nu object
167  unsigned int nancestors = rndm->Integer(12) + 1; // at least one entry
168  for (unsigned int janc = 0; janc < nancestors; ++janc ) {
169  int xpdg = rndm->Integer(100);
170  bsim::Ancestor anc;
171  anc.pdg = janc*10000+xpdg;
172  dk2nu->ancestor.push_back(anc);
173  }
174 
175  // push a couple of user defined values for each entry
176  dk2nu->vint.push_back(42);
177  dk2nu->vint.push_back(ipot);
178 
179 #ifdef ADD_NONSTD
180  // fill non-standard extension to tree with user additions
181  nonstd->foo = ptype + 1000000;
182  nonstd->bar = ipot + ptype;
183 #endif
184 
185  // push entry out to tree
186  dk2nuTree->Fill();
187 
188  } // end of fill loop
189 
190  ///-----------------------------------------------------------------------
191  ///
192  /// equivalent to NumiAnalysis::finish() in g4numi
193  ///
194  ///-----------------------------------------------------------------------
195 
196  /// fill the rest of the metadata (locations filled above)
197  //no! would clear location info // dkmeta->Clear();
198  dkmeta->job = myjob; // needs to match the value in each dk2nu entry
199  dkmeta->pots = 50000; // ntuple represents this many protons-on-target
200  dkmeta->beamsim = "test_fill_dk2nu.C";
201  dkmeta->physics = "bogus";
202  dkmeta->vintnames.push_back("mytemp_42");
203  dkmeta->vintnames.push_back("mytemp_ipot");
204  // push entry out to meta-data tree
205  dkmetaTree->Fill();
206 
207  // finish and clean-up
208  treeFile->cd();
209  dk2nuTree->Write();
210  dkmetaTree->Write();
211  treeFile->Close();
212  delete treeFile; treeFile=0;
213  dk2nuTree=0;
214  dkmetaTree=0;
215 
216 }
std::vector< bsim::Location > location
locations
Definition: dkmeta.h:119
Int_t job
identifying job #
Definition: dk2nu.h:323
std::vector< std::string > vintnames
names of elements for user defined vector of integers
Definition: dkmeta.h:126
TTree * dk2nuTree
void clear(const std::string &opt="")
reset everything
std::vector< bsim::NuRay > nuray
rays through detector fixed points
Definition: dk2nu.h:328
Int_t potnum
proton # processed by simulation
Definition: dk2nu.h:324
TTree * dkmetaTree
double bar
more data
void readWeightLocations(std::string locfilename, std::vector< bsim::Location > &locations)
bsim::Decay decay
basic decay information
Definition: dk2nu.h:327
int myjob
double foo
data for my one-off test
Int_t job
identifying job # (keep files distinct)
Definition: dkmeta.h:86
bsim::Dk2Nu * dk2nu
TRandom3 * rndm
Double_t pots
protons-on-target
Definition: dkmeta.h:87
Int_t pdg
ancestor species
Definition: dk2nu.h:188
TFile * treeFile
void test_fill_dk2nu(unsigned int nentries=1000)
Long64_t nentries
std::string physics
e.g. "fluka08", "g4.9.3p01"
Definition: dkmeta.h:96
std::string beamsim
e.g. "flugg" or "g4numi/<tag>"
Definition: dkmeta.h:95
Int_t ntype
% neutrino flavor (PDG? code)
Definition: dk2nu.h:128
std::vector< Int_t > vint
user defined vector of integers
Definition: dk2nu.h:348
virtual ~NonStd()
OStream cout
Definition: OStream.cxx:6
bsim::DkMeta * dkmeta
void Clear()
void calcLocationWeights(const bsim::DkMeta *dkmeta, bsim::Dk2Nu *dk2nu)
user interface
std::vector< bsim::Ancestor > ancestor
chain from proton to neutrino
Definition: dk2nu.h:329
Int_t ptype
% nu parent species (PDG? code)
Definition: dk2nu.h:144
enum BeamMode string