ROOTPWA
evtweight.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 <complex>
27 #include <iostream>
28 
29 #include <fstream>
30 #include <cstdlib>
31 #include <string>
32 #include <vector>
33 #include <unistd.h>
34 #include <stdlib.h>
35 
36 #include "TClonesArray.h"
37 #include "TFile.h"
38 #include "TH1.h"
39 #include "TH1D.h"
40 #include "TLorentzVector.h"
41 #include "TRandom3.h"
42 #include "TString.h"
43 #include "TTree.h"
44 
45 #include "integral.h"
46 
47 #include "fitResult.h"
48 #include "event.h"
49 
50 #include "../amplitude/particle.h"
51 #include "particleDataTable.h"
52 #include "libconfig.h++"
53 
54 
55 using namespace std;
56 using namespace libconfig;
57 using namespace rpwa;
58 
59 void printUsage(char* prog, int errCode=0) {
60  cerr << "usage:" << endl
61  << prog
62  << " -e <file> -o <file> -w <file> -i <file> -n samples -m mass"
63  << " where:" << endl
64  << " -e <file> acc or ps events in .evt or .root format"<< endl
65  << " -o <file> ROOT output file"<< endl
66  << " -w <file.root> use TFitBin tree as input"<< endl
67  << " -i <file> integral file"<< endl
68  << " -m mass center of mass bin"<< endl
69  << " -n samples sampling of the model parameters (default=1)"<< endl
70  << " -b width width of mass bin"<< endl
71  << endl;
72  exit(errCode);
73 }
74 
75 
76 int getReflectivity(const TString& waveName)
77 {
78  int refl = 0;
79  unsigned int reflIndex = 6; // position of reflectivity in wave
80  // check whether it is parameter or wave name
81  if(waveName[0] == 'V') {
82  reflIndex = 9;
83  }
84  if(waveName[reflIndex] == '-') {
85  refl = -1;
86  } else if(waveName[reflIndex] == '+') {
87  refl = +1;
88  } else {
89  printErr << "Cannot parse parameter/wave name '" << waveName << "'. Cannot not determine reflectivity. Aborting." << endl;
90  throw;
91  }
92  return refl;
93 }
94 
95 
96 void parseWaveList(const string& waveListFileName,
97  vector<string>& waveNames,
98  vector<double>& waveThresholds,
99  vector<int>& refl)
100 {
101 
102  unsigned int _nmbWavesPosRefl = 0;
103  unsigned int _nmbWavesNegRefl = 0;
104  bool _debug = false;
105  printInfo << "Reading amplitude names from wave list file '" << waveListFileName << "'." << endl;
106  ifstream waveListFile(waveListFileName.c_str());
107  if(!waveListFile) {
108  printErr << "Cannot open file '" << waveListFileName << "'. Exiting." << endl;
109  throw;
110  }
111  string line;
112  int lineNmb = 0;
113  while(getline(waveListFile, line)) {
114  if(line[0] == '#') { // comments start with #
115  continue;
116  }
117  stringstream lineStream;
118  lineStream.str(line);
119  string waveName;
120  if(lineStream >> waveName) {
121  double threshold;
122  // !!! it would be safer to make the threshold value in the wave list file mandatory
123  if(!(lineStream >> threshold))
124  threshold = 0;
125  if(_debug)
126  cout << " reading line " << setw(3) << lineNmb + 1 << ": " << waveName << ", threshold = " << setw(4) << threshold << " MeV/c^2" << endl;
127  if(getReflectivity(waveName) > 0){
128  ++_nmbWavesPosRefl; refl.push_back(1);
129  } else {
130  ++_nmbWavesNegRefl; refl.push_back(-1);
131  }
132  waveNames.push_back(waveName);
133  waveThresholds.push_back(threshold);
134  } else {
135  printWarn << "Cannot parse line '" << line << "' in wave list file '" << waveListFileName << "'." << endl;
136  }
137  ++lineNmb;
138  }
139  waveListFile.close();
140  printInfo << "Read " << lineNmb << " lines from wave list file '" << waveListFileName << "'." << endl;
141 }
142 
143 
144 int main(int argc, char** argv)
145 {
146 
147  if(argc < 3) {
148  printUsage(argv[0],1);
149  }
150 
151  string output_file = "genpw.root";
152  string integrals_file;
153  string waveListFileName = "wavelist";
154  string evtfilename;
155  double binCenter = 0;
156  double binWidth = 60; // MeV
157  unsigned int nsamples = 1;
158 
159  int c;
160  while ((c = getopt(argc, argv, "e:o:w:i:m:n:h")) != -1) {
161  switch (c) {
162  case 'e':
163  evtfilename = optarg;
164  break;
165  case 'o':
166  output_file = optarg;
167  break;
168  case 'w':
169  waveListFileName = optarg;
170  break;
171  case 'i':
172  integrals_file = optarg;
173  break;
174  case 'm':
175  binCenter = atof(optarg);
176  break;
177  case 'n':
178  nsamples = atoi(optarg);
179  break;
180  case 'b':
181  binWidth = atof(optarg);
182  break;
183 
184  case 'h':
185  printUsage(argv[0]);
186  break;
187  default:
188  printUsage(argv[0]);
189  break;
190  }
191  }
192 
193  binCenter += 0.5 * binWidth;
194 
195  TFile* outfile = TFile::Open(output_file.c_str(), "RECREATE");
196  TH1D* hWeights = new TH1D("hWeights", "PW Weights", 100, 0, 100);
197  TTree* outtree = new TTree("pwevents", "pwevents");
198  double weight;
199  TClonesArray* p = new TClonesArray("TLorentzVector");
200  TLorentzVector beam;
201  double qbeam;
202  std::vector<int> q; // array of charges
203 
204  outtree->Branch("weight", &weight, "weight/d");
205  outtree->Branch("p", &p);
206  outtree->Branch("beam", &beam);
207  outtree->Branch("q", &q);
208  outtree->Branch("qbeam", &qbeam, "qbeam/i");
209 
210  // load integrals ---------------------------------------------------
211  integral normInt;
212  ifstream intFile(integrals_file.c_str());
213  if(!intFile) {
214  printErr << "Cannot open file '"
215  << integrals_file << "'. Exiting." << endl;
216  throw;
217  }
218  // !!! integral.scan() performs no error checks!
219  normInt.scan(intFile);
220  intFile.close();
221 
222  // load production amplitudes ------------------------------------------
223  // read TFitResult is used as input
224  TFile* fitresults = TFile::Open(waveListFileName.c_str(), "READ");
225  fitResult* Bin = NULL;
226  if(!fitresults || fitresults->IsZombie()) {
227  cerr << "Cannot open start fit results file " << waveListFileName << endl;
228  return 1;
229  }
230  // get tree with start values
231  bool hasfit = true;
232  TTree* tree;
233  fitresults->GetObject("pwa", tree);
234  if(!tree) {
235  cerr << "Cannot find fitbin tree '"<< "pwa" << "' "<< endl;
236  } else {
237  Bin = new fitResult();
238  tree->SetBranchAddress("fitResult_v2", &Bin);
239  // find entry which is closest to mass bin center
240  // and has the lowest likelihood
241  unsigned int iBest = 0;
242  double mBest = 0;
243  //double loglike = 0;
244  for(unsigned int i = 0; i < tree->GetEntriesFast(); ++i) {
245  tree->GetEntry(i);
246  if (fabs(binCenter - Bin->massBinCenter()) <= fabs(binCenter - mBest)) {
247  // check also if this bin is more likely in case of many fits per bin
248  //if (loglike == 0 || Bin->logLikelihood() < loglike){
249  iBest = i;
250  mBest = Bin->massBinCenter();
251  //loglike = Bin->logLikelihood();
252  //}
253  //else cerr << "This is a redundant fit" << endl;
254  }
255  } // end loop over TFitBins
256  cerr << "mBest= " << mBest << endl;
257 
258 
259  if((mBest < (binCenter-binWidth / 2.)) || (mBest > (binCenter + binWidth / 2.))) {
260  cerr << "No fit found for Mass bin m=" << binCenter << endl;
261  Bin->reset();
262  hasfit = false;
263  return 1;
264  } else {
265  cerr << "Using data from Mass bin m=" << mBest << " bin " << iBest << endl;
266  tree->GetEntry(iBest);
267  }
268  // write wavelist files for generator
269  // string tmpname("genamps.txt");
270  // create slightly modified versions of the model according to covariances
271  for(unsigned int isamples = 0; isamples < nsamples; ++isamples) {
272  // enumerate genamps files
273  TString tmpname("genamps");
274  tmpname+=isamples;
275  tmpname+=".txt";
276  ofstream tmpfile(tmpname.Data());
277  if(isamples == 0) {
278  Bin->printAmpsGenPW(tmpfile);
279  } else {
280  fitResult* clone = Bin->cloneVar();
281  clone->printAmpsGenPW(tmpfile);
282  delete clone;
283  }
284  tmpfile.close();
285  waveListFileName = "genamps0.txt";
286  } // end loop over model versions
287  }
288 
289  vector<string> waveNames;
290 
291  vector<vector<complex<double> >*> prodAmps(nsamples); //production amplitudes
292  for(unsigned int isamples = 0; isamples < nsamples; ++isamples) {
293  prodAmps[isamples] = new vector<complex<double> >();
294  }
295 
296  vector<int> reflectivities;
297  vector<int> ms;
298  vector<int> ranks;
299  int maxrank = 0;
300  // if we have read in a TFitResult the the name of the file has been changed!
301  // so it is ok to use the same variable here. See above!
302 
303  // read in a list of waveListFiles
304  for(unsigned int isamples = 0; isamples < nsamples; ++isamples) {
305  TString tmpname("genamps");
306  tmpname += isamples;
307  tmpname += ".txt";
308  cerr << "Reading back file " << tmpname << endl;
309  ifstream wavefile(tmpname);
310  while(wavefile.good()) {
311  TString wavename;
312  double RE, IM;
313  int rank = 0;
314  int refl = 0;
315  int m = 0;
316  wavefile >> wavename >> RE >> IM;
317  //cerr << wavename << endl;
318 
319  if(wavename.Length() < 2) {
320  continue;
321  }
322  if((RE == 0) && (IM == 0)) {
323  continue;
324  }
325  // check if there is rank information
326  if(wavename(0) == 'V') {
327  // we multiply rank by to to make space for refl+- production vectors
328  rank = 2 * atoi(wavename(1, 1).Data());
329  // check reflecitivity to sort into correct production vector
330  refl = wavename(9)=='+' ? 0 : 1;
331  m = wavename(8)=='0' ? 0 : 1;
332  //cerr << wavename(9) << endl;
333  wavename = wavename(3, wavename.Length());
334  } else if(wavename != "flat") {
335  refl = wavename(6)=='+' ? 0 : 1;
336  m = wavename(5)=='0' ? 0 : 1;
337  //cerr << wavename(6) << endl;
338  }
339 
340  std::complex<double> amp(RE, IM);
341  prodAmps[isamples]->push_back(amp);
342  cerr << wavename << " " << amp << " r=" << rank/2
343  << " eps=" << refl
344  << " m=" << m << endl;
345  wavefile.ignore(256, '\n');
346  // for first file store info on waves
347  if(isamples == 0) {
348  waveNames.push_back(wavename.Data());
349  reflectivities.push_back(refl);
350  ms.push_back(m);
351  ranks.push_back(rank);
352  if(maxrank < rank) {
353  maxrank = rank;
354  }
355  }
356  } // loop over file
357 
358  cerr << "Rank of fit was:" << maxrank + 1 << endl;
359  } // end loop over samples
360 
361  unsigned int nmbWaves = waveNames.size();
362  // TODO reserve list of wheight branches!
363 
364  vector<ifstream*> ampfiles;
365  // reserve vector beforehand because Branch
366  // will take a pointer onto the elements
367  //vector<double> weights((nmbWaves+1)*nmbWaves/2);
368  vector<double> weights(nmbWaves);
369  //unsigned int wcount=0;
370  //create wheight vectors for individual intensities and interference terms
371  // for(unsigned int iw=0;iw<nmbWaves;++iw){
372  // for(unsigned int jw=iw;jw<nmbWaves;++jw){
373  // TString weightname("W_");
374  // if(iw==jw)weightname+=waveNames[iw];
375  // else weightname+=waveNames[iw] +"_"+ waveNames[jw];
376 
377  // outtree->Branch(weightname.Data(),&weights[wcount++],(weightname+"/d").Data());
378  // }
379  // }
380  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
381  TString weightname("Wintens_");
382  weightname += waveNames[iw];
383  outtree->Branch(weightname.Data(), &weights[iw], (weightname + "/d").Data());
384  }
385 
386  // create branches for the weights of the different model variants
387  vector<double> modelweights(nsamples);
388  for(unsigned int isamples = 0; isamples < nsamples; ++isamples) {
389  TString weightname("W");
390  weightname += isamples;
391  outtree->Branch(weightname.Data(), &modelweights[isamples], (weightname + "/d").Data());
392  }
393 
394  // open decay amplitude files --------------------------------------------
395  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
396  ampfiles.push_back(new ifstream(waveNames[iw].c_str()));
397  }
398 
399 
400  // event loop ------------------------------------------------------------
401  unsigned int counter = 0;
402  if(evtfilename.find(".evt") > 0) {
403  event e;
404  list< ::particle> f_mesons;
405  ifstream evtfile(evtfilename.c_str());
406 
407  while(!evtfile.eof() && evtfile.good()) {
408  // evt2tree
409  if(counter % 1000 == 0) {
410  std::cout << ".";
411  }
412  if(counter++ % 10000 == 0) {
413  std::cout << counter;
414  }
415  evtfile >> e;
416  p->Delete(); // clear output arrays
417  q.clear();
418  f_mesons = e.f_mesons();
419  fourVec pX;
420  list< ::particle>::iterator it = f_mesons.begin();
421  while(it != f_mesons.end()) {
422  pX = it->get4P();
423  new ((*p)[p->GetEntries()]) TLorentzVector(pX.x(), pX.y(), pX.z(), pX.t());
424  q.push_back(it->Charge());
425  ++it;
426  }
427  fourVec evtbeam = e.beam().get4P();
428  beam.SetPxPyPzE(evtbeam.x(), evtbeam.y(), evtbeam.z(), evtbeam.t());
429  qbeam = e.beam().Charge();
430 
431 
432  // read decay amps for this event
433  vector<complex<double> > decayamps(nmbWaves);
434  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
435  complex<double> decayamp;
436  ampfiles[iw]->read((char*) &decayamp, sizeof(complex<double> ));
437  decayamps[iw] = decayamp;
438  }
439 
440  // weighting - do this for each model-sample
441  for(unsigned int isample=0; isample<nsamples; ++isample){
442 
443  vector<complex<double> > posm0amps(maxrank + 1); // positive refl vector m=0
444  vector<complex<double> > posm1amps(maxrank + 1); // positive refl vector m=1
445 
446  vector<complex<double> > negm0amps(maxrank + 1); // negative refl vector m=0
447  vector<complex<double> > negm1amps(maxrank + 1); // negative refl vector m=1
448 
449  complex<double> flatamp;
450 
451  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
452  complex<double> decayamp = decayamps[iw];
453  string w1 = waveNames[iw];
454  if(w1 == "flat"){
455  flatamp = prodAmps[isample]->at(iw);
456  continue;
457  }
458  //cerr << w1 << " " << decayamp << endl;
459  double nrm = sqrt(normInt.val(w1, w1).real());
460  complex<double> amp = decayamp / nrm * prodAmps[isample]->at(iw);
461  if(isample == 0) {// fill wheights of individual waves
462  weights[iw] = norm(amp);
463  }
464 
465  if (reflectivities[iw] == 0) {
466  if(ms[iw] == 0) {
467  posm0amps[ranks[iw]] += amp;
468  } else if(ms[iw] == 1) {
469  posm1amps[ranks[iw]] += amp;
470  }
471  } else if(reflectivities[iw] == 1) {
472  if(ms[iw] == 0) {
473  negm0amps[ranks[iw]] += amp;
474  } else if(ms[iw] == 1) {
475  negm1amps[ranks[iw]] += amp;
476  }
477  }
478  }
479  // end loop over waves
480 
481  // incoherent sum over rank and diffrerent reflectivities
482  weight = 0;
483  if(hasfit) {
484  for(int ir = 0; ir < maxrank + 1; ++ir) {
485  weight += norm(posm0amps[ir] + posm1amps[ir]);
486  //weight += norm(posm1amps[ir]);
487  weight += norm(negm0amps[ir] + negm1amps[ir]);
488  // weight += norm(negm1amps[ir]);
489  }
490  }
491  // add flat
492  weight += norm(flatamp);
493 
494 
495  if(isample == 0) {
496  hWeights->Fill(weight);
497  }
498  modelweights[isample] = weight;
499 
500  }// end loop over model samples
501  outtree->Fill();
502 
503  } // end loop over events
504  } // end we have an eventfile with decay amplitudes
505  else {
506  // Initialize PDGTable
507  rpwa::particleDataTable& pdt = rpwa::particleDataTable::instance();
508  pdt.readFile();
509 
510  TClonesArray prodKinPar("TObjString");
511  TClonesArray prodKinMom("TVector3");
512  TClonesArray decayKinPar("TObjString");
513  TClonesArray decayKinMom("TVector3");
514  TClonesArray *pprodKinPar = &prodKinPar;
515  TClonesArray *pprodKinMom = &prodKinMom;
516  TClonesArray *pdecayKinPar = &decayKinPar;
517  TClonesArray *pdecayKinMom = &decayKinMom;
518 
519  TFile* infile = TFile::Open(evtfilename.c_str(), "READ");
520  TTree *evt_tree;
521  infile->GetObject("rootPwaEvtTree;1", evt_tree);
522 
523  evt_tree->SetBranchAddress("prodKinParticles", &pprodKinPar);
524  evt_tree->SetBranchAddress("prodKinMomenta", &pprodKinMom);
525  evt_tree->SetBranchAddress("decayKinParticles", &pdecayKinPar);
526  evt_tree->SetBranchAddress("decayKinMomenta", &pdecayKinMom);
527 
528  unsigned int counter = 0;
529  unsigned long entries = evt_tree->GetEntries();
530  for(unsigned long i = 0; i < entries; i++) {
531  // evt2tree
532  evt_tree->GetEntry(i);
533  if(counter % 1000 == 0) {
534  std::cout << ".";
535  }
536  if(counter++ % 10000 == 0) {
537  std::cout << counter;
538  }
539  p->Delete(); // clear output arrays
540  q.clear();
541 
542  // we need 4 momentum and charge of beam and decay particles
543  // in rootPwaEvtTree the particle name plus the 3 momentum is stored
544  // so first get charge and particle mass from pdg table via particle name
545  // then construct the 4 momenta
546 
547 
548  // beam particle
549  rpwa::particle beampar(((TString*)prodKinPar[0])->Data()); // beam particle should be at first place -> 0
550  qbeam = beampar.charge();
551  beam.SetVectM(*(TVector3*)prodKinMom[0], pdt.entry(beampar.bareName())->mass());
552 
553  // decay particles
554  for(int dp = 0; dp < decayKinPar.GetEntries(); dp++) {
555  rpwa::particle dpar(((TString*)prodKinPar[dp])->Data());
556  q.push_back(dpar.charge());
557  new ((*p)[p->GetEntries()]) TLorentzVector(*(TVector3*)decayKinMom[dp],
558  pdt.entry(dpar.bareName())->mass());
559  }
560 
561  // weighting
562 
563  vector<complex<double> > posm0amps(maxrank + 1); // positive refl vector m=0
564  vector<complex<double> > posm1amps(maxrank + 1); // positive refl vector m=1
565 
566  vector<complex<double> > negm0amps(maxrank + 1); // negative refl vector m=0
567  vector<complex<double> > negm1amps(maxrank + 1); // negative refl vector m=1
568 
569  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
570  complex<double> decayamp;
571  ampfiles[iw]->read((char*) &decayamp, sizeof(complex<double> ));
572  string w1 = waveNames[iw];
573  //cerr << w1 << " " << decayamp << endl;
574  double nrm = sqrt(normInt.val(w1, w1).real());
575  complex<double> amp = decayamp / nrm * prodAmps[0]->at(iw);
576  if(reflectivities[iw] == 1) {
577  if(ms[iw] == 0) {
578  posm0amps[ranks[iw]] += amp;
579  } else if(ms[iw] == 1) {
580  posm1amps[ranks[iw]] += amp;
581  }
582  } else {
583  if(ms[iw] == 0) {
584  negm0amps[ranks[iw]] += amp;
585  } else if(ms[iw] == 1) {
586  negm1amps[ranks[iw]] += amp;
587  }
588  }
589  } // end loop over waves
590 
591  // incoherent sum:
592  weight = 0;
593  if(hasfit) {
594  for(int ir = 0; ir < maxrank + 1; ++ir) {
595  weight += norm(posm0amps[ir]+posm1amps[ir]);
596  //weight += norm(posm1amps[ir]);
597  weight += norm(negm0amps[ir]+negm1amps[ir]);
598  //weight += norm(negm1amps[ir]);
599  }
600  }
601  hWeights->Fill(weight);
602  outtree->Fill();
603  }
604  } // end of event loop
605  cout << endl << "Processed " << counter << " events" << endl;
606 
607  outfile->cd();
608  hWeights->Write();
609  outtree->Write();
610  outfile->Close();
611 
612  for(unsigned int iw = 0; iw < nmbWaves; ++iw) {
613  ampfiles[iw]->close();
614  delete ampfiles[iw];
615  }
616  ampfiles.clear();
617  for(unsigned int isamples=0; isamples < nsamples; ++isamples) {
618  delete prodAmps[isamples];
619  }
620  prodAmps.clear();
621 
622  return 0;
623 
624 }
625 
Generated by   doxygen 1.8.1.2