File indexing completed on 2024-09-28 07:02:50
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "eicsmear/erhic/Kinematics.h"
0011
0012 #include <algorithm>
0013 #include <cmath>
0014 #include <functional>
0015 #include <iostream>
0016 #include <list>
0017 #include <memory>
0018 #include <numeric> // For std::accumulate
0019 #include <stdexcept>
0020 #include <utility>
0021 #include <vector>
0022
0023 #include <TDatabasePDG.h>
0024 #include <TParticlePDG.h>
0025 #include <TVector3.h>
0026
0027 #include "eicsmear/erhic/EventDis.h"
0028 #include "eicsmear/erhic/ParticleIdentifier.h"
0029
0030 bool erhic::DisKinematics::BoundaryWarning=true;
0031 namespace {
0032
0033 const double chargedPionMass =
0034 TDatabasePDG::Instance()->GetParticle(211)->Mass();
0035
0036 typedef erhic::VirtualParticle Particle;
0037
0038
0039
0040
0041 double computeW2FromXQ2M(double x, double Q2, double m) {
0042 if (x > 0.) {
0043 return std::pow(m, 2.) + (1. - x) * Q2 / x;
0044 }
0045 return 0.;
0046 }
0047
0048
0049
0050
0051 double bounded(double x, double minimum, double maximum) {
0052
0053
0054
0055
0056
0057
0058
0059
0060 if ( erhic::DisKinematics::BoundaryWarning && ( x< minimum || x > maximum ) ){
0061 std::cerr << "Warning in Kinematics, bounded(): x (or y) = " << x
0062 << " is outside [" << minimum << "," << maximum << "]" << std::endl;
0063 std::cerr << "To disable this warning, set erhic::DisKinematics::BoundaryWarning=false;" << std::endl;
0064 }
0065 return x;
0066
0067 }
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078 using erhic::ParticleMC;
0079 using erhic::VirtualParticle;
0080 typedef std::vector<const VirtualParticle*>::iterator VirtPartIter;
0081 typedef std::vector<const VirtualParticle*>::const_iterator cVirtPartIter;
0082
0083 class MeasuredParticle {
0084 public:
0085 static ParticleMC* Create(const erhic::VirtualParticle* particle) {
0086 if (!particle) {
0087 throw std::invalid_argument("MeasuredParticle given NULL pointer");
0088 }
0089 ParticleMC* measured = new ParticleMC;
0090
0091 measured->SetId(CalculateId(particle));
0092
0093 TParticlePDG* pdg = measured->Id().Info();
0094 if (pdg) {
0095 measured->SetM(pdg->Mass());
0096 }
0097 std::pair<double, double> ep =
0098 CalculateEnergyMomentum(particle, pdg->Mass());
0099 TLorentzVector vec(0., 0., ep.second, ep.first);
0100 vec.SetTheta(particle->GetTheta());
0101 vec.SetPhi(particle->GetPhi());
0102 measured->Set4Vector(vec);
0103 return measured;
0104 }
0105
0106
0107
0108
0109
0110
0111 static int CalculateId(const erhic::VirtualParticle* particle) {
0112 int id(0);
0113 TParticlePDG* pdg = particle->Id().Info();
0114
0115 if (pdg && particle->Id() != 0) {
0116 id = particle->Id();
0117 } else if (particle->GetP() > 0.) {
0118
0119
0120
0121
0122
0123
0124
0125
0126 id = 211;
0127 } else {
0128 id = 22;
0129 }
0130 return id;
0131 }
0132
0133
0134
0135
0136
0137
0138
0139
0140 static std::pair<double, double> CalculateEnergyMomentum(
0141 const erhic::VirtualParticle* particle, double mass = -1.) {
0142 if (mass < 0.) {
0143 int id = CalculateId(particle);
0144 TParticlePDG* pdg = TDatabasePDG::Instance()->GetParticle(id);
0145 if (pdg) {
0146 mass = pdg->Mass();
0147 } else {
0148 mass = 0.;
0149 }
0150 }
0151
0152
0153
0154 std::pair<double, double> ep(0., 0.);
0155 if (particle->GetP() > 0.) {
0156 ep.first = sqrt(pow(particle->GetP(), 2.) + pow(mass, 2.));
0157 ep.second = particle->GetP();
0158 } else if (particle->GetE() > 0.) {
0159 ep.first = particle->GetE();
0160
0161
0162 auto E = particle->GetE();
0163 if ( E >= mass ){
0164 ep.second = sqrt(pow(E, 2.) - pow(mass, 2.));
0165 } else {
0166 ep.second = 0;
0167 }
0168 }
0169 return ep;
0170 }
0171 };
0172
0173 }
0174
0175 namespace erhic {
0176
0177 DisKinematics::DisKinematics()
0178 : mX(0.)
0179 , mQ2(0.)
0180 , mW2(0.)
0181 , mNu(0.)
0182 , mY(0) {
0183 }
0184
0185 DisKinematics::DisKinematics(double x, double y, double nu,
0186 double Q2, double W2)
0187 : mX(x)
0188 , mQ2(Q2)
0189 , mW2(W2)
0190 , mNu(nu)
0191 , mY(y) {
0192 }
0193
0194
0195
0196 LeptonKinematicsComputer::LeptonKinematicsComputer(const EventDis& event) {
0197
0198 mBeams.push_back(event.BeamLepton());
0199 mBeams.push_back(event.BeamHadron());
0200 mBeams.push_back(event.ExchangeBoson());
0201 mBeams.push_back(event.ScatteredLepton());
0202 }
0203
0204
0205
0206 DisKinematics* LeptonKinematicsComputer::Calculate() {
0207
0208 DisKinematics* kin = new DisKinematics;
0209 try {
0210
0211
0212
0213 std::unique_ptr<const VirtualParticle> scattered( MeasuredParticle::Create(mBeams.at(3)));
0214
0215
0216
0217 if (scattered->GetTheta() > 0. && scattered->GetP() > 0.) {
0218 const TLorentzVector& l = mBeams.at(0)->Get4Vector();
0219 const TLorentzVector& h = mBeams.at(1)->Get4Vector();
0220
0221
0222
0223
0224 double Q2 = 2. * ( l.E() * scattered->GetE() + scattered->GetP()*l.P()*cos(scattered->GetTheta()) - l.M2()*l.M2() );
0225 kin->mQ2 = std::max(0., Q2);
0226
0227
0228 double gamma = mBeams.at(1)->GetE() / mBeams.at(1)->GetM();
0229 double beta = mBeams.at(1)->GetP() / mBeams.at(1)->GetE();
0230 double ELeptonInNucl = gamma * (l.P() - beta * l.Pz());
0231 double ELeptonOutNucl = gamma *
0232 (scattered->GetP() - beta * scattered->GetPz());
0233 kin->mNu = std::max(0., ELeptonInNucl - ELeptonOutNucl);
0234
0235
0236
0237
0238 const TLorentzVector q = l - scattered->Get4Vector();
0239 const double ldoth = l.Dot(h);
0240 double y = q.Dot(h) / ldoth ;
0241 if ( y<0) y=0;
0242 kin->mY = bounded(y, 0., 1.);
0243 double x = 0;
0244
0245
0246
0247
0248 if ( y>0 && std::abs(ldoth) > 0) x = kin->mQ2 / 2 / y / ldoth;
0249 if ( x<0 ) x=0;
0250 kin->mX = bounded(x, 0., 1.);
0251 kin->mW2 = computeW2FromXQ2M(kin->mX, kin->mQ2, h.M());
0252
0253
0254 }
0255 }
0256 catch(std::invalid_argument& except) {
0257
0258 std::cerr << "No lepton for kinematic calculations" << std::endl;
0259 }
0260 return kin;
0261 }
0262
0263
0264
0265 JacquetBlondelComputer::~JacquetBlondelComputer() {
0266
0267 for (VirtPartIter i = mParticles.begin(); i != mParticles.end(); ++i) {
0268 if (*i) {
0269 delete *i;
0270 *i = NULL;
0271 }
0272 }
0273 mParticles.clear();
0274 }
0275
0276
0277
0278 JacquetBlondelComputer::JacquetBlondelComputer(const EventDis& event)
0279 : mEvent(event) {
0280
0281
0282 std::vector<const erhic::VirtualParticle*> final;
0283 mEvent.HadronicFinalState(final);
0284
0285
0286
0287
0288
0289
0290
0291 for (VirtPartIter it = final.begin() ; it != final.end(); ++it){
0292 mParticles.push_back ( MeasuredParticle::Create ( *it ) );
0293 }
0294 }
0295
0296
0297
0298 DisKinematics* JacquetBlondelComputer::Calculate() {
0299
0300 DisKinematics* kin = new DisKinematics;
0301 kin->mY = ComputeY();
0302 kin->mQ2 = ComputeQSquared();
0303 kin->mX = ComputeX();
0304 kin->mW2 = computeW2FromXQ2M(kin->mX, kin->mQ2,
0305 mEvent.BeamHadron()->GetM());
0306 return kin;
0307 }
0308
0309
0310
0311 Double_t JacquetBlondelComputer::ComputeY() const {
0312 double y(0.);
0313
0314 const VirtualParticle* hadron = mEvent.BeamHadron();
0315 const VirtualParticle* lepton = mEvent.BeamLepton();
0316 if (hadron && lepton) {
0317
0318 std::vector<double> E;
0319 for (cVirtPartIter it = mParticles.begin() ; it != mParticles.end(); ++it){
0320 E.push_back ( (*it)->GetE() );
0321 }
0322 const double sumEh = std::accumulate(E.begin(), E.end(), 0.);
0323
0324
0325 std::vector<double> pz;
0326 for (cVirtPartIter it = mParticles.begin() ; it != mParticles.end(); ++it){
0327 pz.push_back ( (*it)->GetPz() );
0328 }
0329 const double sumPzh = std::accumulate(pz.begin(), pz.end(), 0.);
0330
0331
0332
0333
0334
0335 y = (hadron->GetE() * sumEh -
0336 hadron->GetPz() * sumPzh -
0337 pow(hadron->GetM(), 2.)) /
0338 (lepton->GetE() * hadron->GetE() - lepton->GetPz() * hadron->GetPz());
0339
0340
0341
0342 if (y<0) y=0;
0343 }
0344
0345
0346
0347
0348
0349
0350 return y;
0351 }
0352
0353
0354
0355
0356
0357 Double_t JacquetBlondelComputer::ComputeQSquared() const {
0358 double Q2(0.);
0359
0360 const VirtualParticle* hadron = mEvent.BeamHadron();
0361 if (hadron) {
0362
0363 std::vector<double> px;
0364 for (cVirtPartIter it = mParticles.begin() ; it != mParticles.end(); ++it){
0365 px.push_back ( (*it)->GetPx() );
0366 }
0367
0368
0369 std::vector<double> py;
0370 for (cVirtPartIter it = mParticles.begin() ; it != mParticles.end(); ++it){
0371 py.push_back ( (*it)->GetPy() );
0372 }
0373
0374 double sumPx = std::accumulate(px.begin(), px.end(), 0.);
0375 double sumPy = std::accumulate(py.begin(), py.end(), 0.);
0376 double y = ComputeY();
0377 if (y < 1.) {
0378 Q2 = (pow(sumPx, 2.) + pow(sumPy, 2.)) / (1. - y);
0379 }
0380 }
0381 return std::max(0., Q2);
0382 }
0383
0384
0385
0386
0387 Double_t JacquetBlondelComputer::ComputeX() const {
0388 double x(0.);
0389
0390 const VirtualParticle* hadron = mEvent.BeamHadron();
0391 const VirtualParticle* lepton = mEvent.BeamLepton();
0392 if (hadron && lepton) {
0393 double y = ComputeY();
0394 double s = (hadron->Get4Vector() + lepton->Get4Vector()).M2();
0395
0396 if (y > 0.0) {
0397 x = ComputeQSquared() / y / s;
0398 }
0399 }
0400
0401
0402
0403
0404
0405 return bounded(x, 0., 100.);
0406 }
0407
0408
0409
0410 DoubleAngleComputer::~DoubleAngleComputer() {
0411
0412 typedef std::vector<const VirtualParticle*>::iterator Iter;
0413 for (Iter i = mParticles.begin(); i != mParticles.end(); ++i) {
0414 if (*i) {
0415 delete *i;
0416 *i = NULL;
0417 }
0418 }
0419 mParticles.clear();
0420 }
0421
0422
0423
0424 DoubleAngleComputer::DoubleAngleComputer(const EventDis& event)
0425 : mEvent(event) {
0426
0427
0428 std::vector<const erhic::VirtualParticle*> final;
0429 mEvent.HadronicFinalState(final);
0430
0431
0432 for (VirtPartIter it = final.begin() ; it != final.end(); ++it){
0433 mParticles.push_back ( MeasuredParticle::Create ( *it ) );
0434 }
0435 }
0436
0437
0438
0439
0440 DisKinematics* DoubleAngleComputer::Calculate() {
0441 mHasChanged = true;
0442 DisKinematics* kin = new DisKinematics;
0443 kin->mQ2 = ComputeQSquared();
0444 kin->mX = ComputeX();
0445 kin->mY = ComputeY();
0446
0447 kin->mW2 = computeW2FromXQ2M(kin->mX, kin->mQ2, mEvent.BeamHadron()->GetM());
0448 return kin;
0449 }
0450
0451
0452
0453
0454
0455
0456
0457
0458
0459 Double_t DoubleAngleComputer::ComputeQuarkAngle() const {
0460
0461
0462 if (!mHasChanged) {
0463 return mAngle;
0464 }
0465 std::vector<TLorentzVector> hadrons;
0466
0467
0468
0469
0470 for (cVirtPartIter it = mParticles.begin() ; it != mParticles.end(); ++it){
0471 hadrons.push_back ( (*it)->Get4Vector() );
0472 }
0473 TLorentzVector h = std::accumulate(hadrons.begin(),
0474 hadrons.end(),
0475 TLorentzVector(0., 0., 0., 0.));
0476 mAngle = 2. * TMath::ATan((h.E() - h.Pz()) / h.Pt());
0477 mHasChanged = false;
0478 return mAngle;
0479 }
0480
0481
0482
0483 Double_t DoubleAngleComputer::ComputeY() const {
0484 double y(0.);
0485 const VirtualParticle* scattered = mEvent.ScatteredLepton();
0486 if (scattered) {
0487 double theta = scattered->GetTheta();
0488 double gamma = ComputeQuarkAngle();
0489 double denominator = tan(theta / 2.) + tan(gamma / 2.);
0490 if (denominator > 0.) {
0491 y = tan(gamma / 2.) / denominator;
0492 }
0493 }
0494
0495 return bounded(y, 0., 1.);
0496 }
0497
0498
0499
0500 Double_t DoubleAngleComputer::ComputeQSquared() const {
0501 double Q2(0.);
0502 const VirtualParticle* lepton = mEvent.BeamLepton();
0503 const VirtualParticle* scattered = mEvent.ScatteredLepton();
0504 if (lepton && scattered) {
0505 double theta = scattered->GetTheta();
0506 double gamma = ComputeQuarkAngle();
0507 double denominator = tan(theta / 2.) + tan(gamma / 2.);
0508
0509 if (denominator > 0. ) {
0510 Q2 = 4. * pow(lepton->GetE(), 2.) / tan(theta / 2.) / denominator;
0511 }
0512 }
0513
0514 return std::max(0., Q2);
0515 }
0516
0517
0518
0519 Double_t DoubleAngleComputer::ComputeX() const {
0520 double x(0.);
0521 const VirtualParticle* lepton = mEvent.BeamLepton();
0522 const VirtualParticle* hadron = mEvent.BeamHadron();
0523 if (lepton && hadron) {
0524 double s = (lepton->Get4Vector() + hadron->Get4Vector()).M2();
0525 double y = ComputeY();
0526 if (s > 0. && y > 0.) {
0527 x = ComputeQSquared() / y / s;
0528 }
0529 }
0530
0531 return bounded(x, 0., 1.);
0532 }
0533
0534 }