estarlight/src/eventchannel.cpp

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 255                         $: revision of last commit
0024 // $Author:: jnystrand                $: author of last commit
0025 // $Date:: 2016-04-06 14:01:46 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //    Class needed for root output
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032
0033
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cmath>
0037
0038 #include "eventchannel.h"
0039
0040
0041 using namespace std;
0042
0043
0044 //______________________________________________________________________________
0045 eventChannel::eventChannel(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem)
0046     : readLuminosity(inputParametersInstance),
0047       _bbs(bbsystem),
0048       _nmbAttempts(0),
0049       _nmbAccepted(0),
0050       _totalChannelCrossSection(0)
0051 {
0052   _ptCutEnabled  = inputParametersInstance.ptCutEnabled();
0053   _ptCutMin      = inputParametersInstance.ptCutMin();
0054   _ptCutMax      = inputParametersInstance.ptCutMax();
0055   _etaCutEnabled = inputParametersInstance.etaCutEnabled();
0056   _etaCutMin     = inputParametersInstance.etaCutMin();
0057   _etaCutMax     = inputParametersInstance.etaCutMax();
0058   _randy.SetSeed(inputParametersInstance.randomSeed());
0059 }
0060
0061
0062 //______________________________________________________________________________
0063 eventChannel::~eventChannel()
0064 { }
0065
0066
0067 //______________________________________________________________________________
0068 void
0069 eventChannel::transform(const double  betax,
0070                         const double  betay,
0071                         const double  betaz,
0072                         double&       E,
0073                         double&       px,
0074                         double&       py,
0075                         double&       pz,
0076                         int&          iFbadevent)
0077 {
0078   // carries out a lorentz transform of the frame.  (Not a boost!)???
0079   const double E0  = E;
0080   const double px0 = px;
0081   const double py0 = py;
0082   const double pz0 = pz;
0083
0084   const double beta = sqrt(betax * betax + betay * betay + betaz * betaz);
0085   if (beta >= 1)
0086       iFbadevent = 1;
0087   const double gamma = 1. / sqrt(1. - beta * beta);
0088   const double gob   = (gamma - 1) / (beta * beta);
0089
0090   E   = gamma * (E0 - betax * px0 - betay * py0 - betaz*  pz0);
0091   px  = -gamma * betax * E0 + (1. + gob * betax * betax) * px0
0092       + gob * betax * betay * py0 + gob * betax * betaz * pz0;
0093   py  = -gamma * betay * E0 + gob * betay * betax * px0
0094       + (1. + gob * betay * betay) * py0 + gob * betay * betaz  *pz0;
0095   pz  = -gamma * betaz * E0 + gob * betaz * betax * px0
0096       + gob * betaz * betay * py0 + (1. + gob * betaz * betaz) * pz0;
0097 }
0098
0099
0100 //______________________________________________________________________________
0101 double
0102 eventChannel::pseudoRapidity(const double px,
0103                              const double py,
0104                              const double pz)
0105 {
0106   const double pT= sqrt(px * px + py * py);
0107   const double p = sqrt(pz * pz + pT * pT);
0108   double eta = -99.9;  // instead of special value, std::numeric_limits<double>::quiet_NaN() should be used
0109   if ((p - pz) != 0)
0110       eta = 0.5 * log((p + pz)/(p - pz));
0111   return eta;
0112 }
0113