ROOTPWA
testMassDependence.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 // basic test program for mass dependence functions
29 //
30 //
31 // Author List:
32 // Boris Grube TUM (original author)
33 //
34 //
35 //-------------------------------------------------------------------------
36 
37 
38 #include <boost/timer/timer.hpp>
39 
40 #include "TFile.h"
41 #include "TH1.h"
42 #include "TH2.h"
43 #include "TMath.h"
44 
45 #include "particleDataTable.h"
46 #include "massDependence.h"
47 #include "isobarDecayVertex.h"
48 #include "nBodyPhaseSpaceGen.h"
49 #include "reportingUtilsRoot.hpp"
50 
51 #include "massDep.h"
52 #include "particle.h"
53 
54 
55 extern particleDataTable PDGtable;
56 
57 
58 using namespace std;
59 using namespace boost;
60 using namespace boost::timer;
61 using namespace rpwa;
62 
63 
64 // constucts Lorentz vector of parent particle
65 TLorentzVector
66 constructParent(TRandom3& random,
67  const double mass, // [GeV/c^2]
68  const double maxPt = 0.150, // [GeV/c]
69  const double maxEta = 1)
70 {
71  double pt = maxPt * random.Rndm(); // parent transverse momentum in lab frame
72  double eta = 2 * maxEta * random.Rndm() - maxEta; // parent pseudorapidity in lab frame
73  double phi = 2 * pi * random.Rndm(); // parent azimuth in lab frame
74  double tanhEta = tanh(eta);
75  TLorentzVector parent;
76  parent.SetPx(pt * cos(phi));
77  parent.SetPy(pt * sin(phi));
78  parent.SetPz(tanhEta * sqrt((mass * mass + pt * pt) / (1 - tanhEta * tanhEta)));
79  parent.SetVectM(parent.Vect(), mass);
80  return parent;
81 }
82 
83 
84 void
85 compareAmplitudes(const size_t nmbEvents,
86  const unsigned int seed,
87  massDependencePtr& massDep,
88  isobarDecayVertexPtr& vertex,
89  ::massDep* pwa2kMassDep,
90  ::particle& pwa2kParent,
91  const string& name = "",
92  const double maxMass = 3)
93 {
94  vertex->setMassDependence(massDep);
95  pwa2kParent.setMassDep(pwa2kMassDep);
96  printInfo << "comparing '" << name << "' dynamic amplitudes" << endl
97  << " " << *(vertex->massDependence()) << endl
98  << " ";
99  pwa2kMassDep->print();
100  cout << endl
101  << " using decays" << endl
102  << " " << *vertex << endl;
103  pwa2kParent.print();
104 
105  // setup phase-space generator
106  particlePtr daughter1 = vertex->daughter1();
107  particlePtr daughter2 = vertex->daughter2();
108  const double daughterMasses[2] = {daughter1->mass(), daughter2->mass()};
109  const double massRange [2] = {daughterMasses[0] + daughterMasses[1], maxMass}; // [GeV/c^2]
110  nBodyPhaseSpaceGen psGen;
111  psGen.setWeightType(nBodyPhaseSpaceGen::S_U_CHUNG);
112  psGen.setKinematicsType(nBodyPhaseSpaceGen::BLOCK);
113  psGen.setDecay(2, daughterMasses);
114  psGen.setSeed(seed);
115  psGen.setMaxWeight(1.01 * psGen.estimateMaxWeight(massRange[1]));
116 
117  // loop over events
118  TRandom3 random(seed);
119  cpu_timer timer, timerPwa2k;
120  vector<complex<double> > newAmps (nmbEvents, 0);
121  vector<complex<double> > oldAmps (nmbEvents, 0);
122  vector<double> parentMass(nmbEvents, 0);
123  size_t countEvts = 0;
124  bool firstCall = true;
125  while (countEvts < nmbEvents) {
126  const TLorentzVector parentLv
127  = constructParent(random, massRange[0] + random.Rndm() * (massRange[1] - massRange[0]));
128  if (not psGen.generateDecayAccepted(parentLv))
129  continue;
130  // calculate new amplitude
131  daughter1->setMomentum(psGen.daughter(0).Vect());
132  daughter2->setMomentum(psGen.daughter(1).Vect());
133  vertex->calcParentLzVec();
134  if (firstCall)
135  timer.start();
136  else
137  timer.resume();
138  newAmps[countEvts] = vertex->massDepAmplitude();
139  timer.stop();
140  parentMass[countEvts] = parentLv.M();
141  // calculate PWA2000 amplitude
142  assert(pwa2kParent.Decay()->_children.size() == 2);
143  list< ::particle>::iterator child = pwa2kParent.Decay()->_children.begin();
144  ::particle& d1 = *child;
145  ++child;
146  ::particle& d2 = *child;
147  fourVec pwa2kMom;
148  pwa2kMom.set(daughter1->lzVec().E(), daughter1->lzVec().X(),
149  daughter1->lzVec().Y(), daughter1->lzVec().Z());
150  d1.set4P(pwa2kMom);
151  pwa2kMom.set(daughter2->lzVec().E(), daughter2->lzVec().X(),
152  daughter2->lzVec().Y(), daughter2->lzVec().Z());
153  d2.set4P(pwa2kMom);
154  pwa2kMom.set(vertex->parent()->lzVec().E(), vertex->parent()->lzVec().X(),
155  vertex->parent()->lzVec().Y(), vertex->parent()->lzVec().Z());
156  pwa2kParent.set4P(pwa2kMom);
157  if (firstCall)
158  timerPwa2k.start();
159  else
160  timerPwa2k.resume();
161  oldAmps[countEvts] = pwa2kMassDep->val(pwa2kParent);
162  timerPwa2k.stop();
163  ++countEvts;
164  if (firstCall)
165  firstCall = false;
166  }
167  printInfo << "calculated mass dependence for " << nmbEvents << " events" << endl
168  << " this consumed: (new code) "
169  << timer.format(3, "[%u sec user] + [%s sec sys] = [%t sec]")
170  << " vs. (PWA2000) " << timerPwa2k.format(3, "[%t sec] = [%u sec user] + [%s sec sys]")
171  << endl;
172 
173  // compare amplitude values
174  TH2D* hAmp = new TH2D(("h" + name + "Amp" ).c_str(), (name + " Amplitude" ).c_str(), 1000, -1, 1, 1000, -0.5, 1.5);
175  TH1D* hAmpReal = new TH1D(("h" + name + "AmpReal" ).c_str(), (name + " Amplitude (Real Part)" ).c_str(), 1000, -1.2, 1.2);
176  TH1D* hAmpImag = new TH1D(("h" + name + "AmpImag" ).c_str(), (name + " Amplitude (Imag Part)" ).c_str(), 1000, -1.2, 1.2);
177  TH2D* hInt = new TH2D(("h" + name + "Int" ).c_str(), (name + " Intensity" ).c_str(), 1000, massRange[0] - 0.2, massRange[1] + 0.2, 1000, 0, 1.2);
178  TH2D* hPhase = new TH2D(("h" + name + "Phase" ).c_str(), (name + " Phase" ).c_str(), 1000, massRange[0] - 0.2, massRange[1] + 0.2, 1000, -200, 200);
179  TH1D* hAbsDiffReal = new TH1D(("h" + name + "AbsDiffReal").c_str(), (name + " Absolute Diff. (Real Part)").c_str(), 100000, -1e-12, 1e-12);
180  TH1D* hAbsDiffImag = new TH1D(("h" + name + "AbsDiffImag").c_str(), (name + " Absolute Diff. (Imag Part)").c_str(), 100000, -1e-12, 1e-12);
181  TH1D* hRelDiffReal = new TH1D(("h" + name + "RelDiffReal").c_str(), (name + " Relative Diff. (Real Part)").c_str(), 100000, -1e-12, 1e-12);
182  TH1D* hRelDiffImag = new TH1D(("h" + name + "RelDiffImag").c_str(), (name + " Relative Diff. (Imag Part)").c_str(), 100000, -1e-12, 1e-12);
183  double maxAbsDiffReal = 0;
184  double maxRelDiffReal = 0;
185  double maxAbsDiffImag = 0;
186  double maxRelDiffImag = 0;
187  for (unsigned int i = 0; i < newAmps.size(); ++i) {
188  const complex<double> absDiff = oldAmps[i] - newAmps[i];
189  const complex<double> relDiff = complex<double>(absDiff.real() / oldAmps[i].real(),
190  absDiff.imag() / oldAmps[i].imag());
191  hAmpReal->Fill(newAmps[i].real());
192  hAmpImag->Fill(newAmps[i].imag());
193  hAmp->Fill(newAmps[i].real(), newAmps[i].imag());
194  hInt->Fill(parentMass[i], abs(newAmps[i]));
195  hPhase->Fill(parentMass[i], TMath::RadToDeg() * std::arg(newAmps[i]));
196  hAbsDiffReal->Fill(absDiff.real());
197  hAbsDiffImag->Fill(absDiff.imag());
198  hRelDiffReal->Fill(relDiff.real());
199  hRelDiffImag->Fill(relDiff.imag());
200  if (abs(absDiff.real()) > maxAbsDiffReal)
201  maxAbsDiffReal = abs(absDiff.real());
202  if (abs(absDiff.imag()) > maxAbsDiffImag)
203  maxAbsDiffImag = abs(absDiff.imag());
204  if ((oldAmps[i].real() != 0) and (abs(relDiff.real()) > maxRelDiffReal))
205  maxRelDiffReal = abs(relDiff.real());
206  if ((oldAmps[i].imag() != 0) and (abs(relDiff.imag()) > maxRelDiffImag))
207  maxRelDiffImag = abs(relDiff.imag());
208  }
209  printSucc << endl
210  << " real part: maximum absolute deviation is " << maxAbsDiffReal << "; "
211  << "maximum relative deviation is " << maxRelDiffReal << endl
212  << " imag part: maximum absolute deviation is " << maxAbsDiffImag << "; "
213  << "maximum relative deviation is " << maxRelDiffImag << endl;
214 }
215 
216 
217 int
218 main(int argc, char** argv)
219 {
220  printCompilerInfo();
221  printSvnVersion();
222  cout << endl;
223 
224  // define parameters
225  const size_t nmbEvents = 1000000;
226  //const size_t nmbEvents = 10;
227  const unsigned int seed = 123456789;
228 
229  rpwa::particleDataTable& pdt = rpwa::particleDataTable::instance();
230  pdt.readFile();
231  PDGtable.initialize("../../keyfiles/key5pi/pdgTable.txt");
232 
233  // massDependence::setDebug(true);
234  // ::particle::debug();
235 
236  // define decay
237  // particlePtr daughter1 = createParticle("pi+");
238  // ::particle pwa2kDaughter1(PDGtable.get("pi"), +1);
239  // particlePtr daughter2 = createParticle("pi-");
240  // ::particle pwa2kDaughter2(PDGtable.get("pi"), -1);
241  // const unsigned int L = 2;
242  // const unsigned int S = 0;
243  // particlePtr parent = createParticle("rho(770)0");
244  // ::particle pwa2kParent(PDGtable.get("rho(770)"), 0);
245  // const unsigned int L = 0;
246  // const unsigned int S = 0;
247  // particlePtr parent = createParticle("sigma");
248  // ::particle pwa2kParent(PDGtable.get("sigma"), 0);
249 
250  // particlePtr daughter1 = createParticle("rho(770)0");
251  // ::particle pwa2kDaughter1(PDGtable.get("rho(770)"), 0);
252  // particlePtr daughter2 = createParticle("sigma0");
253  // ::particle pwa2kDaughter2(PDGtable.get("sigma"), 0);
254  particlePtr daughter1 = createParticle("a1(1260)+");
255  ::particle pwa2kDaughter1(PDGtable.get("a1(1260)"), +1);
256  particlePtr daughter2 = createParticle("pi-");
257  ::particle pwa2kDaughter2(PDGtable.get("pi"), -1);
258  const unsigned int L = 0;
259  const unsigned int S = 2;
260  particlePtr parent = createParticle("rho(1700)0");
261  ::particle pwa2kParent(PDGtable.get("rho(1700)"), 0);
262 
263  // setup data structures
264  isobarDecayVertexPtr vertex = createIsobarDecayVertex(parent, daughter1, daughter2, L, S);
265  ::decay pwa2kDecay;
266  pwa2kDecay.addChild(pwa2kDaughter1);
267  pwa2kDecay.addChild(pwa2kDaughter2);
268  pwa2kDecay.setL(L);
269  pwa2kDecay.setS(S);
270  pwa2kParent.setDecay(pwa2kDecay);
271 
272  // compare amplitudes
273  TFile* outFile = TFile::Open("testMassDependence.root", "RECREATE");
274  const size_t nmbMassDep = 2;
275  const string names[nmbMassDep] = {
276  // "FLAT",
277  "BW",
278  // "AMP_M",
279  // "AMP_VES",
280  // "AMP_KACH",
281  "RHO_PRIME"
282  };
283  massDependencePtr massDep[nmbMassDep] = {
284  // createFlatMassDependence(),
285  createRelativisticBreitWigner(),
286  // createPiPiSWaveAuMorganPenningtonM(),
287  // createPiPiSWaveAuMorganPenningtonVes(),
288  // createPiPiSWaveAuMorganPenningtonKachaev(),
289  createRhoPrimeMassDep()
290  };
291  ::massDep* pwa2kMassDep[nmbMassDep] = {
292  // new flat (),
293  new breitWigner(),
294  // new AMP_M (),
295  // new AMP_ves (),
296  // new AMP_kach (),
297  new rhoPrime ()
298  };
299  for (size_t i = 0; i < nmbMassDep; ++i) {
300  compareAmplitudes(nmbEvents, seed, massDep[i], vertex, pwa2kMassDep[i], pwa2kParent, names[i]);
301  cout << endl;
302  }
303  outFile->Write();
304  outFile->Close();
305 }
Generated by   doxygen 1.8.1.2