Index: trunk/src/TableCollection.cpp =================================================================== --- trunk/src/TableCollection.cpp (revision 373) +++ trunk/src/TableCollection.cpp (revision 374) @@ -1,569 +1,569 @@ //============================================================================== // TableCollection.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 that we do not take the lambda_A table into account when calculating +// Note that we do not take the lambda_real table into account when calculating // the range since there is a fall back solution to calculate lambda if the // table is not present. See class CrossSection. // //============================================================================== #include "TSystemDirectory.h" #include "TSystem.h" #include "TList.h" #include "EventGeneratorSettings.h" #include "TableCollection.h" #include "Table.h" #include #include #include #include #include #define PR(x) cout << #x << " = " << (x) << endl; TableCollection::TableCollection() {/* no op */} TableCollection::TableCollection(int A, DipoleModelType typ, DipoleModelParameterSet set, int vmID) { init(A, typ, set, vmID); } TableCollection& TableCollection::operator=(const TableCollection& tc) { if (this != &tc) { for (unsigned int i=0; iWorkingDirectory(); // // Build directory path // stringstream pathstream; pathstream << getenv("SARTRE_DIR") << "/tables/" << A << '/'; if (type == bSat) pathstream << "bSat"; else if (type == bNonSat) pathstream << "bNonSat"; else pathstream << "bCGC"; pathstream << '/' << vmID; string path = pathstream.str(); // // Query list of all files in directory // and create tables for each file ending // in ".root", ignore others. // TSystemDirectory directory; directory.SetDirectory(path.c_str()); TList *list = directory.GetListOfFiles(); if (!list) { cout << "TableCollection::init(): Error, cannot find directory '" << path.c_str() << "' holding tables." << endl; return false; } TIter next(list); unsigned int numberOfTablesRead = 0; bool upcMode = EventGeneratorSettings::instance()->UPC(); while (TSystemFile* file = dynamic_cast(next()) ) { if (file->IsDirectory()) continue; // ignore directories string name(file->GetName()); size_t pos = name.find(".root"); if (pos == string::npos || name.substr(pos) != string(".root")) continue; // ignore files not ending in .root string fullpath = path + '/' + name; Table *table = new Table; if (table->read(fullpath.c_str()) && table->dipoleModelParameterSet() == set) { if ( (!upcMode && !table->isUPC()) || (upcMode && table->isUPC())) mTables.push_back(table); numberOfTablesRead++; if (EventGeneratorSettings::instance()->verboseLevel() > 1) cout << "Loaded table from file '" << fullpath.c_str() << "'." << endl; } } // // Cleanup // // Change ROOT directory back to directory we were before reading the tables. // Otherwise reading tables interferes with user application. // list->Delete(); gSystem->ChangeDirectory(saveCWD.c_str()); if (!numberOfTablesRead) { cout << "TableCollection::init(): Error, could not find any tables at '" << path.c_str() << "'" << endl; cout << " that fit the requested parameters: A=" << A << ", type=" << type << ", set=" << set << " , vmID=" << vmID << endl; return false; } return true; } bool TableCollection::tableExists(GammaPolarization pol, AmplitudeMoment mom) const { Table* currentTable; for (unsigned int i=0; ipolarization() != pol) continue; if (currentTable->moment() != mom) continue; return true; } return false; } // // UPC version // bool TableCollection::tableExists(AmplitudeMoment mom) const { Table* currentTable; for (unsigned int i=0; imoment() != mom) continue; return true; } return false; } bool TableCollection::available(double Q2, double W2, double t, GammaPolarization pol, AmplitudeMoment mom) const { // // Check if table can provide this value // unsigned short nTables = 0; Table* currentTable; for (unsigned int i=0; ipolarization() != pol) continue; if (currentTable->moment() != mom) continue; if (t >= currentTable->minT() && t <= currentTable->maxT()) { if (Q2 >= currentTable->minQ2() && Q2 <= currentTable->maxQ2()) { if (W2 >= currentTable->minW2() && W2 <= currentTable->maxW2()) { nTables++; } } } } if (nTables) return true; else return false; } bool TableCollection::available(double xpom, double t, AmplitudeMoment mom) const { // // Check if table can provide this value // unsigned short nTables = 0; Table* currentTable; for (unsigned int i=0; imoment() != mom) continue; if (t >= currentTable->minT() && t <= currentTable->maxT()) { if (xpom >= currentTable->minX() && xpom <= currentTable->maxX()) { nTables++; } } } if (nTables) return true; else return false; } double TableCollection::get(double Q2, double W2, double t, GammaPolarization pol, AmplitudeMoment mom) const { Table *table; return get(Q2, W2, t, pol, mom, table); } // // UPC version // double TableCollection::get(double xpom, double t, AmplitudeMoment mom) const { Table *table; return get(xpom, t, mom, table); } double TableCollection::get(double Q2, double W2, double t, GammaPolarization pol, AmplitudeMoment mom, Table *&table) const { static unsigned int errorCount = 0; const unsigned int maxErrorCount = 10; // // First get the tables that contain the necessary info. // Later this should be a bit refined, here we simply // loop over all tables to collect the relevant one(s). // vector associatedTables; Table* currentTable; for (unsigned int i=0; ipolarization() != pol) continue; if (currentTable->moment() != mom) continue; if (t >= currentTable->minT() && t <= currentTable->maxT()) { if (Q2 >= currentTable->minQ2() && Q2 <= currentTable->maxQ2()) { if (W2 >= currentTable->minW2() && W2 <= currentTable->maxW2()) { associatedTables.push_back(currentTable); } } } } if (associatedTables.size() == 0) { table = 0; - if (mom != lambda_A && mom != lambda_skew) { // no warnings needed (can be calculated w/o tables) + if (mom != lambda_real && mom != lambda_skew) { // no warnings needed (can be calculated w/o tables) string txt; if (mom == mean_A) txt = "mean_A"; else if (mom == mean_A2) txt = "mean_A2"; else if (mom == variance_A) txt = "variance_A"; else txt = "unknown"; cout << "TableCollection::get(): Warning, could not find any table containing t=" << t << ", Q2=" << Q2 << ", W2=" << W2 << endl; cout << " Tables searched were for moment = " << txt.c_str() << ", polarization = " << (pol == transverse ? 'T' : 'L') << endl; errorCount++; if (errorCount > maxErrorCount) { cout << "TableCollection::get(): Error: Too many warnings (>" << maxErrorCount << "), possibly due to missing table(s)." << endl; cout << " Stop execution now. Please check the installation of tables." << endl; exit(1); } } return 0; } // // In case of overlap of tables the following // policy applies: // 1. Use the table with the highest priority. // 2. If there's more than one high priority table // we average their values (if > 0). // unsigned int maxPriority = 0; for (unsigned int i=0; ipriority() > maxPriority) maxPriority = associatedTables[i]->priority(); } double result = 0; int validCounter = 0; table = 0; for (unsigned int i=0; ipriority() == maxPriority) { double value = associatedTables[i]->get(Q2, W2, t); if (value > 0) { validCounter++; result += value; table = associatedTables[i]; } } } if (validCounter) result /= validCounter; return result; } // // UPC version // double TableCollection::get(double xpom, double t, AmplitudeMoment mom, Table *&table) const { static unsigned int errorCount = 0; const unsigned int maxErrorCount = 10; // // First get the tables that contain the necessary info. // Later this should be a bit refined, here we simply // loop over all tables to collect the relevant one(s). // vector associatedTables; Table* currentTable; for (unsigned int i=0; imoment() != mom) continue; if (t >= currentTable->minT() && t <= currentTable->maxT()) { if (xpom >= currentTable->minX() && xpom <= currentTable->maxX()) { associatedTables.push_back(currentTable); } } } if (associatedTables.size() == 0) { table = 0; - if ( mom != lambda_skew && mom != lambda_A ) { // no warnings needed (can be calculated w/o tables) + if ( mom != lambda_skew && mom != lambda_real ) { // no warnings needed (can be calculated w/o tables) string txt; if (mom == mean_A) txt = "mean_A"; else if (mom == mean_A2) txt = "mean_A2"; else if (mom == variance_A) txt = "variance_A"; else txt = "unknown"; cout << "TableCollection::get(): Warning, could not find any table containing t=" << t << ", xp=" << xpom << ". Moment = " << txt << ", " << mom << endl; errorCount++; if (errorCount > maxErrorCount) { cout << "TableCollection::get(): Error: Too many warnings (>" << maxErrorCount << "), possibly due to missing table(s)." << endl; cout << " Stop execution now. Please check the installation of tables." << endl; exit(1); } } return 0; } // // In case of overlap of tables the following // policy applies: // 1. Use the table with the highest priority. // 2. If there's more than one high priority table // we average their values (if > 0). // unsigned int maxPriority = 0; for (unsigned int i=0; ipriority() > maxPriority) maxPriority = associatedTables[i]->priority(); } double result = 0; int validCounter = 0; table = 0; double value = 0; for (unsigned int i=0; ipriority() == maxPriority) { value = associatedTables[i]->get(xpom, t); if (value > 0) { validCounter++; result += value; table = associatedTables[i]; } } } if (validCounter) result /= validCounter; return result; } void TableCollection::list(ostream& os, bool opt) const { for (unsigned int i=0; ilist(os, opt); } double TableCollection::minQ2() const { if (EventGeneratorSettings::instance()->UPC()) return 0; return minimumValue(0); } double TableCollection::maxQ2() const { if (EventGeneratorSettings::instance()->UPC()) return 0; return maximumValue(0); } double TableCollection::minW2() const { if (EventGeneratorSettings::instance()->UPC()) return 0; return minimumValue(1); } double TableCollection::maxW2() const { if (EventGeneratorSettings::instance()->UPC()) return 0; return maximumValue(1); } double TableCollection::minW() const {return sqrt(minW2());} double TableCollection::maxW() const {return sqrt(maxW2());} double TableCollection::minT() const { if (EventGeneratorSettings::instance()->UPC()) return minimumValue(1); else return minimumValue(2); } double TableCollection::maxT() const { if (EventGeneratorSettings::instance()->UPC()) return maximumValue(1); else return maximumValue(2); } double TableCollection::minX() const { if (EventGeneratorSettings::instance()->UPC()) return minimumValue(0); else return 0; } double TableCollection::maxX() const { if (EventGeneratorSettings::instance()->UPC()) return maximumValue(0); else return 0; } // // For regular tables: kind: Q2=0, W2=1, T=2 // For UPC tables: kind: xpom = 0, t=1 // double TableCollection::minimumValue(unsigned int kind) const { double minPerTableType[4]; // L, L2, T, T2 fill(minPerTableType, minPerTableType+4, numeric_limits::max()); for (unsigned int i=0; iUPC()) { switch (kind) { case (0): val = mTables[i]->minX(); break; case (1): val = mTables[i]->minT(); break; } } else { switch (kind) { case (0): val = mTables[i]->minQ2(); break; case (1): val = mTables[i]->minW2(); break; default: val = mTables[i]->minT(); break; } } if (mTables[i]->isLongitudinal()) { // L or L2 if (mTables[i]->isMeanA()) minPerTableType[0] = min(minPerTableType[0],val); // L else minPerTableType[1] = min(minPerTableType[1],val); // L2 } else { // T or T2 if (mTables[i]->isMeanA()) minPerTableType[2] = min(minPerTableType[2],val); // T else minPerTableType[3] = min(minPerTableType[3],val); // T2 } } int startElement = EventGeneratorSettings::instance()->UPC() ? 2 : 0; double largestMin = *max_element(minPerTableType+startElement, minPerTableType+4); return largestMin; } double TableCollection::maximumValue(unsigned int kind) const { double maxPerTableType[4]; // L, L2, T, T2 fill(maxPerTableType, maxPerTableType+4, -numeric_limits::max()); for (unsigned int i=0; iUPC()) { switch (kind) { case (0): val = mTables[i]->maxX(); break; case (1): val = mTables[i]->maxT(); break; } } else { switch (kind) { case (0): val = mTables[i]->maxQ2(); break; case (1): val = mTables[i]->maxW2(); break; default: val = mTables[i]->maxT(); break; } } if (mTables[i]->isLongitudinal()) { // L or L2 if (mTables[i]->isMeanA()) maxPerTableType[0] = max(maxPerTableType[0],val); // L else maxPerTableType[1] = max(maxPerTableType[1],val); // L2 } else { // T or T2 if (mTables[i]->isMeanA()) maxPerTableType[2] = max(maxPerTableType[2],val); // T else maxPerTableType[3] = max(maxPerTableType[3],val); // T2 } } int startElement = EventGeneratorSettings::instance()->UPC() ? 2 : 0; double smallestMax = *min_element(maxPerTableType+startElement, maxPerTableType+4); return smallestMax; } Index: trunk/src/Sartre.h =================================================================== --- trunk/src/Sartre.h (revision 373) +++ trunk/src/Sartre.h (revision 374) @@ -1,111 +1,113 @@ //============================================================================== // Sartre.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$ //============================================================================== #ifndef Sartre_h #define Sartre_h #include "Event.h" #include "EventGeneratorSettings.h" #include "ExclusiveFinalStateGenerator.h" #include "CrossSection.h" #include "TableCollection.h" #include "FrangibleNucleus.h" #include "Enumerations.h" #include "TLorentzVector.h" #include "Math/Functor.h" #include "TUnuran.h" #include #include #include using namespace std; class TUnuranMultiContDist; class Sartre { public: Sartre(); virtual ~Sartre(); virtual bool init(const char* = 0); virtual bool init(const string&); virtual Event* generateEvent(); virtual double totalCrossSection(); // in kinematic limits used for generation virtual double totalCrossSection(double lower[3], double upper[3]); // t, Q2, W EventGeneratorSettings* runSettings(); const FrangibleNucleus* nucleus() const; void listStatus(ostream& os=cout) const; time_t runTime() const; vector > kinematicLimits(); // t, Q2, W private: virtual double calculateTotalCrossSection(double lower[3], double upper[3]); // t, Q2, W2 - + virtual bool tableFitsKinematicRange(TableCollection*, AmplitudeMoment, GammaPolarization); + virtual bool tableFitsKinematicRange(TableCollection*, AmplitudeMoment); // UPC version + private: Sartre(const Sartre&); Sartre operator=(const Sartre&); bool mIsInitialized; time_t mStartTime; unsigned long mEvents; unsigned long mTries; double mTotalCrossSection; TLorentzVector mElectronBeam; TLorentzVector mHadronBeam; double mS; unsigned int mA; int mVmID; DipoleModelType mDipoleModelType; DipoleModelParameterSet mDipoleModelParameterSet; Event *mCurrentEvent; FrangibleNucleus *mNucleus; FrangibleNucleus *mUpcNucleus; TableCollection *mTableCollection; TableCollection *mProtonTableCollection; CrossSection *mCrossSection; EventGeneratorSettings *mSettings; double mLowerLimit[3]; // t, Q2, W2 (t, xpom for UPC) double mUpperLimit[3]; // t, Q2, W2 ROOT::Math::Functor *mPDF_Functor; TUnuran *mUnuran; TUnuranMultiContDist *mPDF; unsigned long mEventCounter; unsigned long mTriesCounter; ExclusiveFinalStateGenerator mFinalStateGenerator; }; ostream& operator<<(ostream& os, const TLorentzVector&); #endif Index: trunk/src/Enumerations.h =================================================================== --- trunk/src/Enumerations.h (revision 373) +++ trunk/src/Enumerations.h (revision 374) @@ -1,33 +1,33 @@ //============================================================================== // Enumerations.h // // Copyright (C) 2010-2017 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$ //============================================================================== #ifndef Enumerations_h #define Enumerations_h enum DipoleModelType {bSat, bNonSat, bCGC}; enum GammaPolarization {transverse, longitudinal}; -enum AmplitudeMoment {mean_A, mean_A2, lambda_A, variance_A, lambda_skew}; +enum AmplitudeMoment {mean_A, mean_A2, variance_A, lambda_real, lambda_skew}; enum DiffractiveMode {coherent, incoherent}; enum DipoleModelParameterSet {KMW, HMPZ, CUSTOM}; enum TableSetType {total_and_coherent, coherent_and_incoherent}; #endif Index: trunk/src/CrossSection.cpp =================================================================== --- trunk/src/CrossSection.cpp (revision 373) +++ trunk/src/CrossSection.cpp (revision 374) @@ -1,911 +1,921 @@ //============================================================================== // CrossSection.cpp // // Copyright (C) 2010-2018 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; } 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; // // 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); + if (mSettings->correctForRealAmplitude()) { + lambda = logDerivateOfAmplitude(t, xpom); result *= realAmplitudeCorrection(lambda); - } - if (mSettings->correctSkewedness()){ - lambda = logDerivateOfGluonDensity(t, xpom); + } + 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 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); + 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); + 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); + 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;} 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_A, table); + 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 time per point. + // 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. - // In all cases where something goes wrong we set lambda to - // its most probably value (0.2) - // + // - const double lambdaInCaseOfError = 0.2; - - 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 (lambda < 0) { - if (mSettings->verboseLevel() > 2) { - cout << "CrossSection::logDerivateOfAmplitude(): "; - cout << "Warning, lambda is negative (" << lambda << "). " ; - cout << "Set to " << (-lambda) << "." << endl; - } - lambda = -lambda; - } - */ 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; - lambda = lambdaInCaseOfError; - } - 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; - lambda = lambdaInCaseOfError; - } - 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; + cout << " Set to 0." << 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; - } - // - // 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. - // In all cases where something goes wrong we set lambda to - // its most probably value (0.2) - // - - const double lambdaInCaseOfError = 0.2; - - - 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; - lambda = lambdaInCaseOfError; - } if (std::isnan(lambda)) { - cout << "CrossSection::logDerivateOfGluonDensity(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl; + cout << "CrossSection::logDerivateOfAmplitude(): error, lambda is NaN for pol=" << (pol == transverse ? 'T' : 'L') << endl; cout << " t=" << t << ", Q2=" << Q2 << ", W=" << sqrt(W2) << endl; - lambda = lambdaInCaseOfError; + 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()); + 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_A, table); - + 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 + (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); - derivative = log(a/b)/(hplus+hminus); - + 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; + cout << "Info, lambda taken from table." << endl; else - cout << "Info, lambda derived numerically from proton amplitude table" << endl; - cout << " t=" << t << ", xpom=" << xpom << endl; + 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. - // In all cases where something goes wrong we set lambda to - // its most probably value (0.2) // - const double lambdaInCaseOfError = 0.2; - 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; - lambda = lambdaInCaseOfError; + 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; - lambda = lambdaInCaseOfError; + 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); + + lambda = logDerivateOfAmplitude(t, xpom); } if (mSettings->verboseLevel() > 3) { cout << "CrossSection::logDerivateOfGluonDensity(): "; if (lambdaFromTable) - cout << "Info, lambda taken from table." << endl; + cout << "Info, lambda taken from table." << endl; else - cout << "Info, lambda taken from logDerivateOfAmplitude as approximation." << endl; + 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. - // In all cases where something goes wrong we set lambda to - // its most probably value (0.2) // - const double lambdaInCaseOfError = 0.2; - 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; - lambda = lambdaInCaseOfError; + 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; - lambda = lambdaInCaseOfError; + 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 373) +++ trunk/src/Sartre.cpp (revision 374) @@ -1,1012 +1,1105 @@ //============================================================================== // 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; + 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; + 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 << "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) + // 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; } } // - // Check if the lambda table covers the kinematic range - // or if we have to calculate lambda from the p mean_A table - // directly. - // This is just to inform the user, no action is taken. - // Here we assume that transverse and longitudinal tables - // always have the same range and vheck only one. - // - // If even the p mean_A table is not big enough we switch - // corrections off and inform the user. + // 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) { - bool lambdaTableExists, tooSmallLambda, tooSmallMean; - if (mSettings->UPC()) { - lambdaTableExists = mProtonTableCollection->tableExists(lambda_A); - - tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), lambda_A); - - tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mean_A); - + 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 { - lambdaTableExists = mProtonTableCollection->tableExists(transverse, lambda_A); + 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; - tooSmallLambda = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, lambda_A); - - tooSmallMean = !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mLowerLimit[1]+numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mLowerLimit[0]+numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mLowerLimit[2]+numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A) || - !mProtonTableCollection->available(mUpperLimit[1]-numeric_limits::epsilon(), mUpperLimit[2]-numeric_limits::epsilon(), mUpperLimit[0]-numeric_limits::epsilon(), transverse, mean_A); - } - - if (!lambdaTableExists || tooSmallLambda) { - - if (tooSmallMean) { - - cout << "Warning: the kinematic coverage of the table containing the lambda values " << endl; - cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl; - cout << " then the required range (or does not exist). The proton amplitude " << endl; - cout << " tables that could be used to calculate them on the fly are too" << endl; - cout << " small as well. Corrections will be switched off." << endl; - - mSettings->setCorrectForRealAmplitude(false); - mSettings->setCorrectSkewedness(false); - delete mProtonTableCollection; - mProtonTableCollection = 0; - } - else { - if (lambdaTableExists) { - cout << "Info: the kinematic coverage of the table containing the lambda values " << endl; - cout << " needed for skewedness and/or real amplitude corrections is smaller" << endl; - cout << " then the required range. The missing lambda values will be " << endl; - cout << " calculated using the proton amplitude tables directly." << endl; - } - else { - cout << "Info: the table containing the lambda values needed for skewedness and/or real " << endl; - cout << " amplitude corrections is not available. The missing lambda values will be" << endl; - cout << " calculated using the referring proton amplitude tables directly." << endl; + 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; } } 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(); 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(); 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; }