ROOTPWA
testNBodyPhaseSpaceGen.C
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 // tests nBodyPhaseSpaceGen
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <iostream>
39 #include <string>
40 #include <algorithm>
41 
42 #include "TStopwatch.h"
43 #include "TFile.h"
44 #include "TLorentzVector.h"
45 #include "TRandom3.h"
46 #include "TGenPhaseSpace.h"
47 #include "TH1.h"
48 #include "TH2.h"
49 #include "TCanvas.h"
50 #include "TF1.h"
51 #include "TF2.h"
52 #include "TColor.h"
53 
54 #include "reportingUtils.hpp"
55 #include "physUtils.hpp"
56 #include "conversionUtils.hpp"
57 #include "nBodyPhaseSpaceGen.h"
58 
59 
60 using namespace std;
61 using namespace rpwa;
62 
63 
64 const double gChargedPionMass = 0.13957018; // charged pion rest mass [GeV/c^2] [PDG08]
65 
66 
67 // constucts Lorentz vector of mother particle
68 TLorentzVector constructMother(TRandom3& random,
69  const double mass, // [GeV/c^2]
70  const double maxPt = 0.150, // [GeV/c]
71  const double maxEta = 1)
72 {
73  double pt = maxPt * random.Rndm(); // n-body transverse momentum in lab frame [0, 150] MeV/c
74  double eta = 2 * maxEta * random.Rndm() - maxEta; // n-body pseudorapidity in lab frame [-1, 1]
75  double phi = 4 * asin(1.) * random.Rndm(); // n-body azimuth in lab frame [0, 2 pi]
76  double tanhEta = tanh(eta);
77  TLorentzVector mother;
78  mother.SetPx(pt * cos(phi));
79  mother.SetPy(pt * sin(phi));
80  mother.SetPz(tanhEta * sqrt((mass * mass + pt * pt) / (1 - tanhEta * tanhEta)));
81  mother.SetVectM(mother.Vect(), mass);
82  return mother;
83 }
84 
85 
86 // signature for TF1
87 double
88 dalitzKinematicBorder(double* var,
89  double* par)
90 {
91  // variables
92  const double mass_2 = var[0]; // 2-body mass squared on x-axis
93 
94  // parameters
95  const double M = par[0]; // 3-body mass
96  const double* m = &par[1]; // array with the 3 daughter masses
97  const bool min = (par[4] < 0) ? true : false; // switches between curves for minimum and maximum mass squared on y-axis
98 
99  return dalitzKinematicBorder(mass_2, M, m, min);
100 }
101 
102 
103 void setResidualPalette(const unsigned int nmbSteps = 99)
104 {
105  const unsigned int nmbCol = 3;
106  Double_t rChannel[nmbCol] = {0, 1, 1};
107  Double_t gChannel[nmbCol] = {0, 1, 0};
108  Double_t bChannel[nmbCol] = {1, 1, 0};
109  Double_t zPos [nmbCol] = {0, 0.5, 1};
110  TColor::CreateGradientColorTable(nmbCol, zPos, rChannel, gChannel, bChannel, nmbSteps);
111 }
112 
113 
114 // fit function for flat Dalitz plot
115 double
116 dalitzFitFunc(double* var,
117  double* par)
118 {
119  const double m01_2 = var[0];
120  const double m12_2 = var[1];
121 
122  // parameters
123  const double M = par[1]; // 3-body mass
124  const double* m = &par[2]; // array with the 3 daughter masses
125 
126  if ( (m12_2 > dalitzKinematicBorder(m01_2, M, m, true ))
127  && (m12_2 < dalitzKinematicBorder(m01_2, M, m, false)))
128  return par[0];
129  else {
130  TF2::RejectPoint();
131  return 0;
132  }
133 }
134 
135 
136 // generalization of TF2::Integral
137 double
138 nDimIntegral(TF1* f,
139  int n, // number of variables
140  const double* min, // array of lower integration borders
141  const double* max, // array of upper integration borders
142  const double* funcPar, // function parameters
143  int maxNmbOfPoints = 10000, // maximum allowed number of function evaluations
144  double epsilon = 0.00001) // relative accuracy
145 {
146  if (n == 1)
147  return f->Integral(min[0], max[0], funcPar);
148  f->SetParameters(funcPar);
149  double resRelErr; // estimate of the relative accuracy of the result
150  int nmbOfFuncEval; // number of function evaluations performed
151  int errCode; // 0 = normal exit, 1 = maxNmbOfPoints too small, 3 = n < 2 or n > 15
152  double result = f->IntegralMultiple(n, min, max, -1, maxNmbOfPoints, epsilon, resRelErr, nmbOfFuncEval, errCode);
153  return result;
154 }
155 
156 
157 double
158 nBodyPhaseSpaceElement(double* var, // effective masses of the intermediate (i + 2)-body systems
159  double* par) // # of daughters, n daughter masses, M_n
160 {
161  const int n = static_cast<int>(par[0]); // number of decay daughters > 2
162  const double* m = &par[1]; // masses of the daughter particles
163  double M[n]; // (i + 1)-body invariant mass
164  M[0] = m[0]; // 1-body mass is mass of first daughter
165  for (int i = 1; i < (n - 1); ++i) // intermediate (i + 1)-body effective masses
166  M[i] = var[i - 1];
167  M[n - 1] = par[n + 1]; // n-body mass
168 
169  // reject points outside of kinematical range
170  for (int i = 1; i < n; ++i)
171  if (M[i] < M[i - 1] + m[i]) {
172  TF1::RejectPoint();
173  return 0;
174  }
175 
176  double momProd = 1; // product of breakup momenta
177  for (int i = 1; i < n; ++i) // loop over 2- to n-bodies
178  momProd *= breakupMomentum(M[i], M[i - 1], m[i]);
179  return momProd; // the term that needs to be intergrated; normalization in integration routine
180 }
181 
182 
183 TF1* dLips; // n-body Lorentz-invariant phase space element
184 
185 
186 double
187 nBodyPhaseSpace(double* var, // n-body mass
188  double* par) // # of daughters, n daughter masses
189 {
190  const double M_n = var[0]; // n-body mass
191 
192  const int n = static_cast<int>(par[0]); // number of decay daughters > 2
193  const double* m = &par[1]; // masses of the daughter particles
194  const double A = par[n + 1]; // scale factor
195 
196  // reject points outside of kinematical range
197  double mSum = 0;
198  for (int i = 0; i < n; ++i)
199  mSum += m[i];
200  if (M_n < mSum)
201  return 0;
202 
203  // normalization
204  const double f_n = 1 / (2 * pow(twoPi, 2 * n - 3)); // S. U. Chung's normalization
205  //const double f_n = pi * pow(twoPi, n - 2); // genbod's normalization
206  //const double f_n = 1;
207  const double norm = f_n / M_n;
208 
209  if (n == 2)
210  return A * norm * breakupMomentum(M_n, m[0], m[1]);
211 
212  // prepare parameter array for phase space element
213  double lipsPar[n + 2];
214  lipsPar[0] = n;
215  for (int i = 1; i <= n; ++i)
216  lipsPar[i] = m[i - 1];
217  lipsPar[n + 1] = M_n;
218 
219  // set integration limits
220  double min[n - 2]; // lower integration limits for effective (i + 2)-body masses
221  min[0] = m[0] + m[1];
222  for (int i = 1; i < n - 2; ++i)
223  min[i] = min[i - 1] + m[i + 1];
224  double max[n - 2]; // upper integration limits for effective (i + 2)-body masses
225  mSum = 0;
226  for (int i = n - 3; i >= 0; --i) {
227  mSum += m[i + 2];
228  max[i] = M_n - mSum;
229  }
230 
231  // perform integration
232  return A * norm * nDimIntegral(dLips, n - 2, min, max, lipsPar, 500000, 0.00001);
233 }
234 
235 
236 void testNBodyPhaseSpaceGen(const unsigned int nmbEvents = 1000000,
237  const string& outFileName = "./testNBodyPhaseSpaceGen.root")
238 {
239  TStopwatch timer;
240  timer.Start();
241 
242  // set parameters
243  const unsigned int n = 3; // number of daughters
244  double m[n] = {gChargedPionMass, gChargedPionMass, gChargedPionMass}; // daughter masses
245  const double massRange[2] = {m[0] + m[1] + m[2], 2}; // [GeV/c^2]
246  const double mass = (massRange[0] + massRange[1]) / 2;
247  const unsigned int seed = 1234567890;
248  const bool hitMissMode = false;
249 
250  // create output file
251  printInfo << "creating output file '" << outFileName << "' ... " << endl;
252  TFile* outFile = TFile::Open(outFileName.c_str(), "RECREATE");
253  if (!outFile) {
254  printErr << "could not create output file. aborting..." << endl;
255  return;
256  }
257  cout << " ... success" << endl;
258 
259  // setup nBodyPhaseSpaceGen generators
260  const unsigned int nmbGen = 5;
261  nBodyPhaseSpaceGen* gen[5];
262  for (unsigned int i = 0; i < nmbGen; ++i)
263  gen[i] = new nBodyPhaseSpaceGen();
264  gen[0]->setWeightType(nBodyPhaseSpaceGen::S_U_CHUNG);
265  gen[1]->setWeightType(nBodyPhaseSpaceGen::NUPHAZ);
266  gen[2]->setWeightType(nBodyPhaseSpaceGen::GENBOD);
267  gen[3]->setWeightType(nBodyPhaseSpaceGen::FLAT);
268  gen[4]->setWeightType(nBodyPhaseSpaceGen::S_U_CHUNG);
269  for (unsigned int i = 0; i < 4; ++i)
270  gen[i]->setKinematicsType(nBodyPhaseSpaceGen::RAUBOLD_LYNCH);
271  gen[4]->setKinematicsType(nBodyPhaseSpaceGen::BLOCK);
272  for (unsigned int i = 0; i < nmbGen; ++i) {
273  gen[i]->setDecay(n, m);
274  gen[i]->setSeed(seed);
275  gen[i]->setMaxWeight(1.01 * gen[i]->estimateMaxWeight(mass));
276  }
277  // setup TGenPhaseSpace generator
278  TGenPhaseSpace TGen;
279  gRandom->SetSeed(seed); // TGenPhaseSpace uses gRandom
280  double maxTGenWeight = 0;
281  for (unsigned int i = 0; i < 10000; ++i) {
282  TLorentzVector mother(0, 0, 0, mass);
283  TGen.SetDecay(mother, (int)n, m);
284  const double weight = TGen.Generate();
285  maxTGenWeight = (weight > maxTGenWeight) ? weight : maxTGenWeight;
286  }
287  maxTGenWeight *= 1.01;
288  double maxTGenWeightObserved = 0;
289  TRandom3 random(seed);
290 
291  // book histograms
292  TH2F* hDalitz [nmbGen + 1];
293  TH1F* hMass [nmbGen + 1];
294  for (unsigned int i = 0; i < nmbGen; ++i) {
295  hDalitz[i] = new TH2F(("hDalitz" + toString(i)).c_str(),
296  ("Dalitz Plot Residuum nBodyPhaseSpaceGen " + toString(i) + ";m_{01}^{2} [(GeV/c^{2})^{2}];m_{12}^{2} [(GeV/c^{2})^{2}]").c_str(),
297  100, 0, 1.5, 100, 0, 1.5);
298  hMass[i] = new TH1F(("hMass" + toString(i)).c_str(),
299  ("Mother Mass nBodyPhaseSpaceGen " + toString(i) + ";m [GeV/c^{2}]").c_str(),
300  1000, 0, 2.5);
301  }
302  hDalitz [nmbGen] = new TH2F("hDalitzTGen", "Dalitz Plot Residuum TGen;m_{01}^{2} [(GeV/c^{2})^{2}];m_{12}^{2} [(GeV/c^{2})^{2}]", 100, 0, 1.5, 100, 0, 1.5);
303  hMass [nmbGen] = new TH1F("hMassTGen", "Mother Mass TGen;m [GeV/c^{2}]", 1000, 0, 2.5);
304  TH1F* hMassRaw = new TH1F("hMassRaw", "Mother Mass Raw;m [GeV/c^{2}]", 1000, 0, 2.5);
305 
306  if (0) {
307  printInfo << "testing generators for constant 3-body mass = " << mass << " GeV/c^2" << endl
308  << " generating " << nmbEvents << " events ... " << endl;
309  unsigned int countAtt[nmbGen + 1];
310  unsigned int countAcc[nmbGen + 1];
311  bool done = false;
312  for (unsigned int i = 0; i <= nmbGen; ++i)
313  countAtt[i] = countAcc[i] = 0;
314  while (!done) {
315  TLorentzVector mother = constructMother(random, mass);
316  for (unsigned int i = 0; i < nmbGen; ++i)
317  if (countAcc[i] < nmbEvents){
318  ++countAtt[i];
319  if (gen[i]->generateDecayAccepted(mother)) {
320  ++countAcc[i];
321  const TLorentzVector m01 = gen[i]->daughter(0) + gen[i]->daughter(1);
322  const TLorentzVector m12 = gen[i]->daughter(1) + gen[i]->daughter(2);
323  hDalitz[i]->Fill(m01.M2(), m12.M2());
324  }
325  }
326  TGen.SetDecay(mother, (int)n, m);
327  const double TGenWeight = TGen.Generate();
328  maxTGenWeightObserved = (TGenWeight > maxTGenWeightObserved) ? TGenWeight : maxTGenWeightObserved;
329  if (countAcc[nmbGen] < nmbEvents) {
330  ++countAtt[nmbGen];
331  if ((TGenWeight / maxTGenWeight) > random.Rndm()) {
332  ++countAcc[nmbGen];
333  const TLorentzVector m01 = *TGen.GetDecay(0) + *TGen.GetDecay(1);
334  const TLorentzVector m12 = *TGen.GetDecay(1) + *TGen.GetDecay(2);
335  hDalitz[nmbGen]->Fill(m01.M2(), m12.M2());
336  }
337  }
338  done = true;
339  for (unsigned int i = 0; i <= nmbGen; ++i)
340  if (countAcc[i] < nmbEvents)
341  done = false;
342  }
343  cout << " ... done." << endl << endl;
344  for (unsigned int i = 0; i < nmbGen; ++i)
345  printInfo << "gen[" << i << "]: " << countAtt[i] << " attempts, " << ((double)nmbEvents / countAtt[i]) * 100 << " % efficiency" << endl
346  << *gen[i] << endl;
347  printInfo << "TGenPhaseSpace: " << countAtt[nmbGen] << " attempts, " << ((double)nmbEvents / countAtt[nmbGen]) * 100 << " % efficiency" << endl
348  << "TGenPhaseSpace maximum weight = " << maxTGenWeight
349  << " vs. maximum weight observed = " << maxTGenWeightObserved << endl << endl;
350 
351  // plot residuum of Dalitz plots after subtracting expected flat component
352  setResidualPalette(99);
353  for (unsigned int i = 0; i <= nmbGen; ++i) {
354  new TCanvas();
355  const double xMin = (m[0] + m[1]) * (m[0] + m[1]);
356  const double xMax = (mass - m[2]) * (mass - m[2]);
357  const double yMin = (m[1] + m[2]) * (m[1] + m[2]);
358  const double yMax = (mass - m[0]) * (mass - m[0]);
359  TF2* dalitzFit = new TF2("dalitzFit", dalitzFitFunc, xMin, xMax, yMin, yMax, 5);
360  dalitzFit->SetParameter(0, 1);
361  dalitzFit->FixParameter(1, mass);
362  dalitzFit->FixParameter(2, m[0]);
363  dalitzFit->FixParameter(3, m[1]);
364  dalitzFit->FixParameter(4, m[2]);
365  const double dalitzInt = dalitzFit->Integral(xMin, xMax, yMin, yMax, 1e-6)
366  / (hDalitz[i]->GetXaxis()->GetBinWidth(1) * hDalitz[i]->GetYaxis()->GetBinWidth(1));
367  dalitzFit->SetParameter(0, nmbEvents / dalitzInt);
368  hDalitz[i]->Sumw2();
369  hDalitz[i]->Add(dalitzFit, -1, "I"); // ROOT does not support I flag for 2D functions; increased residual at borders
370  // center z-range
371  const double maxVal = max(fabs(hDalitz[i]->GetMinimum()), fabs(hDalitz[i]->GetMaximum()));
372  hDalitz[i]->SetMinimum(-maxVal);
373  hDalitz[i]->SetMaximum( maxVal);
374  hDalitz[i]->SetContour(99);
375  hDalitz[i]->Draw("COLZ");
376  // draw border
377  TF1* dalitzBorder = new TF1("dalitzBorder", dalitzKinematicBorder, xMin, xMax, 5);
378  dalitzBorder->SetParameter(0, mass);
379  dalitzBorder->SetParameter(1, m[0]);
380  dalitzBorder->SetParameter(2, m[1]);
381  dalitzBorder->SetParameter(3, m[2]);
382  dalitzBorder->SetParameter(4, -1);
383  dalitzBorder->SetNpx(1000);
384  dalitzBorder->SetLineWidth(1);
385  dalitzBorder->SetLineColor(3);
386  dalitzBorder->DrawCopy("SAME");
387  dalitzBorder->SetParameter(4, +1);
388  dalitzBorder->DrawCopy("SAME");
389  gPad->Update();
390  gPad->RedrawAxis();
391  }
392  }
393 
394  if (1) {
395  printInfo << "testing generators for 3-body in mass range [" << massRange[0] << ", " << massRange[1] << "] GeV/c^2" << endl
396  << " generating " << nmbEvents << " events ... " << endl;
397  // recalculate maximum weights
398  gen[0]->setMaxWeight(1.01 * gen[0]->estimateMaxWeight(massRange[1]));
399  gen[1]->setMaxWeight(1.01 * gen[1]->estimateMaxWeight(massRange[1]));
400  gen[2]->setMaxWeight(1.01 * gen[2]->estimateMaxWeight(massRange[0] + 0.0001));
401  gen[3]->setMaxWeight(1.01 * gen[3]->estimateMaxWeight(massRange[1]));
402  gen[4]->setMaxWeight(1.01 * gen[4]->estimateMaxWeight(massRange[1]));
403  maxTGenWeight = 0;
404  for (unsigned int i = 0; i < 10000; ++i) {
405  TLorentzVector mother(0, 0, 0, massRange[0] + 0.0001);
406  TGen.SetDecay(mother, (int)n, m);
407  const double weight = TGen.Generate();
408  maxTGenWeight = (weight > maxTGenWeight) ? weight : maxTGenWeight;
409  }
410  maxTGenWeight *= 1.01;
411  maxTGenWeightObserved = 0;
412 
413  unsigned int countAtt[nmbGen + 1];
414  unsigned int countAcc[nmbGen + 1];
415  bool done = false;
416  for (unsigned int i = 0; i <= nmbGen; ++i)
417  countAtt[i] = countAcc[i] = 0;
418  while (!done) {
419  TLorentzVector mother = constructMother(random, massRange[0] + random.Rndm() * (massRange[1] - massRange[0]));
420  hMassRaw->Fill(mother.M());
421  for (unsigned int i = 0; i < nmbGen; ++i)
422  if (countAcc[i] < nmbEvents) {
423  ++countAtt[i];
424  if (hitMissMode) {
425  if (gen[i]->generateDecayAccepted(mother)) {
426  ++countAcc[i];
427  hMass[i]->Fill(mother.M());
428  }
429  } else {
430  hMass[i]->Fill(mother.M(), gen[i]->generateDecay(mother));
431  if (gen[i]->eventAccepted())
432  ++countAcc[i];
433  }
434  }
435  TGen.SetDecay(mother, (int)n, m);
436  const double TGenWeight = TGen.Generate();
437  maxTGenWeightObserved = (TGenWeight > maxTGenWeightObserved) ? TGenWeight : maxTGenWeightObserved;
438  if (countAcc[nmbGen] < nmbEvents) {
439  ++countAtt[nmbGen];
440  if (hitMissMode) {
441  if ((TGenWeight / maxTGenWeight) > random.Rndm())
442  hMass[nmbGen]->Fill(mother.M());
443  } else
444  hMass[nmbGen]->Fill(mother.M(), TGenWeight);
445  if ((TGenWeight / maxTGenWeight) > random.Rndm())
446  ++countAcc[nmbGen];
447  }
448  done = true;
449  for (unsigned int i = 0; i <= nmbGen; ++i)
450  if (countAcc[i] < nmbEvents)
451  done = false;
452  }
453  cout << " ... done." << endl << endl;
454  for (unsigned int i = 0; i < nmbGen; ++i)
455  printInfo << "gen[" << i << "]: " << countAtt[i] << " attempts, " << ((double)nmbEvents / countAtt[i]) * 100 << " % efficiency" << endl
456  << *gen[i] << endl;
457  printInfo << "TGenPhaseSpace: " << countAtt[nmbGen] << " attempts, " << ((double)nmbEvents / countAtt[nmbGen]) * 100 << " % efficiency" << endl
458  << "maximum weight = " << maxTGenWeight
459  << " vs. maximum weight observed = " << maxTGenWeightObserved << endl << endl;
460 
461  // draw histograms
462  dLips = new TF1("nBodyPhaseSpaceElement", nBodyPhaseSpaceElement, 0, 100, n + 2);
463  TF1* Lips = new TF1("nBodyPhaseSpace", nBodyPhaseSpace, massRange[0], massRange[1], n + 2);
464  Lips->SetParameter(0, n);
465  for (unsigned int i = 1; i <= n; ++i)
466  Lips->SetParameter(i, m[i - 1]);
467  Lips->SetLineColor(2);
468  Lips->SetLineWidth(2);
469  for (unsigned int i = 0; i <= nmbGen; ++i) {
470  Lips->SetParameter(n + 1, 1);
471  double A;
472  if (hitMissMode)
473  // scale function to last bin
474  A = hMass[i]->GetBinContent(hMass[i]->FindBin(massRange[1]) - 1) / Lips->Eval(massRange[1]);
475  else {
476  // scale histograms by number of events per bin
477  for (int j = 1; j < hMassRaw->GetNbinsX(); ++j) {
478  const double nmbEvBin = hMassRaw->GetBinContent(j);
479  if (nmbEvBin > 0)
480  hMass[i]->SetBinContent(j, hMass[i]->GetBinContent(j) / nmbEvBin);
481  }
482  A = 1;
483  }
484  Lips->SetParameter(n + 1, A);
485  new TCanvas();
486  hMass[i]->Draw();
487  Lips->DrawCopy("LSAME");
488  if (hitMissMode)
489  hMass[i]->Draw("SAME E");
490  }
491  new TCanvas;
492  hMassRaw->Draw();
493  }
494 
495  if (0) {
496  const unsigned int n2 = 4; // number of daughters
497  double m2[n2] = {gChargedPionMass, gChargedPionMass,
498  gChargedPionMass, gChargedPionMass}; // daughter masses
499  const double massRange2[2] = {m2[0] + m2[1] + m2[2] + m2[3], 2}; // [GeV/c^2]
500  printInfo << "comparing run time of generators for 4-body in mass range [" << massRange2[0] << ", " << massRange2[1] << "] GeV/c^2" << endl
501  << " generating " << nmbEvents << " events using TGenPhaseSpace ... " << endl;
502 
503  // setup nBodyPhaseSpaceGen generator
504  nBodyPhaseSpaceGen gen2;
505  gen2.setWeightType(nBodyPhaseSpaceGen::GENBOD);
506  //gen2.setWeightType(nBodyPhaseSpaceGen::S_U_CHUNG);
507  gen2.setKinematicsType(nBodyPhaseSpaceGen::BLOCK);
508  gen2.setDecay(n2, m2);
509  gen2.setSeed(seed);
510  gen2.setMaxWeight(1.01 * gen2.estimateMaxWeight(massRange2[0] + 0.0001));
511  //gen2.setMaxWeight(1.01 * gen2.estimateMaxWeight(massRange2[1]));
512  // recalculate maximum weight for TGenPhaseSpace generator
513  maxTGenWeight = 0;
514  for (unsigned int i = 0; i < 10000; ++i) {
515  TLorentzVector mother(0, 0, 0, massRange2[0] + 0.0001);
516  TGen.SetDecay(mother, (int)n2, m2);
517  const double weight = TGen.Generate();
518  maxTGenWeight = (weight > maxTGenWeight) ? weight : maxTGenWeight;
519  }
520  maxTGenWeight *= 1.01;
521  maxTGenWeightObserved = 0;
522 
523  TStopwatch timer2;
524  timer2.Start();
525  unsigned int countAtt = 0;
526  unsigned int countAcc = 0;
527  random.SetSeed(seed);
528  while (countAcc < nmbEvents) {
529  TLorentzVector mother = constructMother(random, massRange2[0] + random.Rndm() * (massRange2[1] - massRange2[0]));
530  TGen.SetDecay(mother, (int)n2, m2);
531  const double TGenWeight = TGen.Generate();
532  maxTGenWeightObserved = (TGenWeight > maxTGenWeightObserved) ? TGenWeight : maxTGenWeightObserved;
533  ++countAtt;
534  if ((TGenWeight / maxTGenWeight) > random.Rndm())
535  ++countAcc;
536  }
537  timer2.Stop();
538  const double timeTGen = timer2.RealTime();
539  cout << " ... done." << endl
540  << " " << countAtt << " attempts, " << ((double)nmbEvents / countAtt) * 100 << " % efficiency" << endl
541  << " consumed time: ";
542  timer2.Print();
543  timer2.Reset();
544  cout << " TGenPhaseSpace maximum weight = " << maxTGenWeight
545  << " vs. maximum weight observed = " << maxTGenWeightObserved << endl;
546 
547  cout << endl << " generating " << nmbEvents << " events using nBodyPhaseSpaceGen ... " << endl;
548  timer2.Start();
549  countAtt = 0;
550  countAcc = 0;
551  random.SetSeed(seed);
552  while (countAcc < nmbEvents) {
553  TLorentzVector mother = constructMother(random, massRange2[0] + random.Rndm() * (massRange2[1] - massRange2[0]));
554  ++countAtt;
555  if (gen2.generateDecayAccepted(mother))
556  ++countAcc;
557  }
558  timer2.Stop();
559  const double timeNBody = timer2.RealTime();
560  cout << " ... done." << endl
561  << " " << countAtt << " attempts, " << ((double)nmbEvents / countAtt) * 100 << " % efficiency" << endl
562  << " consumed time: ";
563  timer2.Print();
564  timer2.Reset();
565  printInfo << gen2 << endl;
566  printInfo << "nBodyPhaseSpaceGen speed gain w.r.t. TGenPhaseSpace = " << 100 * timeTGen / timeNBody << " %" << endl << endl;
567 
568  printInfo << "checking mass dependence of nBodyPhaseSpaceGen 4-body phase space in mass range [" << massRange2[0] << ", " << massRange2[1] << "] GeV/c^2" << endl;
569  gen2.setWeightType(nBodyPhaseSpaceGen::S_U_CHUNG);
570  gen2.setKinematicsType(nBodyPhaseSpaceGen::BLOCK);
571  gen2.setDecay(n2, m2);
572  gen2.setSeed(seed);
573  gen2.setMaxWeight(1.01 * gen2.estimateMaxWeight(massRange2[1]));
574  TH1F* hMass2 = new TH1F("hMass4Body", "4-Body Mass nBodyPhaseSpaceGen;m [GeV/c^{2}]", 1000, 0, 2.5);
575  countAtt = 0;
576  countAcc = 0;
577  random.SetSeed(seed);
578  while (countAcc < nmbEvents) {
579  TLorentzVector mother = constructMother(random, massRange2[0] + random.Rndm() * (massRange2[1] - massRange2[0]));
580  ++countAtt;
581  if (gen2.generateDecayAccepted(mother)) {
582  ++countAcc;
583  hMass2->Fill(mother.M());
584  }
585  }
586  cout << " ... done." << endl
587  << " " << countAtt << " attempts, " << ((double)nmbEvents / countAtt) * 100 << " % efficiency" << endl
588  << gen2 << endl;
589  dLips = new TF1("nBodyPhaseSpaceElement", nBodyPhaseSpaceElement, 0, 100, n2 + 2);
590  TF1* Lips = new TF1("nBodyPhaseSpace", nBodyPhaseSpace, massRange2[0], massRange2[1], n2 + 2);
591  Lips->SetParameter(0, n2);
592  for (unsigned int i = 1; i <= n2; ++i)
593  Lips->SetParameter(i, m2[i - 1]);
594  Lips->SetLineColor(2);
595  Lips->SetLineWidth(2);
596  Lips->SetParameter(n2 + 1, 1);
597  const double A = hMass2->GetBinContent(hMass2->FindBin(massRange2[1]) - 1) / Lips->Eval(massRange2[1]);
598  Lips->SetParameter(n2 + 1, A);
599  new TCanvas();
600  hMass2->Draw();
601  Lips->DrawCopy("LSAME");
602  hMass2->Draw("SAME E");
603  }
604 
605  timer.Stop();
606  printInfo << endl << "this job consumed: ";
607  timer.Print();
608 }
Generated by   doxygen 1.8.1.2