Index: trunk/src/CrossSection.h =================================================================== --- trunk/src/CrossSection.h (revision 436) +++ trunk/src/CrossSection.h (revision 437) @@ -1,104 +1,106 @@ //============================================================================== // CrossSection.h // // Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich // // This file is part of Sartre. // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation. // This program is distributed in the hope that it will be useful, // but without any warranty; without even the implied warranty of // merchantability or fitness for a particular purpose. See the // GNU General Public License for more details. // You should have received a copy of the GNU General Public License // along with this program. If not, see . // // Author: Thomas Ullrich // Last update: // $Date$ // $Author$ //============================================================================== // // Functor class. // // operator() returns d^3sig/(dt dQ2 dW2) in nb/GeV^6. // //=============================================================================== #ifndef CrossSection_h #define CrossSection_h #include "Enumerations.h" #include "PhotonFlux.h" class TRandom3; class EventGeneratorSettings; class TableCollection; class CrossSection { public: CrossSection(TableCollection* = 0, TableCollection* = 0); ~CrossSection(); double operator()(double t, double Q2, double W2); double operator()(double t, double xpom); // UPC version double operator()(const double*); // array of t, Q2, W2 or t, xpom for UPC double unuranPDF(const double*); // for UNU.RAN using log(Q2) (or log(xpom) for UPC) // and returning log of cross-section void setTableCollection(TableCollection*); void setProtonTableCollection(TableCollection*); GammaPolarization polarizationOfLastCall() const; DiffractiveMode diffractiveModeOfLastCall() const; - + double crossSectionRatioLTOfLastCall() const; + void setCheckKinematics(bool); double dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization) const; protected: friend class Sartre; // mostly for debugging and QA double dsigdtdQ2dW2_total_checked(double t, double Q2, double W2); double dsigdtdxp_total_checked(double t, double xpom); double dsigdt_total(double t, double Q2, double W2, GammaPolarization) const; // modified double dsigdt_coherent(double t, double Q2, double W2, GammaPolarization) const; double dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization) const; // new double dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization) const; double dsigdt_total(double t, double xpom) const; // UPC version double dsigdt_coherent(double t, double xpom) const; // UPC version double dsigdt_incoherent(double t, double xpom) const; // UPC version double dsigdtdxp_total(double t, double xpom) const; // UPC double dsigdtdxp_coherent(double t, double xpom) const; // UPC double logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization) const; double logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization) const; double realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const; double skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const; double logDerivateOfAmplitude(double t, double xpom) const; // UPC version double logDerivateOfGluonDensity(double t, double xpom) const; // UPC version double realAmplitudeCorrection(double t, double xpom) const; // UPC version double skewednessCorrection(double t, double xpom) const; // UPC version double skewednessCorrection(double lambda) const; double realAmplitudeCorrection(double lambda) const; double UPCPhotonFlux(double t, double xpom) const; private: TRandom3* mRandom; GammaPolarization mPolarization; DiffractiveMode mDiffractiveMode; PhotonFlux mPhotonFlux; TableCollection* mTableCollection; TableCollection* mProtonTableCollection; EventGeneratorSettings* mSettings; double mS; - double mVmMass; - + double mVmMass; + double mCrossSectionRatioLT; + bool mCheckKinematics; }; #endif Index: trunk/src/CrossSection.cpp =================================================================== --- trunk/src/CrossSection.cpp (revision 436) +++ trunk/src/CrossSection.cpp (revision 437) @@ -1,921 +1,926 @@ //============================================================================== // CrossSection.cpp // // Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich // // This file is part of Sartre. // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation. // This program is distributed in the hope that it will be useful, // but without any warranty; without even the implied warranty of // merchantability or fitness for a particular purpose. See the // GNU General Public License for more details. // You should have received a copy of the GNU General Public License // along with this program. If not, see . // // Author: Thomas Ullrich // Last update: // $Date$ // $Author$ //============================================================================== #include "CrossSection.h" #include "EventGeneratorSettings.h" #include "Kinematics.h" #include "TableCollection.h" #include "Table.h" #include "Constants.h" #include "TH2F.h" #include "TFile.h" #include "DglapEvolution.h" #include #include #include #define PR(x) cout << #x << " = " << (x) << endl; CrossSection::CrossSection(TableCollection* tc, TableCollection* ptc) { mSettings = EventGeneratorSettings::instance(); mRandom = mSettings->randomGenerator(); mS = Kinematics::s(mSettings->eBeam(), mSettings->hBeam()); mPhotonFlux.setS(mS); mTableCollection = tc; mProtonTableCollection = ptc; TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mSettings->vectorMesonId()); mVmMass = vectorMesonPDG->Mass(); - mCheckKinematics = true; + mCheckKinematics = true; + mCrossSectionRatioLT = 0; } CrossSection::~CrossSection() { /* no op */ } void CrossSection::setTableCollection(TableCollection* tc) {mTableCollection = tc;} void CrossSection::setProtonTableCollection(TableCollection* ptc) {mProtonTableCollection = ptc;} void CrossSection::setCheckKinematics(bool val) {mCheckKinematics = val;} double CrossSection::operator()(double t, double Q2, double W2) { return dsigdtdQ2dW2_total_checked(t, Q2, W2); } // UPC Version double CrossSection::operator()(double t, double xpom) { return dsigdtdxp_total_checked(t, xpom); } double CrossSection::operator()(const double* array) { if (mSettings->UPC()) return dsigdtdxp_total_checked(array[0], array[1]); else return dsigdtdQ2dW2_total_checked(array[0], array[1], array[2]); } // // PDF passed to UNURAN // Array holds: t, log(Q2), W2 for e+p/A running // t, log(xpom) for UPC // double CrossSection::unuranPDF(const double* array) // array is t, log(Q2), W2 { // or t and log(xpom) double result = 0; if (mSettings->UPC()) { double xpom = exp(array[1]); result = dsigdtdxp_total_checked(array[0], xpom); // t, xpom result *= xpom; // Jacobian } else { double Q2 = exp(array[1]); result = dsigdtdQ2dW2_total_checked(array[0], Q2, array[2]); // t, Q2, W2 result *= Q2; // Jacobian } return log(result); } double CrossSection::dsigdtdQ2dW2_total_checked(double t, double Q2, double W2) { double result = 0; // // Check if in kinematically allowed region // if (mCheckKinematics && !Kinematics::valid(mS, t, Q2, W2, mVmMass, false, (mSettings->verboseLevel() > 2) )) { if (mSettings->verboseLevel() > 2) cout << "CrossSection::dsigdtdQ2dW2_total_checked(): warning, t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << " is outside of kinematically allowed region. Return 0." << endl; return result; } // // Total cross-section dsig2/(dQ2 dW2 dt) // This is the probability density function needed for UNU.RAN // double csT = dsigdtdQ2dW2_total(t, Q2, W2, transverse); double csL = dsigdtdQ2dW2_total(t, Q2, W2, longitudinal); - result = csT + csL; + result = csT + csL; + mCrossSectionRatioLT = csL/csT; // // Polarization // if (mRandom->Uniform(result) <= csT) mPolarization = transverse; else mPolarization = longitudinal; // // Diffractive Mode // double sampleRange = (mPolarization == transverse ? csT : csL); double sampleDivider = dsigdtdQ2dW2_coherent(t, Q2, W2, mPolarization); if (mRandom->Uniform(sampleRange) <= sampleDivider) mDiffractiveMode = coherent; else mDiffractiveMode = incoherent; // // Print-out at high verbose levels // if (mSettings->verboseLevel() > 10) { // Spinal Tap ;-) cout << "CrossSection::dsigdtdQ2dW2_total_checked(): " << result; cout << " at t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2); cout << " (" << (mPolarization == transverse ? "transverse" : "longitudinal"); cout << " ," << (mDiffractiveMode == coherent ? "coherent" : "incoherent"); cout << ')' << endl; } // // Check validity of return value // if (std::isnan(result)) { cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = NaN at" << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; result = 0; } if (std::isinf(result)) { cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, return value = inf at" << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; result = 0; } if (result < 0) { cout << "CrossSection::dsigdtdQ2dW2_total_checked(): Error, negative cross-section at" << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; result = 0; } return result; } double CrossSection::dsigdtdxp_total_checked(double t, double xpom) { double result = 0; // // Check if in kinematically allowed region // if (mCheckKinematics && !Kinematics::validUPC(mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy(), t, xpom, mVmMass, (mSettings->verboseLevel() > 2) )) { if (mSettings->verboseLevel() > 2) cout << "CrossSection::dsigdtdxp_total_checked(): warning, t=" << t << ", xpom=" << xpom << " is outside of kinematically allowed region. Return 0." << endl; return result; } // // Total cross-section dsig2/(dt dxp) // This is the probability density function needed for UNU.RAN // result = dsigdtdxp_total(t, xpom); // // Polarization // mPolarization = transverse; // always // // Diffractive Mode // double sampleRange = result; double sampleDivider = dsigdtdxp_coherent(t, xpom); if (mRandom->Uniform(sampleRange) <= sampleDivider) mDiffractiveMode = coherent; else mDiffractiveMode = incoherent; // // Print-out at high verbose levels // if (mSettings->verboseLevel() > 10) { cout << "CrossSection::dsigdtdxp_total_checked(): " << result; cout << " at t=" << t << ", xp=" << xpom; cout << " (" << (mDiffractiveMode == coherent ? "coherent" : "incoherent"); cout << ')' << endl; } // // Check validity of return value // if (std::isnan(result)) { cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = NaN at" << endl; cout << " t=" << t << ", xp=" << xpom << endl; result = 0; } if (std::isinf(result)) { cout << "CrossSection::dsigdtdxp_total_checked(): Error, return value = inf at" << endl; cout << " t=" << t << ", xp=" << xpom << endl; result = 0; } if (result < 0) { cout << "CrossSection::dsigdtdxp_total_checked(): Error, negative cross-section at" << endl; cout << " t=" << t << ", xp=" << xpom << endl; result = 0; } return result; } double CrossSection::dsigdt_total(double t, double Q2, double W2, GammaPolarization pol) const { double result = 0; if (mSettings->tableSetType() == coherent_and_incoherent) { result = dsigdt_coherent(t, Q2, W2, pol) + dsigdt_incoherent(t, Q2, W2, pol); } else if (mSettings->tableSetType() == total_and_coherent) { result = mTableCollection->get(Q2, W2, t, pol, mean_A2); // units now are fm^4 result /= (16*M_PI); result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()){ lambda = logDerivateOfAmplitude(t, Q2, W2, pol); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()){ lambda = logDerivateOfGluonDensity(t, Q2, W2, pol); result *= skewednessCorrection(lambda); } } return result; } // // UPC Version // double CrossSection::dsigdt_total(double t, double xpom) const { double result = 0; if (mSettings->tableSetType() == coherent_and_incoherent) { result = dsigdt_coherent(t, xpom) + dsigdt_incoherent(t, xpom); } else if (mSettings->tableSetType() == total_and_coherent) { result = mTableCollection->get(xpom, t, mean_A2); // units now are fm^4 result /= (16*M_PI); result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()) { lambda = logDerivateOfAmplitude(t, xpom); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()) { lambda = logDerivateOfGluonDensity(t, xpom); result *= skewednessCorrection(lambda); } } return result; } double CrossSection::dsigdt_coherent(double t, double Q2, double W2, GammaPolarization pol) const { double val = mTableCollection->get(Q2, W2, t, pol, mean_A); // fm^2 double result = val*val/(16*M_PI); // units now are fm^4 result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()) { lambda = logDerivateOfAmplitude(t, Q2, W2, pol); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()) { lambda = logDerivateOfGluonDensity(t, Q2, W2, pol); result *= skewednessCorrection(lambda); } return result; } // // UPC Version // double CrossSection::dsigdt_coherent(double t, double xpom) const { double val = mTableCollection->get(xpom, t, mean_A); // fm^2 double result = val*val/(16*M_PI); // units now are fm^4 result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()) { lambda = logDerivateOfAmplitude(t, xpom); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()) { lambda = logDerivateOfGluonDensity(t, xpom); result *= skewednessCorrection(lambda); } return result; } double CrossSection::dsigdtdQ2dW2_total(double t, double Q2, double W2, GammaPolarization pol) const { double result = dsigdt_total(t, Q2, W2, pol); if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol); return result; } // // UPC version // double CrossSection::dsigdtdxp_total(double t, double xpom) const { double result = dsigdt_total(t, xpom); if (mSettings->applyPhotonFlux()) { result *= UPCPhotonFlux(t, xpom); } return result; } double CrossSection::dsigdt_incoherent(double t, double Q2, double W2, GammaPolarization pol) const { double result = mTableCollection->get(Q2, W2, t, pol, variance_A); // fm^4 result /= (16*M_PI); result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()){ lambda = logDerivateOfAmplitude(t, Q2, W2, pol); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()){ lambda = logDerivateOfGluonDensity(t, Q2, W2, pol); result *= skewednessCorrection(lambda); } return result; } // // UPC Version // double CrossSection::dsigdt_incoherent(double t, double xpom) const { double result = mTableCollection->get(xpom, t, variance_A); // fm^4 result /= (16*M_PI); result /= hbarc2; // in fm^2/GeV^2 result *= 1e7; // in nb/GeV^2 double lambda = 0; if (mSettings->correctForRealAmplitude()) { lambda = logDerivateOfAmplitude(t, xpom); result *= realAmplitudeCorrection(lambda); } if (mSettings->correctSkewedness()) { lambda = logDerivateOfGluonDensity(t, xpom); result *= skewednessCorrection(lambda); } return result; } double CrossSection::dsigdtdQ2dW2_coherent(double t, double Q2, double W2, GammaPolarization pol) const { double result = dsigdt_coherent(t, Q2, W2, pol); if (mSettings->applyPhotonFlux()) result *= mPhotonFlux(Q2,W2,pol); return result; } // // UPC Version // double CrossSection::dsigdtdxp_coherent(double t, double xpom) const { double result = dsigdt_coherent(t, xpom); if (mSettings->applyPhotonFlux()) { result *= UPCPhotonFlux(t, xpom); } return result; } double CrossSection::UPCPhotonFlux(double t, double xpom) const{ double Egamma = Kinematics::Egamma(xpom, t, mVmMass, mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy()); double result = mPhotonFlux(Egamma); // // Jacobian = dEgamma/dN, such that dsig/dtdxp = dsigd/tdEgam*dEgam/dxp // double h = xpom*1e-3; double Egamma_p = Kinematics::Egamma(xpom+h, t, mVmMass, mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy()); double Egamma_m = Kinematics::Egamma(xpom-h, t, mVmMass, mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy()); double jacobian = (Egamma_p-Egamma_m)/(2*h); result *= fabs(jacobian); return result; } GammaPolarization CrossSection::polarizationOfLastCall() const {return mPolarization;} -DiffractiveMode CrossSection::diffractiveModeOfLastCall() const {return mDiffractiveMode;} +DiffractiveMode CrossSection::diffractiveModeOfLastCall() const {return mDiffractiveMode;} + +double CrossSection::crossSectionRatioLTOfLastCall() const {return mCrossSectionRatioLT;} + double CrossSection::logDerivateOfAmplitude(double t, double Q2, double W2, GammaPolarization pol) const { double lambda = 0; bool lambdaFromTable = true; Table *table = 0; if (!mProtonTableCollection) { cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl; cout << " Corrections not available. Should be off." << endl; return 0; } // // If the lambda table is present we use the more accurate and numerically // stable table value. Otherwise we calculate it from the table(s). // // Note (TU): if the lambda value from a table is not valid and not > 0, // get() returns lambda=0 and table=0. This enforces a renewed // calculation. // lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_real, table); if (!table) { // no lambda value from correct table, use fallback solution lambdaFromTable = false; double value = mProtonTableCollection->get(Q2, W2, t, pol, mean_A, table); // use obtained table from here on if (value <= 0) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << value << '.' << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << ", pol=" << (pol == transverse ? 'T' : 'L') << endl; return 0; } if (!table) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl; return 0; } // // Note: the derivate taken from values in the table is a delicate issue. // The standard interpolation method(s) used in Table::get() are // at times not accurate and causes ripples on lambda(W2). // However, interpolation methods or robust fitting is by far too // expensive in terms of CPU time per point. // // // Calculate the derivative using simple method. // double derivative; double hplus, hminus; double dW2 = table->binWidthW2(); hplus = hminus = 0.5*dW2; // half bin width is found to be the best choice after some testing hplus = min(hplus, table->maxW2()-W2); hminus = min(hminus, W2-table->minW2()); hminus -= numeric_limits::epsilon(); hplus -= numeric_limits::epsilon(); if (hminus < 0) hminus = 0; if (hplus < 0) hplus = 0; if (hplus == 0 && hminus == 0) { cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl; return 0; } double a = table->get(Q2, W2+hplus, t); if (a <= 0) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2+hplus) << ", pol=" << (pol == transverse ? 'T' : 'L') << endl; return 0; } double b = table->get(Q2, W2-hminus, t); if (b <= 0) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2-hminus) << ", pol=" << (pol == transverse ? 'T' : 'L') << endl; return 0; } derivative = (a-b)/(hplus+hminus); // // Finally calculate lambda // double jacobian = (W2-protonMass2+Q2)/value; lambda = jacobian*derivative; } // end fall back solution if (mSettings->verboseLevel() > 3) { cout << "CrossSection::logDerivateOfAmplitude(): "; if (lambdaFromTable) cout << "Info, lambda taken from table." << endl; else cout << "Info, lambda derived numerically from proton amplitude table" << endl; } // // Check lambda value. // At a lambda of ~0.6 both corrections have equal value // of around 2.9. This will yield excessive large (unphysical) // corrections. Large values are typically caused by fluctuations // and glitches in the tables and should be rare. // double maxLambda = mSettings->maxLambdaUsedInCorrections(); if (fabs(lambda) > maxLambda) { if (mSettings->verboseLevel() > 2) { cout << "CrossSection::logDerivateOfAmplitude(): "; cout << "Warning, lambda is excessively large (" << lambda << ") at " ; cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl; cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl; } lambda = lambda > 0 ? maxLambda : -maxLambda; } if (std::isinf(lambda)) { cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; cout << " Set to 0." << endl; return 0; } if (std::isnan(lambda)) { cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; cout << " Set to 0." << endl; return 0; } return lambda; } // // UPC only version // double CrossSection::logDerivateOfAmplitude(double t, double xpom) const { double lambda = 0; bool lambdaFromTable = true; Table *table = 0; if (!mProtonTableCollection) { cout << "CrossSection::logDerivateOfAmplitude(): no proton table defined to obtain lambda." << endl; cout << " Corrections not available. Should be off." << endl; return 0; } // // Usual numeric issues at boundaries. // Subtracting an eps in log(xpom) does the trick. // if (xpom > mProtonTableCollection->maxX()) { xpom = exp(log(xpom)-numeric_limits::epsilon()); } if (xpom < mProtonTableCollection->minX()) { xpom = exp(log(xpom)+numeric_limits::epsilon()); } // // If the lambda table is present we use the more accurate and numerically // stable table value. Otherwise we calculate it from the table(s). // lambda = mProtonTableCollection->get(xpom, t, lambda_real, table); if (!table) { // no lambda value from correct table, use fallback solution lambdaFromTable = false; (void) mProtonTableCollection->get(xpom, t, mean_A, table); // use obtained table from here on if (!table) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid pointer to lookup table." << endl; return 0; } // // Here's the tricky part (see comments in non-UPC version above) // double theLogxpom = log(xpom); double dlogxpom = table->binWidthX(); // assuming table is in log x ???????? FIX later double maxLogxpom = log(table->maxX()); double derivative; double hplus, hminus; hplus = hminus = 0.5*dlogxpom; hplus = min(hplus, fabs(maxLogxpom-theLogxpom)); hminus = min(hminus, fabs(theLogxpom-log(table->minX()))); hminus -= numeric_limits::epsilon(); hplus -= numeric_limits::epsilon(); if (hminus < 0) hminus = 0; if (hplus < 0) hplus = 0; if (hplus == 0 && hminus == 0) { cout << "CrossSection::logDerivateOfAmplitude(): Warning, cannot find derivative." << endl; return 0; } double a = table->get(exp(theLogxpom+hplus), t); if (a <= 0) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << a << '.' << endl; cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom+hplus)) << endl; return 0; } double b = table->get(exp(theLogxpom-hminus), t); if (a <= 0) { cout << "CrossSection::logDerivateOfAmplitude(): got invalid value from table, value=" << b << '.' << endl; cout << " t=" << t << ", W=" << sqrt(exp(theLogxpom-hminus)) << endl; return 0; } derivative = log(a/b)/(hplus+hminus); // // Finally calculate lambda // Directly dlog(A)/-dlog(xpom) here. // lambda = -derivative; } if (mSettings->verboseLevel() > 3) { cout << "CrossSection::logDerivateOfAmplitude(): "; if (lambdaFromTable) cout << "Info, lambda taken from table." << endl; else cout << "Info, lambda derived numerically from proton amplitude table" << endl; cout << " t=" << t << ", xpom=" << xpom << endl; } // // Check lambda value. // At a lambda of ~0.6 both corrections have equal value // of around 2.9. This will yield excessive large (unphysical) // corrections. Large values are typically caused by fluctuations // and glitches in the tables and should be rare. // double maxLambda = mSettings->maxLambdaUsedInCorrections(); if (fabs(lambda) > maxLambda) { if (mSettings->verboseLevel() > 2) { cout << "CrossSection::logDerivateOfAmplitude(): "; cout << "Warning, lambda is excessively large (" << lambda << ") at " ; cout << "xpom=" << xpom << ", t=" << t << endl; cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl; } lambda = lambda > 0 ? maxLambda : -maxLambda; } if (std::isinf(lambda)) { cout << "CrossSection::logDerivateOfAmplitude(): error, lambda = infinity for" << endl; cout << " t=" << t << ", xpom=" << xpom << endl; cout << " Set to 0." << endl; return 0; } if (std::isnan(lambda)) { cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for" << endl; cout << " t=" << t << ", xpom=" << xpom << endl; cout << " Set to 0." << endl; return 0; } return lambda; } double CrossSection::logDerivateOfGluonDensity(double t, double Q2, double W2, GammaPolarization pol) const { double lambda = 0; bool lambdaFromTable = true; Table *table = 0; if (!mProtonTableCollection) { cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl; cout << " Corrections not available. Should be off." << endl; return 0; } // // If the lambda table is present we use the more accurate and numerically // stable table value. Otherwise we calculate it from the table(s). // lambda = mProtonTableCollection->get(Q2, W2, t, pol, lambda_skew, table); if (!table) { // // no lambda value from correct table, use fallback solution // lambda_skew ~ lambda_real // lambda = logDerivateOfAmplitude(t, Q2, W2, pol) ; // This is an approximation in this case } // end fall back solution if (mSettings->verboseLevel() > 3) { cout << "CrossSection::logDerivateOfGluonDensity(): "; if (lambdaFromTable) cout << "Info, lambda taken from table." << endl; else cout << "Info, lambda taken from logDerivateOfAmplitude as an approximation" << endl; } // // Checking lambda. // At a lambda of ~0.6 both corrections have equal value // of around 2.9. This will yield excessive large (unphysical) // corrections. Large values are typically caused by fluctuations // and glitches in the tables and should be rare. // double maxLambda = mSettings->maxLambdaUsedInCorrections(); if (fabs(lambda) > maxLambda) { if (mSettings->verboseLevel() > 2) { cout << "CrossSection::logDerivateOfGluonDensity(): "; cout << "Warning, lambda is excessively large (" << lambda << ") at " ; cout << "Q2=" << Q2 << ", W2=" << W2 << ", t=" << t << endl; cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl; } lambda = lambda > 0 ? maxLambda : -maxLambda; } if (std::isinf(lambda)) { cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = inf for pol=" << (pol == transverse ? 'T' : 'L') << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; cout << " Set to 0." << endl; return 0; } if (std::isnan(lambda)) { cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; cout << " Set to 0." << endl; return 0; } return lambda; } // // UPC version // double CrossSection::logDerivateOfGluonDensity(double t, double xpom) const { double lambda = 0; bool lambdaFromTable = true; Table *table = 0; if (!mProtonTableCollection) { cout << "CrossSection::logDerivateOfGluonDensity(): no proton table defined to obtain lambda." << endl; cout << " Corrections not available. Should be off." << endl; return 0; } // // Usual numeric issues at boundaries. // Subtracting an eps in log(xpom) does the trick. // if (xpom > mProtonTableCollection->maxX()) { xpom = exp(log(xpom)-numeric_limits::epsilon()); } if (xpom < mProtonTableCollection->minX()) { xpom = exp(log(xpom)+numeric_limits::epsilon()); } // // If the lambda table is present we use the more accurate and numerically // stable table value. Otherwise we calculate it from the table(s). // lambda = mProtonTableCollection->get(xpom, t, lambda_skew, table); if (!table) { // // no lambda value from correct table, use fallback solution // lambda_skew ~ lambda_real // lambdaFromTable = false; lambda = logDerivateOfAmplitude(t, xpom); } if (mSettings->verboseLevel() > 3) { cout << "CrossSection::logDerivateOfGluonDensity(): "; if (lambdaFromTable) cout << "Info, lambda taken from table." << endl; else cout << "Info, lambda taken from logDerivateOfAmplitude as approximation." << endl; } // // Check lambda value. // At a lambda of ~0.6 both corrections have equal value // of around 2.9. This will yield excessive large (unphysical) // corrections. Large values are typically caused by fluctuations // and glitches in the tables and should be rare. // double maxLambda = mSettings->maxLambdaUsedInCorrections(); if (fabs(lambda) > maxLambda) { if (mSettings->verboseLevel() > 2) { cout << "CrossSection::logDerivateOfGluonDensity(): "; cout << "Warning, lambda is excessively large (" << lambda << ") at " ; cout << "xpom=" << xpom << ", t=" << t << endl; cout << "Set to " << (lambda > 0 ? maxLambda : -maxLambda) << "." << endl; } lambda = lambda > 0 ? maxLambda : -maxLambda; } if (std::isinf(lambda)) { cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda = infinity for" << endl; cout << " t=" << t << ", xpom=" << xpom << endl; cout << " Set to 0." << endl; return 0; } if (std::isnan(lambda)) { cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for" << endl; cout << " t=" << t << ", xpom=" << xpom << endl; cout << " Set to 0." << endl; return 0; } return lambda; } double CrossSection::realAmplitudeCorrection(double lambda) const { // // Correction factor for real amplitude contribution // double beta = tan(lambda*M_PI/2.); double correction = 1 + beta*beta; return correction; } double CrossSection::realAmplitudeCorrection(double t, double Q2, double W2, GammaPolarization pol) const { double lambda = logDerivateOfAmplitude(t, Q2, W2, pol); double correction = realAmplitudeCorrection(lambda); // correction *= exp(-10*Kinematics::x(Q2, W2)); // damped return correction; } // // UPC version // double CrossSection::realAmplitudeCorrection(double t, double xpom) const { double lambda = logDerivateOfAmplitude(t, xpom); double correction = realAmplitudeCorrection(lambda); return correction; } double CrossSection::skewednessCorrection(double lambda) const { // // Skewedness correction // double R = pow(2.,2*lambda+3)*TMath::Gamma(lambda+2.5)/(sqrt(M_PI)*TMath::Gamma(lambda+4)); double correction = R*R; return correction; } double CrossSection::skewednessCorrection(double t, double Q2, double W2, GammaPolarization pol) const { double lambda = logDerivateOfAmplitude(t, Q2, W2, pol); double correction = skewednessCorrection(lambda); // correction *= exp(-10*Kinematics::x(Q2, W2)); // damped return correction; } // // UPC version // double CrossSection::skewednessCorrection(double t, double xpom) const { double lambda = logDerivateOfAmplitude(t, xpom); double correction = skewednessCorrection(lambda); return correction; } Index: trunk/src/Sartre.cpp =================================================================== --- trunk/src/Sartre.cpp (revision 436) +++ trunk/src/Sartre.cpp (revision 437) @@ -1,1106 +1,1108 @@ //============================================================================== // Sartre.cpp // // Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich // // This file is part of Sartre. // // This program is free software: you can redistribute it and/or modify // it under the terms of the GNU General Public License as published by // the Free Software Foundation. // This program is distributed in the hope that it will be useful, // but without any warranty; without even the implied warranty of // merchantability or fitness for a particular purpose. See the // GNU General Public License for more details. // You should have received a copy of the GNU General Public License // along with this program. If not, see . // // Author: Thomas Ullrich // Last update: // $Date$ // $Author$ //============================================================================== // // Note: // When not using runcards, the user must first create an instance of Sartre // then get the settings via one of: // Sartre::runSettings() // EventGeneratorSettings::instance() // Once init() is called settings cannot be changed any more. //============================================================================== #include "Version.h" #include "Kinematics.h" #include "Sartre.h" #include "ModeFinderFunctor.h" #include "Math/BrentMinimizer1D.h" #include "Math/IntegratorMultiDim.h" #include "Math/GSLMinimizer.h" #include "Math/Functor.h" #include "TUnuranMultiContDist.h" #include #include #include #include "TH2D.h" #include "TFile.h" using namespace std; #define PR(x) cout << #x << " = " << (x) << endl; Sartre::Sartre() { mIsInitialized = false; mCurrentEvent = 0; mNucleus = 0; mUpcNucleus= 0; mSettings = 0; mPDF_Functor = 0; mPDF = 0; mEventCounter = 0; mTriesCounter = 0; mTotalCrossSection = 0; mCrossSection = 0; mTableCollection = 0; mProtonTableCollection = 0; mUnuran = 0; mEvents = 0; mTries = 0; mS = 0; mA = 0; } Sartre::~Sartre() { delete mNucleus; delete mUpcNucleus; delete mPDF_Functor; delete mPDF; delete mCrossSection; if (mTableCollection != mProtonTableCollection) { delete mTableCollection; delete mProtonTableCollection; } else delete mTableCollection; delete mUnuran; } bool Sartre::init(const char* runcard) { mStartTime = time(0); bool ok; // // Reset member variables. // Note that one instance of Sartre should be able to get // initialized multiple times. // mEvents = 0; mTries = 0; mTotalCrossSection = 0; // // Print header // string ctstr(ctime(&mStartTime)); ctstr.erase(ctstr.size()-1, 1); cout << "/========================================================================\\" << endl; cout << "| |" << endl; cout << "| Sartre, Version " << setw(54) << left << VERSION << right << '|' << endl; cout << "| |" << endl; cout << "| An event generator for exclusive diffractive vector meson production |" << endl; cout << "| in ep and eA collisions based on the dipole model. |" << endl; cout << "| |" << endl; cout << "| Copyright (C) 2010-2019 Tobias Toll and Thomas Ullrich |" << endl; cout << "| |" << endl; cout << "| This program is free software: you can redistribute it and/or modify |" << endl; cout << "| it under the terms of the GNU General Public License as published by |" << endl; cout << "| the Free Software Foundation, either version 3 of the License, or |" << endl; cout << "| any later version. |" << endl; cout << "| |" << endl; cout << "| Code compiled on " << setw(12) << left << __DATE__; cout << setw(41) << left << __TIME__ << right << '|' << endl; cout << "| Run started at " << setw(55) << left << ctstr.c_str() << right << '|' << endl; cout << "\\========================================================================/" << endl; mSettings = EventGeneratorSettings::instance(); // EventGeneratorSettings is a singleton // // Read runcard if available // if (runcard) { if (!mSettings->readSettingsFromFile(runcard)) { cout << "Error, reading runcard '" << runcard << "'. File doesn't exist or is not readable." << endl; exit(1); } else cout << "Runcard is '" << runcard << "'." << endl; } else cout << "No runcard provided." << endl; // // Set beam particles and center of mass energy // Note, that if we are in UPC mode the electron // is actually treated as the hadron beam, i.e. // the mass is the proton mass. // mElectronBeam = mSettings->eBeam(); mHadronBeam = mSettings->hBeam(); mS = Kinematics::s(mElectronBeam, mHadronBeam); mA = mSettings->A(); bool allowBreakup = mSettings->enableNuclearBreakup(); if (mA == 1) allowBreakup = false; if (allowBreakup) { if (!getenv("SARTRE_DIR")) { cout << "Error, environment variable 'SARTRE_DIR' is not defined. It is required\n" "to locate tables needed for the generation if nuclear breakups." << endl; exit(1); } } if (mNucleus) delete mNucleus; mNucleus = new FrangibleNucleus(mA, allowBreakup); if (mSettings->UPC()) { if (mUpcNucleus) delete mUpcNucleus; mUpcNucleus = new FrangibleNucleus(mSettings->UPCA(), allowBreakup); } string upcNucleusName; if (mSettings->UPC()) { upcNucleusName = mUpcNucleus->name(); cout << "Sartre is running in UPC mode" << endl; cout << "Hadron 1 beam species: " << mNucleus->name() << " (" << mA << ")" << endl; cout << "Hadron 1 beam: " << mHadronBeam << endl; cout << "Hadron 2 beam species: " << upcNucleusName << " (" << mSettings->UPCA() << ")" << endl; cout << "Hadron 2 beam: " << mElectronBeam << endl; } else { cout << "Hadron beam species: " << mNucleus->name() << " (" << mA << ")" << endl; cout << "Hadron beam: " << mHadronBeam << endl; cout << "Electron beam: " << mElectronBeam << endl; } // // Get details about the processes and models // mDipoleModelType = mSettings->dipoleModelType(); mDipoleModelParameterSet = mSettings->dipoleModelParameterSet(); mVmID = mSettings->vectorMesonId(); cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl; cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl; cout << "Process is "; if (mSettings->UPC()) { cout << mNucleus->name() << " + " << upcNucleusName << " -> " << mNucleus->name() << "' + " << upcNucleusName << "' + "; } else { if (mA > 1) cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + "; else cout << "e + p -> e' + p' + "; } TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID); cout << vectorMesonPDG->GetName() << endl; // // Print-out seed for reference // cout << "Random generator seed: " << mSettings->seed() << endl; // // Load in the tables containing the amplitude moments // if (!getenv("SARTRE_DIR")) { cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl; exit(1); } if (mTableCollection) delete mTableCollection; mTableCollection = new TableCollection; ok = mTableCollection->init(mA, mDipoleModelType, mDipoleModelParameterSet, mVmID); if (!ok) { cout << "Error, could not initialize lookup tables for requested process." << endl; return false; } // // Load in the p tables for the lambda lookup tables (or to calculate lambda if not available) // if (mSettings->correctForRealAmplitude() || mSettings->correctSkewedness()) { if (mA == 1) { mProtonTableCollection = mTableCollection; } else { if (mProtonTableCollection) delete mProtonTableCollection; mProtonTableCollection = new TableCollection; ok = mProtonTableCollection->init(1, mDipoleModelType, mDipoleModelParameterSet, mVmID); if (!ok) { cout << "Error: could not initialize proton lookup tables for requested process." << endl; cout << " These tables are needed for skewedness and real amplitude corrections." << endl; return false; } } } else mProtonTableCollection = 0; // // Kinematic limits and generator range // // There are 3 ranges we have to deal with // 1. the kinematic range requested by the user // if given. // The user can only control Q2 and W but not t. // For UPC that's xpom. // 2. the range of the table(s) // 3. the kinematically allowed range // // Of course (3) is more complex than a simple cube/square. // However, we deal with the detailed shape of the kinematic // range later using Kinematics::valid() when we generate the // individual events. // For setting up UNU.RAN we have to get the cubic/square // envelope that satifies (1)-(3). // Note, that they are correlated which makes the order // in which we do things a bit tricky. // // // Step 1: // Set the limits to that of the table(s). // Note, the indices 0-2 refer to t, Q2, and W2. // For UPC we have only t and xpom. // if (mSettings->UPC()) { mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT(); mLowerLimit[1] = mTableCollection->minX(); mUpperLimit[1] = mTableCollection->maxX(); mLowerLimit[2] = mUpperLimit[2] = 0; } else { mLowerLimit[0] = mTableCollection->minT(); mUpperLimit[0] = mTableCollection->maxT(); mLowerLimit[1] = mTableCollection->minQ2(); mUpperLimit[1] = mTableCollection->maxQ2(); mLowerLimit[2] = mTableCollection->minW2(); mUpperLimit[2] = mTableCollection->maxW2(); } // // Step 2: // Kinematic limits might overrule boundaries from step 1 // if (mSettings->UPC()) { double kineXpomMin = Kinematics::xpomMin(vectorMesonPDG->Mass(), mUpperLimit[0], mHadronBeam, mElectronBeam); double kineXpomMax = 1; mLowerLimit[1] = max(kineXpomMin, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineXpomMax); double kineTmax = Kinematics::tmax(kineXpomMin); double kineTmin = Kinematics::tmin(mHadronBeam.E()); mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax); } else { // double kineYmax = Kinematics::ymax(mS, vectorMesonPDG->Mass()); double kineYmin = Kinematics::ymin(mS, vectorMesonPDG->Mass()); double kineQ2min = Kinematics::Q2min(kineYmin); double kineQ2max = Kinematics::Q2max(mS); double kineW2min = Kinematics::W2min(vectorMesonPDG->Mass()); double kineW2max = Kinematics::W2max(mS); kineQ2min = max(kineQ2min, mLowerLimit[1]); kineQ2max = min(kineQ2max, mUpperLimit[1]); kineW2min = max(kineW2min, mLowerLimit[2]); kineW2max = min(kineW2max, mUpperLimit[2]); double kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // first arg (t) set to 0 (recursive) double kineTmax = Kinematics::tmax(kineXPmin); double kineTmin = Kinematics::tmin(mHadronBeam.E()); mLowerLimit[0] = max(kineTmin, mLowerLimit[0]); mUpperLimit[0] = min(mUpperLimit[0], kineTmax); mLowerLimit[1] = max(kineQ2min, mLowerLimit[1]); mUpperLimit[1] = min(mUpperLimit[1], kineQ2max); mLowerLimit[2] = max(kineW2min, mLowerLimit[2]); mUpperLimit[2] = min(mUpperLimit[2], kineW2max); } // // Step 3: // Deal with user provided limits. // User settings are ignored (switched off) if min >= max. // if (mSettings->UPC()) { if (mSettings->xpomMin() < mSettings->xpomMax()) { if (mSettings->xpomMin() < mLowerLimit[1]) { cout << "Warning, requested lower limit of xpomeron (" << mSettings->xpomMin() << ") " << "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). "; cout << "Limit has no effect." << endl; } else { mLowerLimit[1] = mSettings->xpomMin(); } if (mSettings->xpomMax() > mUpperLimit[1]) { cout << "Warning, requested upper limit of xpomeron (" << mSettings->xpomMax() << ") " << "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). "; cout << "Limit has no effect." << endl; } else { mUpperLimit[1] = mSettings->xpomMax(); } } } else { if (mSettings->W2min() < mSettings->W2max()) { // W2 first if (mSettings->W2min() < mLowerLimit[2]) { cout << "Warning, requested lower limit of W (" << mSettings->Wmin() << ") " << "is smaller than limit given by lookup tables and/or kinematic range (" << sqrt(mLowerLimit[2]) << "). "; cout << "Limit has no effect." << endl; } else { mLowerLimit[2] = mSettings->W2min(); } if (mSettings->W2max() > mUpperLimit[2]) { cout << "Warning, requested upper limit of W (" << mSettings->Wmax() << ") " << "exceeds limit given by lookup tables and/or kinematic range (" << sqrt(mUpperLimit[2]) << "). "; cout << "Limit has no effect." << endl; } else { mUpperLimit[2] = mSettings->W2max(); } } if (mSettings->Q2min() < mSettings->Q2max()) { // Q2 if (mSettings->Q2min() < mLowerLimit[1]) { cout << "Warning, requested lower limit of Q2 (" << mSettings->Q2min() << ") " << "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[1] << "). "; cout << "Limit has no effect." << endl; } else { mLowerLimit[1] = mSettings->Q2min(); // ????? kineXPmin = Kinematics::xpomeron(0, kineQ2min, kineW2max, vectorMesonPDG->Mass()); // new Q2min changes tmax // ????? mUpperLimit[0] = min(mUpperLimit[0], Kinematics::tmax(kineXPmin)); } if (mSettings->Q2max() > mUpperLimit[1]) { cout << "Warning, requested upper limit of Q2 (" << mSettings->Q2max() << ") " << "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[1] << "). "; cout << "Limit has no effect." << endl; } else { mUpperLimit[1] = mSettings->Q2max(); } } } // // Check if any phase space is left // if (mLowerLimit[0] >= mUpperLimit[0]) { cout << "Invalid range in t: t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl; exit(1); } if (mLowerLimit[1] >= mUpperLimit[1]) { if (mSettings->UPC()) cout << "Invalid range in xpomeron: xpomeron=["; else cout << "Invalid range in Q2: Q2=["; cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl; exit(1); } if (!mSettings->UPC() && mLowerLimit[2] >= mUpperLimit[2]) { cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl; exit(1); } // // Print-out limits (all verbose levels) // if (mSettings->verbose()) { cout << "Kinematic limits used for event generation:" << endl; if (mSettings->UPC()) { cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl; cout << setw(10) << "xpom=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl; } else { cout << setw(10) << " t=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]" << endl; cout << setw(10) << "Q2=[" << mLowerLimit[1] << ", " << mUpperLimit[1] << "]" << endl; cout << setw(10) << " W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]" << endl; } } // // Skewedness and real part of amplitude corrections rely // on parameters that are stored in tables. Here we check // if they exist and if the tables have the right size. // This is mostly to inform the user what's going on // for verbose level 2 and higher. // // Should they not exist or are too small they are calculated // on the fly (see class CrossSection for details). // If this happens we also check if the kinematic range // of the proton amplitude tables is big enough to allow // calculating lambda for the actual kinematic range. // If not we stop execution. // if (mProtonTableCollection) { if (mSettings->UPC()) { bool skewUPC_TableExists = false; bool realUPC_TableExists = false; bool skewUPC_TableFits = false; bool realUPC_TableFits = false; bool ampUPC_TableFits = false; skewUPC_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_skew); realUPC_TableExists = mProtonTableCollection->tableExists(transverse, lambda_real); if (skewUPC_TableExists) skewUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew); if (realUPC_TableExists) realUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real); ampUPC_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A); if (mSettings->verbose() && mSettings->verboseLevel() > 1) { if (!realUPC_TableExists) { cout << "Info: there is no table with lambda values to perform the real" << endl; cout << " amplitude corrections. The missing lambda values will be" << endl; cout << " calculated on the fly from the ep amplitude tables." << endl; } if (!skewUPC_TableExists) { cout << "Info: there is no table with lambda values to perform the skewedness" << endl; cout << " corrections. We will use the lambda values from the real amplitude" << endl; cout << " corrections instead. This is a good approximation." << endl; } if (realUPC_TableExists && !realUPC_TableFits) { cout << "Info: the table with the lambda values to perform the real part corrections" << endl; cout << " does not have sufficient coverage of the current kinematic range." << endl; cout << " The missing lambda values will be calculated on the fly from the ep" << endl; cout << " amplitude tables." << endl; } if (skewUPC_TableExists && !skewUPC_TableFits) { cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl; cout << " does not have sufficient coverage of the current kinematic range." << endl; cout << " The missing lambda values will be replaced by the lambda values used" << endl; cout << " for the real amplitude corrections. This is a good approximation." << endl; } } if (!ampUPC_TableFits && (!realUPC_TableExists || (realUPC_TableExists && !realUPC_TableFits))) { cout << "Error: due to missing correction tables or tables that do not cover the" << endl; cout << " full kinematic range the needed lambda values have to be calculated" << endl; cout << " on the fly. However, the ep mean amplitude table that is needed for" << endl; cout << " this calculation does not cover the full kinematic range selected." << endl; cout << " Corrections need to be switched off to run with the existing tables" << endl; cout << " and the selected kinematic range." << endl; exit(1); } } else { bool skewT_TableExists, skewL_TableExists; bool realT_TableExists, realL_TableExists; bool skewT_TableFits, skewL_TableFits; bool realT_TableFits, realL_TableFits; bool ampT_TableFits, ampL_TableFits; skewT_TableExists = mProtonTableCollection->tableExists(transverse, lambda_skew); skewL_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_skew); realT_TableExists = mProtonTableCollection->tableExists(transverse, lambda_real); realL_TableExists = mProtonTableCollection->tableExists(longitudinal, lambda_real); if (skewT_TableExists) skewT_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew, transverse); if (skewL_TableExists) skewL_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_skew, longitudinal); if (realT_TableExists) realT_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real, transverse); if (realL_TableExists) realL_TableFits = tableFitsKinematicRange(mProtonTableCollection, lambda_real, longitudinal); ampT_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A, transverse); ampL_TableFits = tableFitsKinematicRange(mProtonTableCollection, mean_A, longitudinal); if (mSettings->verbose() && mSettings->verboseLevel() > 1) { if (!realT_TableExists) { cout << "Info: there is no table with lambda values to perform the real" << endl; cout << " amplitude corrections for tranverse polarized photons. The" << endl; cout << " missing lambda values will be calculated on the fly from" << endl; cout << " the referring ep amplitude tables." << endl; } if (!realL_TableExists) { cout << "Info: there is no table with lambda values to perform the real" << endl; cout << " amplitude corrections for longitudinal polarized photons. " << endl; cout << " The missing lambda values will be calculated on the fly " << endl; cout << " from the referring ep amplitude tables." << endl; } if (!skewT_TableExists) { cout << "Info: there is no table with lambda values to perform the skewedness" << endl; cout << " corrections for tranverse polarized photons. We will use the" << endl; cout << " lambda values from the real amplitude corrections instead. This" << endl; cout << " is a good approximation." << endl; } if (!skewL_TableExists) { cout << "Info: there is no table with lambda values to perform the skewedness" << endl; cout << " corrections for longitudinal polarized photons. We will use the" << endl; cout << " lambda values from the real amplitude corrections instead. This" << endl; cout << " is a good approximation." << endl; } if (realT_TableExists && !realT_TableFits) { cout << "Info: the table with the lambda values to perform the real part corrections" << endl; cout << " for tranverse polarized photons does not have sufficient coverage of" << endl; cout << " the current kinematic range. The missing lambda values will be calculated" << endl; cout << " on the fly from the referring ep amplitude tables." << endl; } if (realL_TableExists && !realL_TableFits) { cout << "Info: the table with the lambda values to perform the real part corrections" << endl; cout << " for longitudinal polarized photons does not have sufficient coverage of" << endl; cout << " the current kinematic range. The missing lambda values will be calculated" << endl; cout << " on the fly from the referring ep amplitude tables." << endl; } if (skewT_TableExists && !skewT_TableFits) { cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl; cout << " for tranverse polarized photons does not have sufficient coverage for " << endl; cout << " the current kinematic range. The missing lambda values will be replaced " << endl; cout << " by the lambda values used for the real amplitude corrections. This is a" << endl; cout << " good approximation." << endl; } if (skewL_TableExists && !skewL_TableFits) { cout << "Info: the table with the lambda values to perform the skewedness corrections" << endl; cout << " for longitudinal polarized photons does not have sufficient coverage for " << endl; cout << " the current kinematic range. The missing lambda values will be replaced " << endl; cout << " by the lambda values used for the real amplitude corrections. This is a" << endl; cout << " good approximation." << endl; } } if ( (!ampT_TableFits && (!realT_TableExists || (realT_TableExists && !realT_TableFits))) || (!ampL_TableFits && (!realL_TableExists || (realL_TableExists && !realL_TableFits)))) { cout << "Error: due to missing correction tables or tables that do not cover the" << endl; cout << " full kinematic range the needed lambda values have to be calculated" << endl; cout << " on the fly. However, the ep mean amplitude table that is needed for" << endl; cout << " this calculation does not cover the full kinematic range selected." << endl; cout << " Corrections need to be switched off to run with the existing tables" << endl; cout << " and the selected kinematic range." << endl; exit(1); } } } // // Setup cross-section functor // It is this functor that is used by all other functors, // functions, and wrappers when dealing with cross-sections. // if (mCrossSection) delete mCrossSection; mCrossSection = new CrossSection; mCrossSection->setTableCollection(mTableCollection); if (mProtonTableCollection) mCrossSection->setProtonTableCollection(mProtonTableCollection); // // UNU.RAN needs the domain (boundaries) and the mode. // The domain is already defined, here we find the mode, which is tricky. // The max. cross-section is clearly at the domain boundary in Q2=Q2min. // The position in W2 and t is not obvious. It sits along a line given // by tmax(W2). The approach here is to use the BrentMinimizer1D that // performs first a scan a then a Brent fit. // double theMode[3]; if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl; if (mSettings->UPC()) { // // For now we keep the option to create a histogram of the // generated cross-sections(t, log(xpom)) for QA purposes. // if (mSettings->verbose() && mSettings->verboseLevel() > 9) { cout << "Creating 2D histogram of cross-section as fct of t and log(xpom)." << endl; cout << "Be patient this might take a while. Histo stored in file 'landscape.root'" << endl; TFile file("landscape.root", "RECREATE"); int nbins_t = 300; int nbins_logx = 300; TH2D histo("histo", "mode finder studies", nbins_t, mLowerLimit[0], mUpperLimit[0], // t nbins_logx, log(mLowerLimit[1]), log(mUpperLimit[1])); // log(xpom) for (int ix = 1; ix <= nbins_logx; ix++) { for (int it = 1; it <= nbins_t; it++) { double logx = histo.GetYaxis()->GetBinCenter(ix); double x = exp(logx); double t = histo.GetXaxis()->GetBinCenter(it); if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy(), t, x, vectorMesonPDG->Mass(), false)) { histo.SetBinContent(it, ix, -5.e9); } else { histo.SetBinContent(it, ix, (*mCrossSection)(t, x)); } } } histo.Write(); file.Close(); cout << "Histogram written to file. All done." << endl; } // // Create functor // Mode is at the max t available. // Note that the functop works with log(xpom). // theMode[0] = mUpperLimit[0]; // t theMode[1] = mUpperLimit[1]; // xpom for UPC theMode[2] = 0; // not used for UPC UPCModeFinderFunctor modeFunctor(mCrossSection, vectorMesonPDG->Mass(), mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy()); ROOT::Math::BrentMinimizer1D minimizer; minimizer.SetFunction(modeFunctor, log(mLowerLimit[1]), log(mUpperLimit[1])); minimizer.SetNpx(100000); ok = minimizer.Minimize(1000000, 1.e-8, 1.e-10); if (! ok) { cout << "Error, failed to find mode of pdf." << endl; return false; } // // Get the result // theMode[1] = exp(minimizer.XMinimum()); // xpom theMode[0] = Kinematics::tmax(theMode[1]); double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1]); if (mSettings->verbose()) { cout << "\tlocation: t=" << theMode[0] << ", xpom=" << theMode[1] << "; value: " << crossSectionAtMode << endl; } // // Consistency check of mode // if (crossSectionAtMode <= 0) { cout << "Error: cross-section at mode value is invalid." << endl; return false; } if (!Kinematics::validUPC(mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy(), theMode[0], theMode[1], vectorMesonPDG->Mass(), true)) { cout << "Error: mode of pdf is outside kinematic limits." << endl; return false; } } // for EIC else { theMode[0] = mUpperLimit[0]; // t theMode[1] = mLowerLimit[1]; // Q2 (xpom for UPC) theMode[2] = mLowerLimit[2]; // W2 (dummy for UPC) ModeFinderFunctor modeFunctor(mCrossSection, theMode[1], vectorMesonPDG->Mass(), mLowerLimit[0], mUpperLimit[0]); ROOT::Math::BrentMinimizer1D minimizer; minimizer.SetFunction(modeFunctor, mLowerLimit[2], mUpperLimit[2]); minimizer.SetNpx(static_cast(mUpperLimit[2]-mLowerLimit[2])); ok = minimizer.Minimize(100000, 0, 1.e-8); if (! ok) { cout << "Error, failed to find mode of pdf." << endl; exit(1); } theMode[2] = minimizer.XMinimum(); // W2 theMode[0] = Kinematics::tmax(0, theMode[1], theMode[2], vectorMesonPDG->Mass()); // first arg (t) must be 0 here if (theMode[0] > mUpperLimit[0]) theMode[0] = mUpperLimit[0]; double crossSectionAtMode = (*mCrossSection)(theMode[0], theMode[1], theMode[2]); if (mSettings->verbose()) { cout << "\tlocation: t=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]); cout << "; value: " << crossSectionAtMode << endl; } } // // Initialize 2D (UPC) or 3D random generator // // Test show that UNU.RAN runs smoother in log(Q2) // and log(cross-section). Functor CrossSection has // a spezialized method for UNU.RAN, unuranPDF(). // // In UPC mode the mPDF_Functor is using a different // method and is only 2D since Q2=0 // // domain and mode for Q2 -> log(Q2) or xpom -> log(xpom) mLowerLimit[1] = log(mLowerLimit[1]); mUpperLimit[1] = log(mUpperLimit[1]); theMode[1] = log(theMode[1]); if (mPDF_Functor) delete mPDF_Functor; if (mPDF) delete mPDF; if (mSettings->UPC()) mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 2); else mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &CrossSection::unuranPDF, 3); mPDF = new TUnuranMultiContDist(*mPDF_Functor, true); // last arg = pdf in log or not mPDF->SetDomain(mLowerLimit, mUpperLimit); mPDF->SetMode(theMode); if (mUnuran) delete mUnuran; mUnuran = new TUnuran; mCrossSection->setCheckKinematics(false); // avoid numeric glitch in Init() mUnuran->Init(*mPDF, "method=hitro"); mCrossSection->setCheckKinematics(true); mUnuran->SetSeed(mSettings->seed()); // // Burn in generator // double xrandom[3]; for (int i=0; i<1000; i++) { mUnuran->SampleMulti(xrandom); } mEventCounter = 0; mTriesCounter = 0; mIsInitialized = true; cout << "Sartre is initialized." << endl << endl; return true; } bool Sartre::init(const string& str) // overloaded version of init() { if (str.empty()) return init(); else return init(str.c_str()); } vector > Sartre::kinematicLimits() { vector > array; array.push_back(make_pair(mLowerLimit[0], mUpperLimit[0])); // t array.push_back(make_pair(exp(mLowerLimit[1]), exp(mUpperLimit[1]))); // Q2 or xpom if (!mSettings->UPC()) array.push_back(make_pair(sqrt(mLowerLimit[2]), sqrt(mUpperLimit[2]))); // W return array; } Event* Sartre::generateEvent() { if (!mIsInitialized) { cout << "Sartre::generateEvent(): Error, Sartre is not initialized yet." << endl; cout << " Call init() before trying to generate events." << endl; return 0; } double xrandom[3]; TParticlePDG *vectorMesonPDG = mSettings->lookupPDG(mVmID); double vmMass = vectorMesonPDG->Mass(); // // Generate one event // while (true) { mTriesCounter++; delete mCurrentEvent; mCurrentEvent = new Event; // // Get t, Q2, W2 from TUnuran and check for kinematics. // Q2 is generated as log(Q2) so we transform it back first. // This is the only place where Kinematics::valid() is called // with the fully correct xpomeron calculation switched on. // mUnuran->SampleMulti(xrandom); xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2 or log(xpom) -> xpom bool isValidEvent; if (mSettings->UPC()) { isValidEvent = Kinematics::validUPC(mSettings->hadronBeamEnergy(), mSettings->electronBeamEnergy(), xrandom[0], xrandom[1], vectorMesonPDG->Mass(), (mSettings->verboseLevel() > 1)); } else { isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], vmMass, true, (mSettings->verboseLevel() > 1)); } if (!isValidEvent) { if (mSettings->verboseLevel() > 2) cout << "Sartre::generateEvent(): event rejected, not within kinematic limits" << endl; continue; } // // Fill beam particles in Event structure // Kinematics for eA is reported as 'per nucleon' // if (mSettings->UPC()) { // // For UPC some of the event variables // are filled in the final state generator. // They are set to 0 here. // mCurrentEvent->eventNumber = mEventCounter; mCurrentEvent->t = xrandom[0]; // t mCurrentEvent->Q2 = 0; mCurrentEvent->x = 0; mCurrentEvent->y = 0; mCurrentEvent->s = mS; // s mCurrentEvent->W = 0; mCurrentEvent->beta = 1; mCurrentEvent->xpom = xrandom[1]; mCurrentEvent->polarization = GammaPolarization::transverse; mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall(); + mCurrentEvent->crossSectionRatioLT = mCrossSection->crossSectionRatioLTOfLastCall(); Particle eIn, hIn; eIn.index = 0; eIn.pdgId = mUpcNucleus->pdgID(); // misuse ebeam spot for this eIn.status = 1; eIn.p = mElectronBeam; hIn.index = 1; hIn.pdgId = mNucleus->pdgID(); hIn.status = 1; hIn.p = mHadronBeam; mCurrentEvent->particles.push_back(eIn); mCurrentEvent->particles.push_back(hIn); } else { mCurrentEvent->eventNumber = mEventCounter; mCurrentEvent->t = xrandom[0]; // t mCurrentEvent->Q2 = xrandom[1]; // Q2 mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y mCurrentEvent->s = mS; // s mCurrentEvent->W = sqrt(xrandom[2]); mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall(); mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall(); + mCurrentEvent->crossSectionRatioLT = mCrossSection->crossSectionRatioLTOfLastCall(); Particle eIn, hIn; eIn.index = 0; eIn.pdgId = 11; // e- eIn.status = 1; eIn.p = mElectronBeam; hIn.index = 1; hIn.pdgId = mNucleus->pdgID(); hIn.status = 1; hIn.p = mHadronBeam; mCurrentEvent->particles.push_back(eIn); mCurrentEvent->particles.push_back(hIn); } // // Generate the final state particles // bool ok; if (mSettings->UPC()) { // also fills some of the undefined event variables ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->xpom, (mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent); } else { ok = mFinalStateGenerator.generate(mVmID, mCurrentEvent->t, mCurrentEvent->y, mCurrentEvent->Q2, (mCurrentEvent->diffractiveMode == incoherent), mA, mCurrentEvent); } if (!ok) { if (mSettings->verboseLevel() > 1) cout << "Sartre::generateEvent(): failed to generate final state" << endl; continue; } break; } mEventCounter++; // // Nuclear breakup // // If the event is incoherent the final state generator does produce a // 'virtual' proton with m > m_p which is used in Nucleus to calculate // the excitation energy and the boost. // int indexOfScatteredHadron = 6; bool allowBreakup = mSettings->enableNuclearBreakup(); if (mA == 1) allowBreakup = false; if (mNucleus) mNucleus->resetBreakup(); // clear previous event in any case if (allowBreakup && mCurrentEvent->diffractiveMode == incoherent && mNucleus) { int nFragments = mNucleus->breakup(mCurrentEvent->particles[indexOfScatteredHadron].p); // // Merge the list of products into the event list. // We loose some information here. The user can always go back to // the nucleus and check the decay products for more details. // In the original list the energy is per nuclei, here we transform it // to per nucleon to stay consistent with Sartre conventions. // const vector& products = mNucleus->breakupProducts(); for (int i=0; iparticles.size(); fragment.pdgId = products[i].pdgId; fragment.status = 1; fragment.p = products[i].p*(1/static_cast(products[i].A)); fragment.parents.push_back(indexOfScatteredHadron); mCurrentEvent->particles.push_back(fragment); } } // // Complete event record // if (!mSettings->UPC()) { mCurrentEvent->xpom = Kinematics::xpomeron(mCurrentEvent->t, mCurrentEvent->Q2, mCurrentEvent->W*mCurrentEvent->W, vmMass); mCurrentEvent->beta = mCurrentEvent->x/mCurrentEvent->xpom; } return mCurrentEvent; } double Sartre::totalCrossSection() { if (mTotalCrossSection == 0) { // // Limits of integration in t, Q2, W2 // or in the case of UPC t, xpom, dummy // double xmin[3]; double xmax[3]; copy(mLowerLimit, mLowerLimit+3, xmin); copy(mUpperLimit, mUpperLimit+3, xmax); // // At this point mLowerLimit[1] and mUpperLimit[1] // are in log(Q2) or log(xpom). // xmin[1] = exp(xmin[1]); // log Q2 limit -> Q2 limit or equivalent xpom for UPC xmax[1] = exp(xmax[1]); // log Q2 limit -> Q2 limit mTotalCrossSection = calculateTotalCrossSection(xmin, xmax); } return mTotalCrossSection; } double Sartre::totalCrossSection(double lower[3], double upper[3]) // t, Q2, W (or t, xpom, dummy for UPC) { lower[2] *= lower[2]; upper[2] *= upper[2]; // W -> W2 double result = calculateTotalCrossSection(lower, upper); return result; } EventGeneratorSettings* Sartre::runSettings() { return EventGeneratorSettings::instance(); } double Sartre::calculateTotalCrossSection(double lower[3], double upper[3]) { double result = 0; if (!mIsInitialized) { cout << "Sartre::calculateTotalCrossSection(): Error, Sartre is not initialized yet." << endl; cout << " Call init() before trying to generate events." << endl; return result; } // // Calculate integral using adaptive numerical method // // Options: ADAPTIVE, kVEGAS, kPLAIN, kMISER // no abs tolerance given -> relative only const double precision = 5e-4; ROOT::Math::Functor wfL((*mCrossSection), mSettings->UPC() ? 2 : 3); ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE, 0, precision, 1000000); ig.SetFunction(wfL); result = ig.Integral(lower, upper); // // If it fails we switch to a MC integration which is usually more robust // although not as accurate. This should happen very rarely if at all. // if (result <= numeric_limits::epsilon()) { cout << "Sartre::calculateTotalCrossSection(): warning, adaptive integration failed - switching to VEGAS method." << endl; ROOT::Math::IntegratorMultiDim igAlt(ROOT::Math::IntegrationMultiDim::kVEGAS); igAlt.SetFunction(wfL); igAlt.SetRelTolerance(precision); igAlt.SetAbsTolerance(0); result = igAlt.Integral(lower, upper); } // // UPC with symmetric beams requires the cross-section // being multiplied by a factor of 2. // if (mSettings->UPC()) { if (mSettings->A() == mSettings->UPCA()) result *= 2; } return result; } const FrangibleNucleus* Sartre::nucleus() const {return mNucleus;} void Sartre::listStatus(ostream& os) const { os << "Event summary: " << mEventCounter<< " events generated, " << mTriesCounter << " tried" << endl; time_t delta = runTime(); os << "Total time used: " << delta/60 << " min " << delta - 60*(delta/60) << " sec" << endl; } time_t Sartre::runTime() const { time_t now = time(0); return now-mStartTime; } bool Sartre::tableFitsKinematicRange(TableCollection* tc, AmplitudeMoment mom, GammaPolarization pol) { if (!tc) return false; bool fitsIn = tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), pol, mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), pol, mom); return fitsIn; } // UPC version bool Sartre::tableFitsKinematicRange(TableCollection* tc, AmplitudeMoment mom) { if (!tc) return false; bool fitsIn = tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) && tc->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mom) && tc->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mom); return fitsIn; } //============================================================================== // // Utility functions and operators (helpers) // //============================================================================== ostream& operator<<(ostream& os, const TLorentzVector& v) { os << v.Px() << '\t' << v.Py() << '\t' << v.Pz() << '\t' << v.E() << '\t'; double m2 = v*v; if (m2 < 0) os << '(' << -sqrt(-m2) << ')'; else os << '(' << sqrt(m2) << ')'; return os; }