ROOTPWA
testLikelihood.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2010
4 //
5 // This file is part of rootpwa
6 //
7 // rootpwa is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
11 //
12 // rootpwa is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
16 //
17 // You should have received a copy of the GNU General Public License
18 // along with rootpwa. If not, see <http://www.gnu.org/licenses/>.
19 //
21 //-------------------------------------------------------------------------
22 // File and Version Information:
23 // $Rev:: $: revision of last commit
24 // $Author:: $: author of last commit
25 // $Date:: $: date of last commit
26 //
27 // Description:
28 // generates tree with log likelihood values and its derivatives
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <map>
39 #include <list>
40 
41 #include "TROOT.h"
42 #include "TRandom3.h"
43 #include "TFile.h"
44 #include "TTree.h"
45 
46 #include "reportingUtils.hpp"
47 #include "TPWALikelihood.h"
48 
49 
50 using namespace std;
51 using namespace rpwa;
52 
53 
54 void
55 usage(const string& progName,
56  const int errCode = 0)
57 {
58  cerr << "usage:" << endl
59  << progName
60  << " -w wavelist [-z # val -d amplitude directory -o outfile -s seed -N -n normfile"
61  << " [-a normfile] -r rank -q -h]" << endl
62  << " where:" << endl
63  << " -w file path to wavelist file" << endl
64  << " -z # number of likelihood values (default: 100)" << endl
65  << " -d dir path to directory with decay amplitude files (default: '.')" << endl
66  << " -o file path to output file (default: 'testLikelihood.root')" << endl
67  << " -s # seed for random start values (default: 1234567)" << endl
68  << " -N use normalization of decay amplitudes (default: false)" << endl
69  << " -n file path to normalization integral file (default: 'norm.int')" << endl
70  << " -a file path to acceptance integral file (default: 'norm.int')" << endl
71  << " -A # number of input events to normalize acceptance to" << endl
72  << " -r # rank of spin density matrix (default: 1)" << endl
73  << " -q run quietly (default: false)" << endl
74  << " -h print help" << endl
75  << endl;
76  exit(errCode);
77 }
78 
79 
80 int
81 main(int argc,
82  char** argv)
83 {
84  // force loading of predefined std::map dictionary
85  gROOT->ProcessLine("#include <map>");
86 
87  const string progName = argv[0];
88  unsigned int nmbLikelihoodVals = 100; // number of likelihood values to calculate
89  int startValSeed = 1234567;
90  string waveListFileName = ""; // wavelist filename
91  string ampDirName = "."; // decay amplitude directory name
92  string outFileName = "testLikelihood.root"; // output filename
93  bool useNormalizedAmps = false; // if true normalized amplitudes are used
94  string normIntFileName = ""; // file with normalization integrals
95  string accIntFileName = ""; // file with acceptance integrals
96  unsigned int numbAccEvents = 0; // number of events used for acceptance integrals
97  unsigned int rank = 1; // rank of spin-density matrix
98  bool quiet = false;
99  extern char* optarg;
100  // extern int optind;
101  int c;
102  while ((c = getopt(argc, argv, "z:w:d:o:s:Nn:a:A:r:qh")) != -1)
103  switch (c) {
104  case 'z':
105  nmbLikelihoodVals = atoi(optarg);
106  break;
107  case 'w':
108  waveListFileName = optarg;
109  break;
110  case 'd':
111  ampDirName = optarg;
112  break;
113  case 'o':
114  outFileName = optarg;
115  break;
116  case 's':
117  startValSeed = atoi(optarg);
118  break;
119  case 'N':
120  useNormalizedAmps = true;
121  break;
122  case 'n':
123  normIntFileName = optarg;
124  break;
125  case 'a':
126  accIntFileName = optarg;
127  break;
128  case 'A':
129  numbAccEvents = atoi(optarg);
130  break;
131  case 'r':
132  rank = atoi(optarg);
133  break;
134  case 'q':
135  quiet = true;
136  break;
137  case 'h':
138  usage(progName);
139  break;
140  }
141  if (normIntFileName.length() <= 1) {
142  normIntFileName = "norm.int";
143  printWarn << "using default normalization integral file '" << normIntFileName << "'." << endl;
144  }
145  if (accIntFileName.length() <= 1) {
146  accIntFileName = "norm.int";
147  printWarn << "using default acceptance normalization integral file '" << accIntFileName << "'." << endl;
148  }
149  if (waveListFileName.length() <= 1) {
150  printErr << "no wavelist file specified! aborting!" << endl;
151  usage(progName, 1);
152  }
153  // report parameters
154  printInfo << "running " << progName << " with the following parameters:" << endl;
155  cout << " # of values to calculate ..................... " << nmbLikelihoodVals << endl
156  << " wave list file ............................... '" << waveListFileName << "'" << endl
157  << " output file .................................. '" << outFileName << "'" << endl
158  << " seed for random parameters ................... " << startValSeed << endl
159  << " use normalization ............................ " << useNormalizedAmps << endl
160  << " file with normalization integral ......... '" << normIntFileName << "'" << endl
161  << " file with acceptance integral ............ '" << accIntFileName << "'" << endl
162  << " number of acceptance norm. events ........ " << numbAccEvents << endl
163  << " rank of spin-density matrix .................. " << rank << endl
164  << " quiet ........................................ " << ((quiet) ? "yes" : "no") << endl;
165 
166  // setup likelihood function
167  printInfo << "creating and setting up likelihood function" << endl;
168  TPWALikelihood<complex<double> > L;
169  if (quiet)
170  L.setQuiet();
171  L.useNormalizedAmps(useNormalizedAmps);
172 #ifdef USE_CUDA
173  //L.enableCuda(false);
174  L.enableCuda(true);
175 #endif
176  L.init(rank, waveListFileName, normIntFileName, accIntFileName, ampDirName, numbAccEvents);
177  if (!quiet)
178  cout << L << endl;
179 
180  // create output file and tree
181  TFile& outFile = *TFile::Open(outFileName.c_str(), "RECREATE");
182  TTree tree("testLikelihoodTree", "testLikelihoodTree");
183 
184  // define tree branches
185  Double_t logLikelihood;
186  map<string, double> derivatives;
187  tree.Branch("logLikelihood", &logLikelihood, "logLikelihood/D");
188  tree.Branch("derivatives", &derivatives);
189 
190  // setup random number generator and production amplitude map
191  TRandom3 random(startValSeed);
192  list<string> prodAmpNames; // needed, to ensure that parameter values do not depend on ordering
193  const unsigned int nmbPar = L.NDim();
194  for (unsigned int i = 0; i < nmbPar; ++i)
195  prodAmpNames.push_back(L.parName(i));
196  prodAmpNames.sort();
197 
198  // calculate some likelihood values and derivatives
199  for (unsigned int i = 0; i < nmbLikelihoodVals; ++i) {
200  // generate random production amplitudes independent of parameter order
201  double prodAmps[nmbPar];
202  bool success = true;
203  for (list<string>::const_iterator j = prodAmpNames.begin(); j != prodAmpNames.end(); ++j) {
204  int parIndex = -1;
205  for (unsigned int k = 0; k < nmbPar; ++k)
206  if (L.parName(k) == *j) {
207  parIndex = k;
208  break;
209  }
210  if (parIndex < 0) {
211  cout << "cannot find parameter '" << *j << "'. skipping." << endl;
212  success = false;
213  continue;
214  }
215  prodAmps[parIndex] = random.Uniform(0, 10);
216  // cout << " [" << *j << "] = [" << parIndex << "] = " << prodAmps[parIndex] << endl;
217  }
218  if (success) {
219  logLikelihood = L.DoEval(prodAmps);
220  // cout << "[" << i << "] = " << logLikelihood << endl;
221  double gradient[nmbPar];
222  L.Gradient(prodAmps, gradient);
223  derivatives.clear();
224  for (unsigned int j = 0; j < nmbPar; ++j) {
225  const string parName = L.parName(j);
226  derivatives[parName] = gradient[j];
227  }
228  tree.Fill();
229  }
230  }
231 
232  tree.Write();
233  outFile.Close();
234 
235  return 0;
236 }
Generated by   doxygen 1.8.1.2