ROOTPWA
leptoProductionVertex.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 // class that describes leptoproduction vertex
29 // the kinematics is defined by the incoming beam lepton, the
30 // scattered lepton and the target; if the recoil particle is not
31 // specified, elastic scattering is assumed
32 //
33 // kinematic quantities are calculated based on COMPASS note
34 // 2009-6 using only Lorentz-scalars and without any
35 // approximations concerning lepton mass or large Q^2
36 //
37 // calculation of photon spin-density matrix is based on
38 // Schilling and Wolf, NPB 61, 381 (1973)
39 //
40 //
41 // Author List:
42 // Boris Grube TUM (original author)
43 //
44 //
45 //-------------------------------------------------------------------------
46 
47 
48 #include <cmath>
49 
50 #include "TClonesArray.h"
51 #include "TClass.h"
52 #include "TObjString.h"
53 #include "TVector3.h"
54 
55 #include "reportingUtils.hpp"
56 #include "reportingUtilsRoot.hpp"
57 #include "leptoProductionVertex.h"
58 
59 
60 using namespace std;
61 using namespace rpwa;
62 
63 
64 bool leptoProductionVertex::_debug = false;
65 
66 
67 leptoProductionVertex::leptoProductionVertex(const particlePtr& beamLepton,
68  const particlePtr& target,
69  const particlePtr& XParticle,
70  const particlePtr& recoil)
71  : productionVertex (),
72  _longPol (0),
73  _beamLeptonMomCache (),
74  _scatteredLeptonMomCache(),
75  _recoilMomCache (),
76  _targetMomCache ()
77 {
78  if (not beamLepton) {
79  printErr << "null pointer to beam lepton particle. aborting." << endl;
80  throw;
81  }
82  if (not target) {
83  printErr << "null pointer to target particle. aborting." << endl;
84  throw;
85  }
86  if (not XParticle) {
87  printErr << "null pointer to particle representing X system. aborting." << endl;
88  throw;
89  }
90  interactionVertex::addInParticle (beamLepton);
91  interactionVertex::addInParticle (target);
92  interactionVertex::addInParticle (createParticle("gamma")); // virtual photon
93  interactionVertex::addOutParticle(XParticle);
94  if (not recoil) {
95  if (_debug)
96  printWarn << "recoil not specified. assuming elastic scattering." << endl;
97  interactionVertex::addOutParticle(createParticle(*target));
98  }
99  interactionVertex::addOutParticle(createParticle(*beamLepton)); // lepton scatters elastically
100  if (_debug)
101  printDebug << "constructed " << *this << endl;
102 }
103 
104 
105 leptoProductionVertex::leptoProductionVertex(const leptoProductionVertex& vert)
106 {
107  *this = vert;
108 }
109 
110 
111 leptoProductionVertex::~leptoProductionVertex()
112 { }
113 
114 
115 leptoProductionVertex&
116 leptoProductionVertex::operator =(const leptoProductionVertex& vert)
117 {
118  if (this != &vert) {
119  interactionVertex::operator =(vert);
120  _beamLeptonMomCache = vert._beamLeptonMomCache;
121  _scatteredLeptonMomCache = vert._scatteredLeptonMomCache;
122  _recoilMomCache = vert._recoilMomCache;
123  _targetMomCache = vert._targetMomCache;
124  }
125  return *this;
126 }
127 
128 
129 leptoProductionVertex*
130 leptoProductionVertex::doClone(const bool cloneInParticles,
131  const bool cloneOutParticles) const
132 {
133  leptoProductionVertex* vertexClone = new leptoProductionVertex(*this);
134  if (cloneInParticles)
135  vertexClone->cloneInParticles();
136  if (cloneOutParticles)
137  vertexClone->cloneOutParticles();
138  if (_debug)
139  printDebug << "cloned " << *this << "; " << this << " -> " << vertexClone << " "
140  << ((cloneInParticles ) ? "in" : "ex") << "cluding incoming particles, "
141  << ((cloneOutParticles) ? "in" : "ex") << "cluding outgoing particles" << std::endl;
142  return vertexClone;
143 }
144 
145 
146 bool
147 leptoProductionVertex::addInParticle(const particlePtr&)
148 {
149  if (_debug)
150  printWarn << "cannot add incoming particle to " << *this << endl;
151  return false;
152 }
153 
154 
155 bool
156 leptoProductionVertex::addOutParticle(const particlePtr&)
157 {
158  if (_debug)
159  printWarn << "cannot add outgoing particle to " << *this << endl;
160  return false;
161 }
162 
163 
164 complex<double>
165 leptoProductionVertex::productionAmp() const
166 {
167  // calculate azimuthal angle between lepton-scattering and
168  // production plane in (virtual photon, target) CM system since
169  TVector3 k1, k2, q, v; // vectors of beam lepton, scattered lepton, virtual photon,
170  // and X particle used to define lepton-scattering and
171  // production plane
172  const TLorentzVector targetLv = target()->lzVec();
173  if (targetLv.Vect() == TVector3(0, 0, 0)) {
174  // fixed target case:
175  // since going from fixed target frame into (virtual photon,
176  // target) CM system involves only a boost along the virtual
177  // photon direction and since the normals of both the
178  // lepton-scattering and the production plane are perpendicular to
179  // the virtual photon, the azimuthal angle can be calculated in
180  // the lab frame as well
181  k1 = beamLepton()->momentum();
182  k2 = scatteredLepton()->momentum();
183  q = virtPhoton()->momentum();
184  v = XParticle()->momentum();
185  } else {
186  // general case
187  // boost vectors to (virtual photon, target) CM system
188  TLorentzVector beamLeptonLv = beamLepton()->lzVec();
189  TLorentzVector scatteredLeptonLv = scatteredLepton()->lzVec();
190  TLorentzVector virtPhotonLv = virtPhoton()->lzVec();
191  TLorentzVector XParticleLv = XParticle()->lzVec();
192  const TLorentzVector photonTargetCM = virtPhotonLv + targetLv;
193  const TVector3 cmBoost = photonTargetCM.BoostVector();
194  beamLeptonLv.Boost(cmBoost);
195  scatteredLeptonLv.Boost(cmBoost);
196  virtPhotonLv.Boost(cmBoost);
197  XParticleLv.Boost(cmBoost);
198  k1 = beamLeptonLv.Vect();
199  k2 = scatteredLeptonLv.Vect(); // use scattered lepton, because virtual photon is not directly measured
200  q = virtPhotonLv.Vect();
201  v = XParticleLv.Vect();
202  }
203  // calculate azimuthal angle in [-pi, +pi]
204  const TVector3 leptonPlaneNormal = k1.Cross(k2);
205  const TVector3 productionPlaneNormal = q.Cross(v);
206  const double parProjection = leptonPlaneNormal.Dot(productionPlaneNormal);
207  const double perpProjection = leptonPlaneNormal.Cross(productionPlaneNormal).Mag();
208  const double phi = atan2(perpProjection, parProjection);
209 
210  // compute some kinematic variables
211  const double epsilon = this->epsilon();
212  const double delta = this->delta(epsilon);
213  printInfo << "phi = " << phi << ", epsilon = " << maxPrecision(epsilon) << ", "
214  << "delta = " << maxPrecision(delta) << ", pol = " << _longPol << endl;
215  // compute some intermediary terms
216  const double xi = sqrt(epsilon * (1 + epsilon + 2 * delta));
217  const double zeta = sqrt(xi * (1 - epsilon) / (1 + epsilon));
218  const double term = _longPol * sqrt(1 - epsilon * epsilon);
219  const complex<double> phase = exp(complex<double>(0, phi));
220  printInfo << "xi = " << xi << ", zeta = " << zeta << ", term = " << term << endl;
221 
222  // define lower triangle of virtual photon spin density matrix
223  // rho_{lambda, lambda'}; lambda = -1, 0, 1 as in Schilling and
224  // Wolf (common factor 1/2 is dropped):
225  // eq. 44 and 57 with 60, where alpha_2 was set to 0 (long. lepton
226  // polarization), into eq. 62
227  complex<double> rho[3][3];
228  // column with lambda' = -1
229  rho[0][0] = 1 + term; // lambda = -1
230  rho[1][0] = (_longPol * zeta + xi) * phase; // lambda = 0
231  rho[2][0] = -epsilon * exp(complex<double>(0, 2 * phi)); // lambda = +1
232  // column with lambda' = 0
233  rho[1][1] = 2 * (epsilon + delta); // lambda = 0
234  rho[2][1] = (_longPol * zeta - xi) * phase; // lambda = +1
235  // column with lambda' = +1
236  rho[2][2] = 1 - term; // lambda = +1
237  // conjugated elements
238  rho[0][1] = conj(rho[1][0]);
239  rho[0][2] = conj(rho[2][0]);
240  rho[1][2] = conj(rho[2][1]);
241  for (unsigned int j = 0; j < 3; ++j)
242  for (unsigned int i = 0; i < 3; ++i)
243  printInfo << "rho[" << i << "][" << j << "] = " << maxPrecisionDouble(rho[i][j]) << endl;
244  const complex<double> detRho = rho[0][0] * rho[1][1] * rho[2][2]
245  + rho[0][1] * rho[1][2] * rho[2][0]
246  + rho[1][0] * rho[2][1] * rho[0][2]
247  - rho[0][2] * rho[1][1] * rho[2][0]
248  - rho[1][2] * rho[0][0] * rho[2][1]
249  - rho[0][1] * rho[2][2] * rho[1][0];
250  printInfo << "det[rho] = " << detRho << " vs. "
251  << ( (epsilon + delta) * (1 - epsilon * epsilon) * (1 - _longPol * _longPol)
252  - (1 + epsilon) * (xi * xi + _longPol * _longPol * zeta * zeta)
253  + 2 * _longPol * _longPol * zeta * xi * sqrt(1 - epsilon * epsilon)) / 4
254  << endl;
255 
256  // perform Cholesky decomposition rho_ij = sum_r V_ir * V_jr^*, where V_ir is a lower
257  // triangle matrix with real diagonal elements
258  complex<double> V[3][3];
259  // first column
260  V[0][0] = sqrt(real(rho[0][0]));
261  V[1][0] = rho[1][0] / real(V[0][0]);
262  V[2][0] = rho[2][0] / real(V[0][0]);
263  // second column
264  V[1][1] = sqrt(real(rho[1][1]) - norm(V[1][0]));
265  V[2][1] = (rho[2][1] - V[2][0] * conj(V[1][0])) / real(V[1][1]);
266  // third column
267  V[2][2] = sqrt(real(rho[2][2]) - norm(V[2][1]) - norm(V[2][0]));
268  // zero elements
269  V[0][1] = 0;
270  V[0][2] = 0;
271  V[1][2] = 0;
272  printInfo << "V[2][2]^2 = " << real(rho[2][2]) - norm(V[2][1]) - norm(V[2][0]) << ": "
273  << real(rho[2][2]) << " - " << norm(V[2][1]) << " - " << norm(V[2][0]) << endl;
274  for (unsigned int j = 0; j < 3; ++j)
275  for (unsigned int i = 0; i < 3; ++i)
276  printInfo << "V[" << i << "][" << j << "] = " << maxPrecisionDouble(V[i][j]) << endl;
277  complex<double> rhoPrime[3][3];
278  for (unsigned int j = 0; j < 3; ++j)
279  for (unsigned int i = 0; i < 3; ++i) {
280  rhoPrime[i][j] = 0;
281  for (unsigned int r = 0; r < 3; ++r)
282  rhoPrime[i][j] += V[i][r] * conj(V[j][r]);
283  }
284  for (unsigned int j = 0; j < 3; ++j)
285  for (unsigned int i = 0; i < 3; ++i)
286  printInfo << "deltaRho[" << i << "][" << j << "] = "
287  << maxPrecisionDouble(rho[i][j] - rhoPrime[i][j]) << endl;
288 
289 
290  // compute production amplitude for given photon helicity
291  complex<double> prodAmp = 0;
292  return prodAmp;
293 }
294 
295 
296 double
297 leptoProductionVertex::epsilon() const
298 {
299  const double Q2 = this->Q2();
300  const double xBj = this->xBj();
301  const double xBj2 = xBj * xBj;
302  const double y = this->y();
303  const double y2 = y * y;
304 
305  // calculate kinematic values based on Lorentz invariants
306  // see COMPASS note 2009-6
307  const double gamma2 = 4 * xBj2 * target()->mass2() / Q2; // eq. 77d^2
308  const double term1 = beamLepton()->mass2() * y2 * gamma2 / Q2;
309  const double term2 = 1 / ((1 - y) * (1 - y));
310  const double u0 = sqrt(1 - term1 * (1 + term2) + term1 * term1 * term2); // eq. 81
311  const double E1E2 = Q2 * (1 - y) / (y2 * gamma2); // eq. 80
312  const double Q2Min = 2 * (E1E2 * (1 - u0) - beamLepton()->mass2()); // eq. 5c and 82
313  const double k1k2 = E1E2 * u0; // eq. 82
314  const double sin2ThetaHalf = (Q2 - Q2Min) / (4 * k1k2); // from eq. 5b
315  const double tan2ThetaHalf = sin2ThetaHalf / (1 - sin2ThetaHalf);
316  const double term3 = 1 - Q2Min / Q2;
317  const double oneOverEpsilon = 1 + 2 * (1 + 1 / gamma2) * tan2ThetaHalf / (term3 * term3); // eq. 31 and 77d
318  return 1 / oneOverEpsilon;
319 
320  // // calculate kinematic values based on Lorentz invariants
321  // // see COMPASS note 2009-6
322  // const double gamma2 = 4 * xBj2 * target()->mass2() / Q2; // eq. 77d^2
323  // const double term1 = beamLepton()->mass2() * y2 * gamma2 / Q2;
324  // const double term2 = 1 / ((1 - y) * (1 - y));
325  // const double u0 = sqrt(1 - term1 * (1 + term2) + term1 * term1 * term2); // eq. 81
326  // const double v0 = 1 + 2 * ( beamLepton()->mass2() / Q2
327  // - (1 - y) * (1 - u0) / (y2 * gamma2)); // eq. 85
328  // const double v02 = v0 * v0;
329  // const double term3 = (1 - y) * u0 * v0;
330  // const double term4 = y2 * gamma2 / 4;
331  // return (term3 - term4 * v02) / (term3 + y2 / 2 + term4 * (2 - v02)); // eq. 87
332 }
333 
334 
335 void
336 leptoProductionVertex::setXFlavorQN()
337 {
339  particle& X = *XParticle();
340  particle& beam = *virtPhoton();
341  X.setBaryonNmb (beam.baryonNmb());
342  X.setStrangeness(beam.strangeness());
343  X.setCharm (beam.charm());
344  X.setBeauty (beam.beauty());
345 }
346 
347 
348 bool
349 leptoProductionVertex::initKinematicsData(const TClonesArray& prodKinPartNames)
350 {
351  _nmbProdKinPart = 0;
352 
353  // check production vertex data
354  const string partClassName = prodKinPartNames.GetClass()->GetName();
355  if (partClassName != "TObjString") {
356  printWarn << "production kinematics particle names are of type '" << partClassName
357  << "' and not TObjString." << endl;
358  return false;
359  }
360  _nmbProdKinPart = prodKinPartNames.GetEntriesFast();
361  if (_nmbProdKinPart < 2) {
362  printWarn << "array of production kinematics particle names has wrong size: "
363  << _nmbProdKinPart << ". need at least beam lepton (index 0) and "
364  << "scattered lepton (index 1); "
365  << "recoil (index 2) target (index 3) are optional." << endl;
366  return false;
367  }
368 
369  // beam lepton at index 0
370  bool success = true;
371  const string beamLeptonName = ((TObjString*)prodKinPartNames[0])->GetString().Data();
372  if (beamLeptonName != beamLepton()->name()) {
373  printWarn << "wrong particle at index 0 in production kinematics input data: "
374  << "read '" << beamLeptonName << "', "
375  << "expected beam lepton '" << beamLepton()->name() << "'" << endl;
376  success = false;
377  }
378 
379  // scattered lepton at index 1
380  const string scatteredLeptonName = ((TObjString*)prodKinPartNames[1])->GetString().Data();
381  if (scatteredLeptonName != scatteredLepton()->name()) {
382  printWarn << "wrong particle at index 1 in production kinematics input data: "
383  << "read '" << scatteredLeptonName << "', "
384  << "expected scattered lepton '" << scatteredLepton()->name() << "'" << endl;
385  success = false;
386  }
387 
388  // recoil at index 2 (optional)
389  if (_nmbProdKinPart >= 3) {
390  const string recoilName = ((TObjString*)prodKinPartNames[2])->GetString().Data();
391  if (recoilName != recoil()->name()) {
392  printWarn << "wrong particle at index 2 in production kinematics input data: "
393  << "read '" << recoilName << "', "
394  << "expected recoil particle '" << recoil()->name() << "'" << endl;
395  success = false;
396  }
397  }
398 
399  // target at index 3 (optional)
400  if (_nmbProdKinPart >= 4) {
401  const string targetName = ((TObjString*)prodKinPartNames[3])->GetString().Data();
402  if (targetName != target()->name()) {
403  printWarn << "wrong particle at index 3 in production kinematics input data: "
404  << "read '" << targetName << "', "
405  << "expected target particle '" << target()->name() << "'" << endl;
406  success = false;
407  }
408  }
409 
410  return success;
411 }
412 
413 
414 bool
415 leptoProductionVertex::readKinematicsData(const TClonesArray& prodKinMomenta)
416 {
417  _beamLeptonMomCache = TVector3();
418  _scatteredLeptonMomCache = TVector3();
419  _recoilMomCache = TVector3();
420  _targetMomCache = TVector3();
421 
422  // check production vertex data
423  const int nmbProdKinMom = prodKinMomenta.GetEntriesFast();
424  if (nmbProdKinMom != _nmbProdKinPart) {
425  printWarn << "array of production kinematics particle momenta has wrong size: "
426  << nmbProdKinMom << " (expected " << _nmbProdKinPart << "). "
427  << "cannot read production kinematics." << endl;
428  return false;
429  }
430 
431  // set beam lepton
432  bool success = true;
433  TVector3* beamLeptonMom = dynamic_cast<TVector3*>(prodKinMomenta[0]);
434  if (beamLeptonMom) {
435  if (_debug)
436  printDebug << "setting momentum of beam lepton '" << beamLepton()->name()
437  << "' to " << *beamLeptonMom << " GeV" << endl;
438  beamLepton()->setMomentum(*beamLeptonMom);
439  _beamLeptonMomCache = beamLepton()->momentum();
440  } else {
441  printWarn << "production kinematics data entry [0] is not of type TVector3. "
442  << "cannot read beam lepton momentum." << endl;
443  success = false;
444  }
445 
446  // set scattered lepton
447  TVector3* scatteredLeptonMom = dynamic_cast<TVector3*>(prodKinMomenta[1]);
448  if (scatteredLeptonMom) {
449  if (_debug)
450  printDebug << "setting momentum of scattered lepton '" << beamLepton()->name()
451  << "' to " << *beamLeptonMom << " GeV" << endl;
452  scatteredLepton()->setMomentum(*scatteredLeptonMom);
453  _scatteredLeptonMomCache = scatteredLepton()->momentum();
454  } else {
455  printWarn << "production kinematics data entry [1] is not of type TVector3. "
456  << "cannot read scattered lepton momentum." << endl;
457  success = false;
458  }
459 
460 
461  // set recoil (optional)
462  if (_nmbProdKinPart >= 3) {
463  TVector3* recoilMom = dynamic_cast<TVector3*>(prodKinMomenta[2]);
464  if (recoilMom) {
465  if (_debug)
466  printDebug << "setting momentum of recoil particle '" << recoil()->name()
467  << "' to " << *recoilMom << " GeV" << endl;
468  recoil()->setMomentum(*recoilMom);
469  _recoilMomCache = recoil()->momentum();
470  } else {
471  printWarn << "production kinematics data entry [2] is not of type TVector3. "
472  << "cannot read recoil particle momentum." << endl;
473  success = false;
474  }
475  }
476 
477  // set target (optional); if not defined fixed target is assumed
478  if (_nmbProdKinPart >= 4) {
479  TVector3* targetMom = dynamic_cast<TVector3*>(prodKinMomenta[3]);
480  if (targetMom) {
481  if (_debug)
482  printDebug << "setting momentum of target particle '" << target()->name()
483  << "' to " << *targetMom << " GeV" << endl;
484  target()->setMomentum(*targetMom);
485  _targetMomCache = target()->momentum();
486  } else {
487  printWarn << "production kinematics data entry [3] is not of type TVector3. "
488  << "cannot read target particle momentum." << endl;
489  success = false;
490  }
491  }
492 
493  // set virtual photon
494  virtPhoton()->setLzVec(beamLepton()->lzVec() - scatteredLepton()->lzVec());
495  return success;
496 }
497 
498 
499 bool
500 leptoProductionVertex::revertMomenta()
501 {
502  if (_debug) {
503  printDebug << "resetting beam lepton momentum to " << _beamLeptonMomCache << " GeV" << endl
504  << " resetting scattered lepton momentum to "
505  << _scatteredLeptonMomCache << " GeV" << endl
506  << " resetting recoil momentum to " << _recoilMomCache << " GeV" << endl
507  << " resetting target momentum to " << _targetMomCache << " GeV" << endl;
508  }
509  beamLepton ()->setMomentum(_beamLeptonMomCache );
510  scatteredLepton()->setMomentum(_scatteredLeptonMomCache);
511  recoil ()->setMomentum(_recoilMomCache );
512  target ()->setMomentum(_targetMomCache );
513  // set virtual photon
514  virtPhoton()->setLzVec(beamLepton()->lzVec() - scatteredLepton()->lzVec());
515  return true;
516 }
517 
518 
519 ostream&
520 leptoProductionVertex::print(ostream& out) const
521 {
522  out << name() << ": "
523  << "beam " << beamLepton()->qnSummary() << " P_L = " << _longPol << " -> "
524  << "scattered " << scatteredLepton()->qnSummary() << " + virtual " << virtPhoton()->qnSummary()
525  << " + target " << target()->qnSummary() << " ---> "
526  << XParticle()->qnSummary() << " + recoil " << recoil()->qnSummary();
527  return out;
528 }
529 
530 
531 ostream&
532 leptoProductionVertex::dump(ostream& out) const
533 {
534  out << name() << ": " << endl
535  << " beam lepton ........ " << *beamLepton() << endl
536  << " target ............. " << *target() << endl
537  << " virtual photon ..... " << *virtPhoton() << endl
538  << " scattered lepton ... " << *scatteredLepton() << endl
539  << " X .................. " << *XParticle() << endl
540  << " recoil ............. " << *recoil() << endl;
541  return out;
542 }
543 
544 
545 ostream&
546 leptoProductionVertex::printPointers(ostream& out) const
547 {
548  out << name() << " " << this << ": "
549  << "beam lepton = " << beamLepton() << "; "
550  << "target particle = " << target() << "; "
551  << "virtual photon = " << virtPhoton() << "; "
552  << "scattered lepton = " << scatteredLepton() << "; "
553  << "X particle = " << XParticle() << "; "
554  << "recoil particle = " << recoil() << endl;
555  return out;
556 }
Generated by   doxygen 1.8.1.2