ROOTPWA
genpw.cc
Go to the documentation of this file.
1 
2 //
3 // Copyright 2009 Sebastian Neubert
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 
26 #include <iostream>
27 #include <fstream>
28 #include <iomanip>
29 #include <complex>
30 #include <cstdlib>
31 #include <string>
32 #include <vector>
33 #include <unistd.h>
34 #include <stdlib.h>
35 #include <event.h>
36 
37 #include "TFile.h"
38 #include "TTree.h"
39 #include "TString.h"
40 #include "TH1.h"
41 #include "TRandom3.h"
42 #include "TLorentzVector.h"
43 #include "TVector3.h"
44 #include "TClonesArray.h"
45 
46 #include "libconfig.h++"
47 
48 #include "reportingUtils.hpp"
49 #include "diffractivePhaseSpace.h"
50 #include "partialWaveWeight.h"
51 #include "productionAmp.h"
52 #include "breitWignerProductionAmp.h"
53 #include "TFitBin.h"
54 
55 #include "primaryVertexGen.h"
56 
57 
58 using namespace std;
59 using namespace libconfig;
60 using namespace rpwa;
61 
62 extern particleDataTable PDGtable;
63 
64 
65 void printUsage(char* prog, int errCode = 0)
66 {
67  cerr << "usage:" << endl
68  << prog
69  << " -n # [-a # -m # -M # -B # -s #] -o <file> -w <file> -k <path> -i <file> -r <file>" << endl
70  << " where:" << endl
71  << " -n # (max) number of events to generate (default: 100)" << endl
72  << " -a # (max) number of attempts to do (default: infinity)" << endl
73  << " -m # maxWeight" << endl
74  << " -o <file> ROOT output file (if not specified, generated automatically)" << endl
75  << " -w <file> wavelist file (contains production amplitudes)" << endl
76  << " -w <file.root> to use TFitBin tree as input" << endl
77  << " -c <0/1> if 1 a comgeant eventfile (.fort.26) is written with same naming as the root file (default 1)" << endl
78  << " -k <path> path to keyfile directory (all keyfiles have to be there)" << endl
79  << " -i <file> integral file" << endl
80  << " -r <file> reaction config file" << endl
81  << " -s # set seed " << endl
82  << " -M # lower boundary of mass range in MeV (overwrites values from config file) " << endl
83  << " -B # width of mass bin in MeV" << endl
84  << endl;
85  exit(errCode);
86 }
87 
88 
89 int main(int argc, char** argv)
90 {
91 
92  unsigned int nevents = 100;
93  unsigned int max_attempts = 0;
94  string output_file = ""; // either given by option or generated automatically by mass range
95  string output_evt = "";
96  string output_wht = "";
97  string output_comgeant = "";
98  string integrals_file;
99  bool hasint = false;
100  string wavelist_file; // format: name Re Im
101  string path_to_keyfiles = "./";
102  string reactionFile;
103  double maxWeight = 0;
104  int seed = 123456;
105  int massLower = 0;
106  int massBinWidth = 0;
107  bool overwriteMass = false;
108  bool writeComGeantout = false;
109 
110  int c;
111  while ((c = getopt(argc, argv, "n:a:o:w:k:i:r:m:s:M:B:h:c")) != -1) {
112  switch (c) {
113  case 'n':
114  nevents = atoi(optarg);
115  break;
116  case 'a':
117  max_attempts = atoi(optarg);
118  break;
119  case 's':
120  seed = atoi(optarg);
121  break;
122  case 'o':
123  output_file = optarg;
124  break;
125  case 'w':
126  wavelist_file = optarg;
127  break;
128  case 'i':
129  integrals_file = optarg;
130  hasint=true;
131  break;
132  case 'k':
133  path_to_keyfiles = optarg;
134  break;
135  case 'r':
136  reactionFile = optarg;
137  break;
138  case 'm':
139  maxWeight = atof(optarg);
140  break;
141  case 'c':
142  writeComGeantout = true;
143  break;
144  case 'M':
145  massLower = atoi(optarg);
146  overwriteMass=true;
147  break;
148  case 'B':
149  massBinWidth = atoi(optarg);
150  overwriteMass=true;
151  break;
152 
153  case 'h':
154  printUsage(argv[0]);
155  break;
156  default:
157  printUsage(argv[0]);
158  break;
159  }
160  }
161 
162  PDGtable.initialize();
163 
164  partialWaveWeight weighter;
165  Config reactConf;
166  reactConf.readFile(reactionFile.c_str());
167 
168  // variable that need to get initialized either by input options
169  // or the config file
170  // will be stored in the tree later
171  double weight, impweight;
172  TClonesArray* p = new TClonesArray("TLorentzVector");
173  TLorentzVector beam;
174  TVector3 vertex;
175  double tprime = 0.;
176  int qbeam;
177  vector<int> q; // array of charges
178 
179  double Mom = reactConf.lookup("beam.momentum");
180  double MomSigma = reactConf.lookup("beam.sigma_momentum");
181  double BeamPartMass = 0.13957018;
182  if(reactConf.exists("beam.mass")){
183  BeamPartMass = reactConf.lookup("beam.mass");
184  }
185  double DxDz = reactConf.lookup("beam.DxDz");
186  double DxDzSigma = reactConf.lookup("beam.sigma_DxDz");
187  double DyDz = reactConf.lookup("beam.DyDz");
188  double DyDzSigma = reactConf.lookup("beam.sigma_DyDz");
189 
190  double targetz = reactConf.lookup("target.pos.z");
191  double targetd = reactConf.lookup("target.length");
192  double targetr = reactConf.lookup("target.radius");
193  double mrecoil = reactConf.lookup("target.mrecoil");
194 
195  double tprime_min(0.);
196  double tprime_max(numeric_limits<double>::max());
197  reactConf.lookupValue("finalstate.t_min", tprime_min);
198  reactConf.lookupValue("finalstate.t_max", tprime_max);
199 
200  double mmin = reactConf.lookup("finalstate.mass_min");
201  double mmax = reactConf.lookup("finalstate.mass_max");
202  if(overwriteMass) {
203  mmin = massLower/1000.0;
204  mmax = (massLower+massBinWidth)/1000.0;
205  }
206  // array of tslopes even when only one is existing
207  double* tslope = NULL;
208  double* inv_m = NULL;
209  int ntslope = 1;
210  if(reactConf.lookup("finalstate.t_slope").isArray()) {
211  ntslope = reactConf.lookup("finalstate.t_slope").getLength();
212  if(reactConf.lookup("finalstate.inv_m").getLength() != ntslope) {
213  cout << " Error: please check number of t' values and the corresponding invariant masses in the Configuration File! " << endl;
214  return 0;
215  }
216  tslope = new double[ntslope];
217  inv_m = new double[ntslope];
218  cout << " found array of t' slopes. Reading " << ntslope << " of values ";
219  for(int i = 0; i < ntslope; i++) {
220  tslope[i] = reactConf.lookup("finalstate.t_slope")[i];
221  inv_m[i] = reactConf.lookup("finalstate.inv_m")[i];
222  //cout << inv_m[i] << " " << tslope[i] << endl;
223  }
224  cout << " done. " << endl;
225  } else {
226  tslope = new double[1];
227  tslope[0] = reactConf.lookup("finalstate.t_slope");
228  //cout << " tslope set to " << tslope[0];
229  }
230  double binCenter = 500 * (mmin + mmax);
231 
232  // check weather to use a primary vertex generator as requested by the config file
233  primaryVertexGen* primaryVtxGen = NULL;
234  string histfilename_primvertex = "";
235  if(reactConf.lookupValue("primvertex.histfilename", histfilename_primvertex)) {
236  primaryVtxGen = new primaryVertexGen(
237  histfilename_primvertex,
238  BeamPartMass,
239  Mom,
240  MomSigma
241  );
242  if (!primaryVtxGen->check()) {
243  cerr << " Error: histogram filename with beam properties not loaded! " << endl;
244  delete primaryVtxGen;
245  primaryVtxGen = NULL;
246  }
247  }
248 
249  if(!reactConf.lookupValue("beam.charge",qbeam)) {
250  qbeam = -1;
251  }
252  // generate the filename automatically if not specified
253  if (output_file == "") {
254  stringstream _filename;
255  _filename << massLower << "." << massLower+massBinWidth << ".genbod.root";
256  cout << " created output filename: " << _filename.str() << endl;
257  output_file = _filename.str();
258  }
259  output_evt = output_file;
260  output_evt.replace(output_evt.find(".root"),5,".evt");
261  output_wht = output_file;
262  output_wht.replace(output_wht.find(".root"),5,".wht");
263  output_comgeant = output_file;
264  output_comgeant.erase(output_comgeant.find(".root"),5);
265  output_comgeant.append(".fort.26");
266 
267  // now create the root file to store the events
268  TFile* outfile = TFile::Open(output_file.c_str(),"RECREATE");
269  TH1D* hWeights = new TH1D("hWeights","PW Weights",100,0,100);
270  TTree* outtree = new TTree("pwevents","pwevents");
271 
272  outtree->Branch("weight",&weight,"weight/d");
273  outtree->Branch("impweight",&impweight,"impweight/d");
274  outtree->Branch("p",&p);
275  outtree->Branch("beam",&beam);
276  outtree->Branch("vertex", &vertex);
277  outtree->Branch("q",&q);
278  outtree->Branch("qbeam",&qbeam,"qbeam/I");
279  outtree->Branch("tprime", &tprime);
280 
281  //string theta_file= reactConf.lookup("finalstate.theta_file");
282 
283  diffractivePhaseSpace difPS;
284  cerr << "Seed=" << seed << endl;
285  difPS.SetSeed(seed);
286  difPS.SetBeam(Mom,MomSigma,DxDz,DxDzSigma,DyDz,DyDzSigma);
287  difPS.SetTarget(targetz,targetd,targetr,mrecoil);
288  difPS.SetTPrimeSlope(tslope, inv_m, ntslope);
289  difPS.SetMassRange(mmin,mmax);
290  difPS.SetPrimaryVertexGen(primaryVtxGen);
291  if(tprime_min >= 0.) {
292  difPS.SettprimeMin(tprime_min);
293  } else {
294  cout << " Error: tprime_min must be positive " << tprime_min << endl;
295  }
296  if(tprime_max >= 0.) {
297  difPS.SettprimeMax(tprime_max);
298  }
299 
300  double impMass;
301  double impWidth;
302 
303  if(reactConf.lookupValue("importance.mass",impMass) &&
304  reactConf.lookupValue("importance.width",impWidth))
305  {
306  int act = reactConf.lookup("importance.active");
307  if(act == 1) {
308  difPS.SetImportanceBW(impMass,impWidth);
309  }
310  }
311 
312  const Setting& root = reactConf.getRoot();
313  const Setting& fspart = root["finalstate"]["particles"];
314  int nparticles = fspart.getLength();
315 
316  for(int ifs=0; ifs<nparticles; ++ifs) {
317  const Setting &part = fspart[ifs];
318  int id;
319  part.lookupValue("g3id",id);
320  int myq;
321  part.lookupValue("charge",myq);
322  double m;
323  part.lookupValue("mass",m);
324  q.push_back(myq);
325  difPS.AddDecayProduct(particleInfo(id,myq,m));
326  }
327 
328  // see if we have a resonance in this wave
329  map<string, breitWignerProductionAmp*> bwAmps;
330  if(reactConf.exists("resonances")) {
331  const Setting &bws = root["resonances"]["breitwigners"];
332  // loop through breitwigners
333  int nbw = bws.getLength();
334  printInfo << "found " << nbw << " Breit-Wigners in config" << endl;
335 
336  for(int ibw = 0; ibw < nbw; ++ibw) {
337  const Setting &bw = bws[ibw];
338  string jpcme;
339  double mass, width;
340  double cRe, cIm;
341  bw.lookupValue("jpcme", jpcme);
342  bw.lookupValue("mass", mass);
343  bw.lookupValue("width", width);
344  bw.lookupValue("coupling_Re", cRe);
345  bw.lookupValue("coupling_Im", cIm);
346  complex<double> coupl(cRe, cIm);
347  cout << " JPCME = " << jpcme << ", mass = " << mass << " GeV/c^2, "
348  << "width = " << width << " GeV/c^2, coupling = " << coupl << endl;
349  bwAmps[jpcme] = new breitWignerProductionAmp(mass, width, coupl);
350  }
351  }
352 
353  // check if TFitBin is used as input
354  if(wavelist_file.find(".root") != string::npos) {
355  cerr << "Using TFitBin as input!" << endl;
356  TFile* fitresults = TFile::Open(wavelist_file.c_str(),"READ");
357  TFitBin* Bin = NULL;
358  if(!fitresults || fitresults->IsZombie()) {
359  cerr << "Cannot open start fit results file " << wavelist_file << endl;
360  return 1;
361  }
362  // get tree with start values
363  TTree* tree;
364  fitresults->GetObject("pwa", tree);
365  if(!tree) {
366  cerr << "Cannot find fitbin tree '" << "pwa" << "' " << endl;
367  } else {
368  Bin = new TFitBin();
369  tree->SetBranchAddress("fitbin", &Bin);
370  // find entry which is closest to mass bin center
371  unsigned int iBest = 0;
372  double mBest = 0;
373  for(unsigned int i = 0; i < tree->GetEntriesFast(); ++i) {
374  tree->GetEntry(i);
375  if(fabs(binCenter - Bin->mass()) <= fabs(binCenter - mBest)) {
376  iBest = i;
377  mBest = Bin->mass();
378  }
379  } // end loop over TFitBins
380  cerr << "Using data from Mass bin m=" << mBest << endl;
381  tree->GetEntry(iBest);
382  // write wavelist file for generator
383  string tmpname = tmpnam(NULL);
384  ofstream tmpfile(tmpname.c_str());
385  Bin->printAmpsGenPW(tmpfile);
386  tmpfile.close();
387  wavelist_file=tmpname;
388  }
389  } // end root file given
390 
391  // read input wavelist and amplitudes
392  ifstream wavefile(wavelist_file.c_str());
393  while(wavefile.good()) {
394  TString wavename;
395  double RE, IM;
396  int rank = 0;
397  int refl = 0;
398  wavefile >> wavename >> RE >> IM;
399 
400  if(wavename.Contains("flat") || wavename.Length()<2) {
401  continue;
402  }
403  if(RE == 0 && IM == 0) {
404  continue;
405  }
406  // check if there is rank information
407  if(wavename(0) == 'V') {
408  // we multiply rank by to to make space for refl+- production vectors
409  rank = 2 * atoi(wavename(1,1).Data());
410  // check reflecitivity to sort into correct production vector
411  refl = wavename(9)=='+' ? 0 : 1;
412  wavename = wavename(3, wavename.Length());
413  }
414 
415  TString jpcme = wavename(2,5);
416 
417  wavename.ReplaceAll(".amp", ".key");
418  wavename.Prepend(path_to_keyfiles.c_str());
419 
420  complex<double> amp(RE, IM);
421  cerr << wavename << " " << amp << " r=" << rank/2
422  << " eps=" << refl << " qn=" << jpcme << endl;
423  wavefile.ignore(256, '\n');
424 
425  productionAmp* pamp;
426  if(bwAmps[jpcme.Data()] != NULL) {
427  pamp = bwAmps[jpcme.Data()];
428  cerr << "Using BW for " << jpcme << endl;
429  // production vector index: rank+refl
430  weighter.addWave(wavename.Data(), pamp, amp, rank+refl);
431  } else {
432  pamp = new productionAmp(amp);
433  weighter.addWave(wavename.Data(), pamp, complex<double>(1, 0), rank+refl);
434  }
435  }
436 
437  if(hasint){
438  // read integral files
439  ifstream intfile(integrals_file.c_str());
440  while(intfile.good()) {
441  string filename;
442  double mass;
443  intfile >> filename >> mass;
444  weighter.loadIntegrals(filename, mass);
445  intfile.ignore(256, '\n');
446  } // loop over integralfile
447  }// endif hasint
448 
449 
450  double maxweight = -1;
451  unsigned int attempts = 0;
452 
453 
454  unsigned int tenpercent = (unsigned int)(nevents / 10);
455  unsigned int i = 0;
456  //difPS.setVerbose(true);
457  ofstream evtout(output_evt.c_str());
458  ofstream evtgeant;
459  if(writeComGeantout) {
460  evtgeant.open(output_comgeant.c_str());
461  }
462  ofstream evtwht(output_wht.c_str());
463  evtwht << setprecision(10);
464 
465  while(i < nevents && ((max_attempts > 0 && attempts < max_attempts) || max_attempts == 0))
466  {
467 
468  ++attempts;
469 
470  p->Delete(); // clear output array
471 
472  stringstream str;
473  if (writeComGeantout) {
474  difPS.event(str, evtgeant);
475  } else {
476  difPS.event(str);
477  }
478  impweight = difPS.impWeight();
479 
480  for(int ip = 0; ip < nparticles; ++ip){
481  new ((*p)[ip]) TLorentzVector(*difPS.GetDecay(ip));
482  }
483 
484  beam = *difPS.GetBeam();
485  vertex = *difPS.GetVertex();
486  tprime = difPS.Gettprime();
487 
488  // calculate weight
489  event e;
490  e.setIOVersion(1);
491 
492  str >> e;
493  evtout << e;
494  evtwht << impweight << endl;
495 
496  //cerr << e << endl;
497 
498  weight = weighter.weight(e);
499  if(weight > maxweight) {
500  maxweight = weight;
501  }
502  hWeights->Fill(weight);
503 
504  if(maxWeight > 0) { // do weighting
505  cout << weight << endl;
506  //if(weight>maxWeight)maxWeight=weight;
507  if(gRandom->Uniform() > weight / maxWeight) {
508  continue;
509  }
510  }
511  //cerr << i << endl;
512 
513  outtree->Fill();
514 
515 
516  if(i > 0 && (i % tenpercent == 0)) {
517  cerr << "\x1B[2K" << "\x1B[0E" << "[" << (double)i/(double)nevents*100. << "%]";
518  }
519 
520  ++i;
521 
522  } // end event loop
523 
524  cerr << endl;
525  cerr << "Maxweight: " << maxweight << endl;
526  cerr << "Attempts: " << attempts << endl;
527  cerr << "Created Events: " << i << endl;
528  cerr << "Efficiency: " << (double)i / (double)attempts << endl;
529 
530 
531  outfile->cd();
532  hWeights->Write();
533  outtree->Write();
534  outfile->Close();
535  evtout.close();
536  evtwht.close();
537  if(writeComGeantout) {
538  evtgeant.close();
539  }
540 
541  delete [] tslope;
542  delete [] inv_m;
543 
544  return 0;
545 
546 }
547 
Generated by   doxygen 1.8.1.2