No OneTemporary
Actions

Size

50 KB

Subscribers

None

View Options

	Index: trunk/src/InclusiveFinalStateGenerator.h
	===================================================================
	--- trunk/src/InclusiveFinalStateGenerator.h (revision 0)
	+++ trunk/src/InclusiveFinalStateGenerator.h (revision 490)
	@@ -0,0 +1,45 @@
	+//==============================================================================
	+// ExclusiveFinalStateGenerator.h
	+//
	+// Copyright (C) 2024 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 <http://www.gnu.org/licenses/>.
	+//
	+// Author: Tobias Toll
	+// Last update:
	+// $Date: $
	+// $Author: ullrich $
	+//==============================================================================
	+#ifndef InclusiveFinalStateGenerator_h
	+#define InclusiveFinalStateGenerator_h
	+#include "Pythia8/Pythia.h"
	+#include "FinalStateGenerator.h"
	+
	+class InclusiveFinalStateGenerator : public FinalStateGenerator {
	+public:
	+ InclusiveFinalStateGenerator();
	+ ~InclusiveFinalStateGenerator();
	+
	+ bool generate(int A, Event *event);
	+
	+ double uiGamma_N(double* var, double* par) const;
	+ double gamma_N(unsigned int i, double val);
	+
	+private:
	+ Pythia8::Pythia mPythia;
	+ Pythia8::Event *mPythiaEvent;
	+ Pythia8::ParticleData *mPdt;
	+
	+
	+};
	+#endif
	Index: trunk/src/InclusiveFinalStateGenerator.cpp
	===================================================================
	--- trunk/src/InclusiveFinalStateGenerator.cpp (revision 0)
	+++ trunk/src/InclusiveFinalStateGenerator.cpp (revision 490)
	@@ -0,0 +1,440 @@
	+//==============================================================================
	+// ExclusiveFinalStateGenerator.cpp
	+//
	+// Copyright (C) 2024 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 <http://www.gnu.org/licenses/>.
	+//
	+// Author: Tobias Toll
	+// Last update:
	+// $Date: $
	+// $Author: ullrich $
	+//==============================================================================
	+#include "InclusiveFinalStateGenerator.h"
	+#include "EventGeneratorSettings.h"
	+#include "Kinematics.h"
	+#include "Event.h"
	+#include "Math/BrentRootFinder.h"
	+#include "Math/GSLRootFinder.h"
	+#include "Math/RootFinderAlgorithms.h"
	+#include "Math/IFunction.h"
	+#include "Math/SpecFuncMathCore.h"
	+#include <cmath>
	+#include <algorithm>
	+#include <limits>
	+#include <iomanip>
	+#include "Math/WrappedTF1.h"
	+#include "Math/GaussIntegrator.h"
	+#include "TF1.h"
	+
	+
	+#define PR(x) cout << #x << " = " << (x) << endl;
	+
	+//-------------------------------------------------------------------------------
	+//
	+// Implementation of ExclusiveFinalStateGenerator
	+//
	+//-------------------------------------------------------------------------------
	+
	+InclusiveFinalStateGenerator::InclusiveFinalStateGenerator()
	+{
	+ //
	+ // Setup Pythia8 which is need to hadronize the dipole (and higher Fock states)
	+ //
	+ mPythiaEvent = &(mPythia.event);
	+ mPdt = &(mPythia.particleData);
	+ mPythia.readString("ProcessLevel:all = off"); // needed
	+ mPythia.readString("Check:event = off");
	+ if (!mPythia.init()) {
	+ cout << "Error initializing Pythia8. Stop here." << endl;
	+ exit(1);
	+ }
	+ else {
	+ cout << "Pythia8 initialized (" << PYTHIA_VERSION_INTEGER << ")" << endl;
	+ }
	+}
	+
	+InclusiveFinalStateGenerator::~InclusiveFinalStateGenerator() {/* no op */}
	+
	+double InclusiveFinalStateGenerator::uiGamma_N(double* var, double* par) const {
	+ double t=var[0];
	+ double s=par[0];
	+ double result = pow(t, s-1)*exp(-t);
	+ return result;
	+}
	+
	+double InclusiveFinalStateGenerator::gamma_N(unsigned int i, double val){
	+ TF1 fGamma_N("fGamma_N", this, &InclusiveFinalStateGenerator::uiGamma_N, -10, 10, 1);
	+ fGamma_N.SetParameter(0, i);
	+ ROOT::Math::WrappedTF1 wfGamma_N(fGamma_N);
	+ ROOT::Math::GaussIntegrator giGamma_N;
	+ giGamma_N.SetFunction(wfGamma_N);
	+ giGamma_N.SetAbsTolerance(0.);
	+ giGamma_N.SetRelTolerance(1e-5);
	+
	+ int sign=1;
	+ double low=0, up=val;
	+ if(val<0){
	+ low=val;
	+ up=0;
	+ sign=-1;
	+ }
	+ double Tfact=1;
	+ for(int j=i-1; j>0; j--){
	+ Tfact*=j;
	+ }
	+ return sign*giGamma_N.Integral(low, up)/Tfact;
	+}
	+
	+bool InclusiveFinalStateGenerator::generate(int A, Event *event)
	+{
	+ //
	+ // Get generator settings and the random generator
	+ //
	+ EventGeneratorSettings *settings = EventGeneratorSettings::instance();
	+ TRandom3 *rndm = settings->randomGenerator();
	+
	+ //
	+ // The beam particles must be present in the event list
	+ //
	+ int ePos = -1;
	+ int hPos = -1;
	+ bool parentsOK = true;
	+ if (event->particles.size() == 2) {
	+ if (abs(event->particles[0].pdgId) == 11) {
	+ ePos = 0;
	+ hPos = 1;
	+ }
	+ else if (abs(event->particles[1].pdgId) == 11) {
	+ ePos = 1;
	+ hPos = 0;
	+ }
	+ else
	+ parentsOK = false;
	+ }
	+ else
	+ parentsOK = false;
	+
	+ if (!parentsOK) {
	+ cout << "ExclusiveFinalStateGenerator::generate(): error, no beam particles in event list." << endl;
	+ return false;
	+ }
	+
	+ //
	+ // Store arguments locally
	+ // (Some could also be obtained from the event structure)
	+ //
	+ mA = A;
	+ double xpom=event->xpom;
	+ mT = Kinematics::tmax(xpom);//event->t;
	+ if (mT > 0) mT = -mT; // ensure t<0
	+ mQ2 = event->Q2;
	+ mY = event->y;
	+ mElectronBeam = event->particles[ePos].p;
	+ mHadronBeam = event->particles[hPos].p;
	+ mMX = event->MX;
	+ mS = Kinematics::s(mElectronBeam, mHadronBeam);
	+ double MX=event->MX;
	+ double MX2=MX*MX;
	+ //
	+ // Constants
	+ //
	+ double const twopi = 2*M_PI;
	+ double const hMass2 = mHadronBeam.M2();
	+ mMY2 = hMass2;
	+
	+
	+ //
	+ // Re-engineer scattered electron
	+ //
	+ // e'=(E', pt', pz') -> 3 unknowns
	+ //
	+ // Three equations:
	+ // 1: meme=E'E'-pt'pt'-pz'pz'
	+ // 2: Q2=-(e-e')^2=-2meme + 2(EE'-pz*pz')
	+ // 3: W2=(P+e-e')^2=mp2+2me2+2(EpE-Pzpz)-2(EpE'-Pzpz')-2(EE'-pzpz')
	+ //
	+
	+ double Ee=mElectronBeam.E();
	+ double Pe=mElectronBeam.Pz();
	+ double Ep=mHadronBeam.E();
	+ double Pp=mHadronBeam.Pz();
	+ double W=event->W;
	+ double W2=W*W;
	+ double Q2=event->Q2;
	+ // Take masses from the beams in case they are not actually electrons or protons
	+ double me2=mElectronBeam.M2();
	+ double mp2=mHadronBeam.M2();
	+ //
	+ // What we want for each particle:
	+ //
	+ double E, pz, pt, px, py, phi;
	+
	+ //
	+ // Equations 2 and 3 yield:
	+ //
	+ E = Pe(W2-mp2-2EeEp) + (Pp+Pe)Q2 + 2PePePp + 2me2*Pp;
	+ E /= 2(EePp-Ep*Pe);
	+ pz = Ee(W2-mp2) + (Ep+Ee)Q2 + 2EePePp + 2Epme2 - 2EeEeEp;
	+ pz /= 2(EePp-Ep*Pe);
	+ //
	+ // Equation 1:
	+ //
	+ pt = sqrt(EE-pzpz-me2);
	+ phi = rndm->Uniform(twopi);
	+ TLorentzVector theScatteredElectron(ptsin(phi), ptcos(phi), pz, E);
	+ //
	+ // Re-engineer virtual photon
	+ //
	+ // gamma=E-E'
	+ E=mElectronBeam.E()-theScatteredElectron.E();
	+ pz=mElectronBeam.Pz()-theScatteredElectron.Pz();
	+ px=mElectronBeam.Px()-theScatteredElectron.Px();
	+ py=mElectronBeam.Py()-theScatteredElectron.Py();
	+ TLorentzVector theVirtualPhoton = TLorentzVector(px, py, pz, E);
	+
	+ //
	+ // Re-engineer scattered proton/dissociated proton
	+ //
	+ E=(1-xpom)mHadronBeam.Pz()mHadronBeam.Pz()-mT/2.+0.5hMass2+0.5mMY2;
	+ E=E/mHadronBeam.E();
	+ pz = (1-xpom)*mHadronBeam.Pz();
	+ pt = sqrt(EE-pzpz-mMY2);
	+ if(pt!=pt){
	+ if(settings->verboseLevel() > 2)
	+ cout<<"No phase-space for scattered proton pt"<<endl;
	+ return false;
	+ }
	+ px = pt*cos(phi);
	+ py = pt*sin(phi);
	+ TLorentzVector theScatteredProton(px, py, pz, E);
	+ //
	+ // Use scattered proton to set up four momentum of pomeron:
	+ //
	+ TLorentzVector thePomeron=mHadronBeam-theScatteredProton;
	+
	+
	+ //
	+ // Finally the diffractive final state,
	+ // treat at first as a pseudo particle
	+ //
	+ TLorentzVector theXparticle((mHadronBeam + mElectronBeam) - (theScatteredElectron + theScatteredProton));
	+ //
	+ //The quark and anti quark
	+ //
	+ double mf=event->quarkMass;
	+
	+ double z=event->z;
	+ double pt2=z(1-z)Q2-mf*mf;
	+ if(4(mfmf+pt2)/MX2>1) pt2=MX2/4.-mf*mf;
	+ //generate angle
	+ phi = rndm->Uniform(twopi);
	+ double pxq=sqrt(pt2)*sin(phi);
	+ double pyq=sqrt(pt2)*cos(phi);
	+ double pxqbar=-pxq;
	+ double pyqbar=-pyq;
	+ double z0=0.5(1-sqrt(1-4(mf*mf+pt2)/MX2));
	+ if(mElectronBeam.Pz()<0 and z<0.5) z0=1-z0;
	+ if(mElectronBeam.Pz()>0 and z>0.5) z0=1-z0;
	+
	+ double qPplus=z0*(theXparticle.E()+theXparticle.Pz());
	+ double qPminus=(1-z0)*(theXparticle.E()-theXparticle.Pz());
	+ double qBarPplus=(1-z0)*(theXparticle.E()+theXparticle.Pz());
	+ double qBarPminus=z0*(theXparticle.E()-theXparticle.Pz());
	+
	+ //
	+ // Generate decorrelation between the qqbar system.
	+ //
	+
	+ //
	+ // First find which twist we are operating at
	+ //
	+ int Twist=0;
	+ double Omega=event->Omega;
	+ double eps=1e-3;
	+ for(int i=0; i<100; i++){
	+ double gammaN=gamma_N(i+1, -Omega/2);
	+ if(fabs(2exp(-Omega/2)gammaN)<eps){
	+ Twist=i+1;
	+ break;
	+ }
	+ }
	+ //
	+ // Now generate 2*Twist number of gluons:
	+ // Do they need a pz as well??
	+ //
	+ double eps2=z(1-z)(Q2+MX2);
	+ double r_average=1./sqrt(eps2);
	+ double Qs=sqrt(2)*Omega/r_average;
	+ vector<double> pom_px, pom_py, pom_pz;
	+ double pom_px_tot=thePomeron.Px();
	+ double pom_py_tot=thePomeron.Py();
	+ double pom_pz_tot=thePomeron.Pz();
	+ for(int i=0; i<2*Twist-1; i++){
	+ //#TT This part does not seem quite right yet.
	+ double p_x=rndm->Gaus(0., Qs);
	+ double p_y=rndm->Gaus(0., Qs);
	+ double p_z=rndm->Gaus(0., Qs);
	+ pom_px.push_back(thePomeron.Px()+p_x);
	+ pom_py.push_back(thePomeron.Py()+p_y);
	+ pom_pz.push_back(thePomeron.Pz()+p_z);
	+ pom_px_tot+=px;
	+ pom_py_tot+=py;
	+ pom_pz_tot+=pz;
	+ }
	+ //conserve the total pomeron momentum:
	+ pom_px.push_back(thePomeron.Px()-pom_px_tot);
	+ pom_py.push_back(thePomeron.Py()-pom_py_tot);
	+ pom_pz.push_back(thePomeron.Pz()-pom_pz_tot);
	+ //Distribute the recoils to the quark and antiquark:
	+ double Eq=(qPplus+qPminus)/2;
	+ double Eqbar=(qBarPplus+qBarPminus)/2;
	+ double pzq=(qPplus-qPminus)/2;
	+ double pzqbar=(qBarPplus-qBarPminus)/2;
	+
	+ for(unsigned int i=0; i<pom_px.size(); i++){
	+ //first calculate the relative velocity between q, qbar and pomeron
	+ double pomP=sqrt(pom_px[i]pom_px[i]+pom_py[i]pom_py[i]+pom_pz[i]*pom_pz[i]);
	+ double qP=sqrt(pxqpxq+pyqpyq+pzq*pzq);
	+ double qbarP=sqrt(pxqbarpxqbar+pyqbarpyqbar+pzqbar*pzqbar);
	+ double deltaVq=sqrt(1+qPqP/Eq/Eq-2(pom_px[i]pxq+pom_py[i]pyq+pom_pz[i]*pzq)/Eq/pomP);
	+ double deltaVqbar=sqrt(1+qbarPqbarP/Eqbar/Eqbar-2(pom_px[i]pxqbar+pom_py[i]pyqbar+pom_pz[i]*pzqbar)/Eqbar/pomP);
	+ //The momentum kick is inversely proportional to relative velocity
	+ double fraction=0.5;
	+ if(deltaVq+deltaVqbar!=0)
	+ fraction=deltaVqbar/(deltaVq+deltaVqbar);
	+ //Update momenta
	+ pxq+=fraction*pom_px[i];
	+ pyq+=fraction*pom_py[i];
	+ pxqbar+=(1-fraction)*pom_px[i];
	+ pyqbar+=(1-fraction)*pom_py[i];
	+ //Update the energies
	+ Eq=sqrt(pxqpxq+pyqpyq+pzqpzq+mfmf);
	+ Eqbar=sqrt(pxqbarpxqbar+pyqbarpyqbar+pzqbarpzqbar+mfmf);
	+ }
	+ //Fill the 4-vectors:
	+
	+ //the Quark:
	+ TLorentzVector theQuark(pxq, pyq, pzq, Eq);
	+
	+ //the Anti-Quark
	+ TLorentzVector theAntiQuark(pxqbar, pyqbar, pzqbar, Eqbar);
	+
	+ theXparticle=theQuark+theAntiQuark;
	+
	+ //
	+ // Add particles to event record
	+ //
	+ event->particles.resize(2+7);
	+ unsigned int eOut = 2;
	+ unsigned int gamma = 3;
	+ unsigned int X = 4;
	+ unsigned int pomeron = 5;
	+ unsigned int hOut = 6;
	+ unsigned int quark = 7;
	+ unsigned int antiquark = 8;
	+
	+ // Global indices
	+ event->particles[eOut].index = eOut;
	+ event->particles[gamma].index = gamma;
	+ event->particles[X].index = X;
	+ event->particles[pomeron].index = pomeron;
	+ event->particles[hOut].index = hOut;
	+ event->particles[quark].index = quark;
	+ event->particles[antiquark].index = antiquark;
	+
	+ // 4-vectors
	+ event->particles[eOut].p = theScatteredElectron;
	+ event->particles[hOut].p = theScatteredProton;
	+ event->particles[gamma].p = theVirtualPhoton;
	+ event->particles[X].p = theXparticle;
	+ event->particles[pomeron].p = thePomeron;
	+ event->particles[quark].p = theQuark;
	+ event->particles[antiquark].p = theAntiQuark;
	+
	+ // PDG Ids
	+ event->particles[eOut].pdgId = event->particles[ePos].pdgId; // same as incoming
	+ event->particles[hOut].pdgId = event->particles[hPos].pdgId; // same as incoming (breakup happens somewhere else)
	+ event->particles[gamma].pdgId = 22;
	+ event->particles[X].pdgId = 90; //pseudoparticle
	+ int iquark=event->quarkSpecies;
	+ if(iquark<=1)
	+ event->particles[pomeron].pdgId = 113;
	+ if(iquark==2)
	+ event->particles[pomeron].pdgId = 333;
	+ if(iquark==3)
	+ event->particles[pomeron].pdgId = 443;
	+ if(iquark==4)
	+ event->particles[pomeron].pdgId = 553;
	+ event->particles[quark].pdgId = iquark+1;
	+ event->particles[antiquark].pdgId = -event->particles[quark].pdgId;
	+
	+ // status
	+ //
	+ // HepMC conventions (February 2009).
	+ // 0 : an empty entry, with no meaningful information
	+ // 1 : a final-state particle, i.e. a particle that is not decayed further by
	+ // the generator (may also include unstable particles that are to be decayed later);
	+ // 2 : a decayed hadron or tau or mu lepton
	+ // 3 : a documentation entry (not used in PYTHIA);
	+ // 4 : an incoming beam particle;
	+ // 11 - 200 : an intermediate (decayed/branched/...) particle that does not
	+ // fulfill the criteria of status code 2
	+
	+ event->particles[ePos].status = 4;
	+ event->particles[hPos].status = 4;
	+ event->particles[eOut].status = 1;
	+ event->particles[hOut].status = mIsIncoherent ? 2 : 1;
	+ event->particles[gamma].status = 2;
	+ event->particles[X].status = 90;
	+ event->particles[pomeron].status = 2;
	+ event->particles[quark].status = 1;
	+ event->particles[antiquark].status = 1;
	+
	+ // parents (ignore dipole)
	+ event->particles[eOut].parents.push_back(ePos);
	+ event->particles[gamma].parents.push_back(ePos);
	+ event->particles[hOut].parents.push_back(hPos);
	+ event->particles[hOut].parents.push_back(pomeron);
	+ event->particles[pomeron].parents.push_back(gamma);
	+ event->particles[pomeron].parents.push_back(gamma);
	+ event->particles[X].parents.push_back(gamma);
	+ event->particles[quark].parents.push_back(gamma);
	+ event->particles[antiquark].parents.push_back(gamma);
	+
	+ // daughters (again ignore dipole)
	+ event->particles[ePos].daughters.push_back(eOut);
	+ event->particles[ePos].daughters.push_back(gamma);
	+ event->particles[gamma].daughters.push_back(X);
	+ event->particles[hPos].daughters.push_back(pomeron);
	+ event->particles[gamma].daughters.push_back(quark);
	+ event->particles[gamma].daughters.push_back(antiquark);
	+ event->particles[pomeron].daughters.push_back(hOut);
	+ event->particles[hPos].daughters.push_back(hOut);
	+
	+ //fill event structure
	+ double y=mHadronBeamtheVirtualPhoton/(mHadronBeammElectronBeam);
	+ W2=(theVirtualPhoton+mHadronBeam).M2();
	+ mQ2=-theVirtualPhoton.M2();
	+ MX2=(theQuark+theAntiQuark).M2();
	+ double Delta=thePomeron.Pt();
	+
	+ event->t=-Delta*Delta;
	+ event->Q2=mQ2;
	+ event->W=sqrt(W2);
	+ event->y=y;
	+ event->MX=sqrt(MX2);
	+
	+ return true;
	+}
	+
	Index: trunk/src/SartreInclusiveDiffraction.h
	===================================================================
	--- trunk/src/SartreInclusiveDiffraction.h (revision 0)
	+++ trunk/src/SartreInclusiveDiffraction.h (revision 490)
	@@ -0,0 +1,110 @@
	+//==============================================================================
	+// SartreInclusiveDiffraction.h
	+//
	+// Copyright (C) 2024 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 <http://www.gnu.org/licenses/>.
	+//
	+// Author: Tobias Toll
	+// Last update:
	+// $Date: $
	+// $Author: ullrich $
	+//==============================================================================
	+#ifndef SartreInclusiveDiffraction_h
	+#define SartreInclusiveDiffraction_h
	+#include "Event.h"
	+#include "EventGeneratorSettings.h"
	+#include "InclusiveFinalStateGenerator.h"
	+#include "InclusiveDiffractiveCrossSections.h"
	+#include "FrangibleNucleus.h"
	+#include "Enumerations.h"
	+#include "TLorentzVector.h"
	+#include "Math/Functor.h"
	+#include "TUnuran.h"
	+#include <ctime>
	+#include <iostream>
	+#include <vector>
	+
	+using namespace std;
	+
	+class TUnuranMultiContDist;
	+
	+class SartreInclusiveDiffraction {
	+public:
	+ SartreInclusiveDiffraction();
	+ virtual ~SartreInclusiveDiffraction();
	+
	+ 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[4], double upper[4]); // beta, Q2, W, z
	+ virtual double integrationVolume();
	+ EventGeneratorSettings* runSettings();
	+
	+ const FrangibleNucleus* nucleus() const;
	+
	+ void listStatus(ostream& os=cout) const;
	+ time_t runTime() const;
	+
	+ vector<pair<double,double> > kinematicLimits(); // t, Q2, W
	+
	+private:
	+ virtual double calculateTotalCrossSection(double lower[4], double upper[4]); //beta, Q2, W, z
	+
	+private:
	+ SartreInclusiveDiffraction(const SartreInclusiveDiffraction&);
	+ SartreInclusiveDiffraction operator=(const SartreInclusiveDiffraction&);
	+
	+ 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;
	+ InclusiveDiffractiveCrossSections *mCrossSection;
	+ EventGeneratorSettings *mSettings;
	+
	+ double mLowerLimit[4]; // beta, Q2, W2, z
	+ double mUpperLimit[4]; // beta, Q2, W2, z
	+
	+ ROOT::Math::Functor *mPDF_Functor;
	+ TUnuran *mUnuran;
	+ TUnuranMultiContDist *mPDF;
	+
	+ unsigned long mEventCounter;
	+ unsigned long mTriesCounter;
	+
	+ InclusiveFinalStateGenerator mFinalStateGenerator;
	+ double cross_section_wrapper(double, double);
	+
	+};
	+
	+ostream& operator<<(ostream& os, const TLorentzVector&);
	+
	+#endif

	Property changes on: trunk/src/SartreInclusiveDiffraction.h
	___________________________________________________________________
	Added: svn:executable
	## -0,0 +1 ##
	+*
	\ No newline at end of property
	Index: trunk/src/SartreInclusiveDiffraction.cpp
	===================================================================
	--- trunk/src/SartreInclusiveDiffraction.cpp (revision 0)
	+++ trunk/src/SartreInclusiveDiffraction.cpp (revision 490)
	@@ -0,0 +1,761 @@
	+//==============================================================================
	+// SartreInclusiveDiffraction.cpp
	+//
	+// Copyright (C) 2024 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 <http://www.gnu.org/licenses/>.
	+//
	+// Author: Tobias Toll
	+// Last update:
	+// $Date: $
	+// $Author: ullrich $
	+//==============================================================================
	+//
	+// Note:
	+// When not using runcards, the user must first create an instance of
	+// SartreInclusiveDiffraction and then get the settings via one of:
	+// SartreInclusiveDiffraction::runSettings()
	+// EventGeneratorSettings::instance()
	+// Once init() is called settings cannot be changed any more.
	+//==============================================================================
	+#include "Version.h"
	+#include "Kinematics.h"
	+#include "SartreInclusiveDiffraction.h"
	+#include "ModeFinderFunctor.h"
	+#include "Math/BrentMinimizer1D.h"
	+#include "Math/IntegratorMultiDim.h"
	+#include "Math/GSLMinimizer.h"
	+#include "Math/Functor.h"
	+#include "TUnuranMultiContDist.h"
	+#include <string>
	+#include <limits>
	+#include <cmath>
	+#include <iomanip>
	+#include "TH2D.h"
	+#include "TFile.h"
	+#include "TF3.h"
	+
	+using namespace std;
	+
	+#define PR(x) cout << #x << " = " << (x) << endl;
	+
	+SartreInclusiveDiffraction::SartreInclusiveDiffraction()
	+{
	+ mIsInitialized = false;
	+ mCurrentEvent = 0;
	+ mNucleus = 0;
	+ mUpcNucleus= 0;
	+ mSettings = 0;
	+ mPDF_Functor = 0;
	+ mPDF = 0;
	+ mEventCounter = 0;
	+ mTriesCounter = 0;
	+ mTotalCrossSection = 0;
	+ mCrossSection = 0;
	+ mUnuran = 0;
	+ mEvents = 0;
	+ mTries = 0;
	+ mS = 0;
	+ mA = 0;
	+}
	+
	+SartreInclusiveDiffraction::~SartreInclusiveDiffraction()
	+{
	+ delete mNucleus;
	+ delete mUpcNucleus;
	+ delete mPDF_Functor;
	+ delete mPDF;
	+ delete mCrossSection;
	+ delete mUnuran;
	+ delete mCurrentEvent;
	+
	+}
	+
	+bool SartreInclusiveDiffraction::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 inclusive diffraction \|" << endl;
	+ cout << "\| in ep and eA collisions based on the dipole model. \|" << endl;
	+ cout << "\| \|" << endl;
	+ cout << "\| Copyright (C) 2010-2024 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
	+ //
	+ 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);
	+
	+ string upcNucleusName;
	+ 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();
	+ cout << "Dipole model: " << mSettings->dipoleModelName().c_str() << endl;
	+ cout << "Dipole model parameter set: " << mSettings->dipoleModelParameterSetName().c_str() << endl;
	+ cout << "Process is ";
	+ if (mA > 1)
	+ cout << "e + " << mNucleus->name() << " -> e' + " << mNucleus->name() << "' + X"<< endl;
	+ else
	+ cout << "e + p -> e' + p' + X"<< endl;
	+
	+ //
	+ // Print-out seed for reference
	+ //
	+ cout << "Random generator seed: " << mSettings->seed() << endl;
	+
	+ //
	+ // Load in the tables containing the amplitude moments
	+ //
	+ if (!getenv("SARTRE_DIR")) {
	+ cout << "Error, required environment variable 'SARTRE_DIR' is not defined." << endl;
	+ exit(1);
	+ }
	+
	+ //
	+ // 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 2:
	+ // Kinematic limits might overrule boundaries from step 1
	+ //
	+ // let:
	+ // 0: beta
	+ // 1: Q2
	+ // 2: W2
	+ // 3: z
	+ //
	+ // 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 InclusiveDiffractiveCrossSections;
	+
	+ //
	+ // Here we need to put inQ2 > mu02.
	+ // The calculation goes haywire for Q2 becoming too small
	+ //
	+ double mu02 = mCrossSection->dipoleModel()->getParameters()->mu02();
	+ double mf = mCrossSection->dipoleModel()->getParameters()->quarkMass(0);
	+
	+ //Now for the diffractive mass:
	+ double kineMX2min = 4mfmf;//max(mu02, 4mfmf);
	+ // This gives the limits for W2 and for inelasticity y:
	+ double kineW2min = kineMX2min+protonMass2;
	+ double kineYmin = Kinematics::ymin(mS, kineMX2min);
	+ // ...and y gives the limit for Q2, unless it's smaller than the cut-off mu02:
	+ double kineQ2min = Kinematics::Q2min(kineYmin);
	+ kineQ2min = max(kineQ2min, mu02);
	+ //The maximum momenta are given by s:
	+ double kineQ2max = Kinematics::Q2max(mS);
	+ double kineW2max = mS-kineQ2min-protonMass2;
	+
	+ //
	+ // Now for all the fractions, making sure they make sense
	+ //
	+ double kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	+ double kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	+ double kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	+ kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	+ double kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	+ double kineZmax=1-kineZmin;
	+ mLowerLimit[0] = kineBetamin;
	+ mUpperLimit[0] = kineBetamax;
	+ mLowerLimit[1] = kineQ2min;
	+ mUpperLimit[1] = kineQ2max;
	+ mLowerLimit[2] = kineW2min;
	+ mUpperLimit[2] = kineW2max;
	+ mLowerLimit[3] = kineZmin;
	+ mUpperLimit[3] = kineZmax;
	+
	+ // Step 3:
	+ // Deal with user provided limits.
	+ // User settings are ignored (switched off) if min >= max.
	+ // Only allow user limits for Q2 and W2,
	+ // not beta (for now) and z.
	+ //
	+ bool W2orQ2changed=false;
	+ if (mSettings->Wmin() < mSettings->Wmax()) { // 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();
	+ W2orQ2changed=true;
	+ }
	+
	+ 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();
	+ W2orQ2changed=true;
	+ }
	+ }
	+ 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();
	+ W2orQ2changed=true;
	+ }
	+
	+ 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();
	+ W2orQ2changed=true;
	+ }
	+ }
	+ if (W2orQ2changed){//mLowerLimit beta, Q2, W2, z
	+ //kineBetamax = kineQ2max/(kineQ2max+kineMX2min);
	+ mUpperLimit[0] = mUpperLimit[1]/(mUpperLimit[1]+kineMX2min);
	+
	+ //kineBetamin = Kinematics::betamin(kineQ2min, mS, sqrt(kineW2min), 0);
	+ mLowerLimit[0] = Kinematics::betamin(mLowerLimit[1], mS, sqrt(mLowerLimit[2]), 0);
	+
	+ // kineMX2max = (1.-kineBetamin)/kineBetamin*kineQ2max;
	+ // kineMX2max=min(kineMX2max, kineW2max-protonMass2);
	+ kineMX2max = (1.-mLowerLimit[0])/mLowerLimit[0]*mUpperLimit[1];
	+ kineMX2max=min(kineMX2max, mUpperLimit[2]-protonMass2);
	+
	+ // kineZmin=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	+ // kineZmax=1-kineZmin;
	+ mLowerLimit[3]=(1-sqrt(1-4mfmf/kineMX2max))/2.;
	+ mUpperLimit[3]=1-mLowerLimit[3];
	+ }
	+
	+ if (mSettings->betamin() < mSettings->betamax()) { // beta
	+
	+ if (mSettings->betamin() < mLowerLimit[0]) {
	+ cout << "Warning, requested lower limit of beta (" << mSettings->betamin() << ") "
	+ << "is smaller than limit given by lookup tables and/or kinematic range (" << mLowerLimit[0] << "). ";
	+ cout << "Limit has no effect." << endl;
	+ }
	+ else {
	+ mLowerLimit[0] = mSettings->betamin();
	+ }
	+
	+ if (mSettings->betamax() > mUpperLimit[0]) {
	+ cout << "Warning, requested upper limit of Q2 (" << mSettings->betamax() << ") "
	+ << "exceeds limit given by lookup tables and/or kinematic range (" << mUpperLimit[0] << "). ";
	+ cout << "Limit has no effect." << endl;
	+ }
	+ else {
	+ mUpperLimit[0] = mSettings->betamax();
	+ }
	+ }
	+
	+ //
	+ // Check if any phase space is left
	+ //
	+ if (mLowerLimit[0] >= mUpperLimit[0]) {
	+ cout << "Invalid range in beta: beta=[" << mLowerLimit[0] << ", " << mUpperLimit[0] << "]." << endl;
	+ exit(1);
	+ }
	+ if (mLowerLimit[1] >= mUpperLimit[1]) {
	+ cout << "Invalid range in Q2: Q2=[";
	+ cout << mLowerLimit[1] << ", " << mUpperLimit[1] << "]." << endl;
	+ exit(1);
	+ }
	+ if (mLowerLimit[2] >= mUpperLimit[2]) {
	+ cout << "Invalid range in W: W=[" << sqrt(mLowerLimit[2]) << ", " << sqrt(mUpperLimit[2]) << "]." << endl;
	+ exit(1);
	+ }
	+
	+ //
	+ // Print-out limits (all verbose levels)
	+ //
	+ if (true){//}(mSettings->verbose()) {
	+ cout << "Sartre was thrown into the world, restricted by kinematics and user inputs alike:" << endl;
	+ cout << setw(10) << " beta=[" << 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;
	+ cout << setw(10) << " z=[" << mLowerLimit[3] << ", " << mUpperLimit[3] << "]" << endl;
	+ cout << setw(10) << " MX2=[" << kineMX2min << ", " << kineMX2max << "]" << endl;
	+ cout << setw(10) << " y=[" << kineYmin << ", 1" << "]" << endl;
	+
	+ }
	+ //
	+ // 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 beta is not obvious. The mode should be at z=0.5.
	+ // The approach here is to use the BrentMinimizer1D that
	+ // performs first a scan a then a Brent fit.
	+ //
	+ mCrossSection->setCheckKinematics(false);
	+ double theMode[4];
	+
	+ //
	+ // Find the mode.
	+ // Assume that the mode is at W2=W2max, Q2=Q2min, z=0.5
	+ // Start with beta=0.5 => MX2=Q2
	+ // It is kinemtaically easier to use MX2 instead of beta here.
	+ //
	+ if (mSettings->verbose()) cout << "Finding mode of pdf:" << endl;
	+ theMode[1] = mLowerLimit[1]; // Q2
	+ theMode[2] = mUpperLimit[2]; // W2
	+ theMode[3] = 0.5; // z
	+
	+ double MX2min=mLowerLimit[1]*(1-mUpperLimit[0])/mUpperLimit[0];
	+ MX2min=max(kineMX2min, MX2min);
	+ double MX2max=mLowerLimit[1]*(1-mLowerLimit[0])/mLowerLimit[0];
	+ MX2max=min(kineMX2max, MX2max);
	+ InclusiveDiffractionModeFinderFunctor modeFunctor(mCrossSection, theMode[1], theMode[2], theMode[3], MX2min, MX2max);
	+
	+ ROOT::Math::BrentMinimizer1D minimizer;
	+ minimizer.SetFunction(modeFunctor, MX2min, MX2max);
	+ minimizer.SetNpx(static_cast<int>(mUpperLimit[0]-mLowerLimit[0]));
	+ ok = minimizer.Minimize(100000, 0, 1.e-8);
	+ if (! ok) {
	+ cout << "Error, failed to find mode of pdf." << endl;
	+ exit(1);
	+ }
	+ theMode[0] = minimizer.XMinimum(); // Mx2
	+ double MX2=theMode[0];
	+ theMode[0]=theMode[1]/(theMode[1]+MX2); //Mx2->beta
	+ double crossSectionAtMode = (*mCrossSection)(MX2, theMode[1], theMode[2], theMode[3]);
	+ if (mSettings->verbose()) {
	+ cout << "\tlocation: beta=" << theMode[0] << ", Q2=" << theMode[1] << ", W=" << sqrt(theMode[2]) << ", z=" << theMode[3];
	+ cout << "; value: " << crossSectionAtMode << endl;
	+ }
	+ mCrossSection->setCheckKinematics(true);
	+
	+ // domain and mode for Q2 -> log(Q2)
	+ 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;
	+
	+ mPDF_Functor = new ROOT::Math::Functor(mCrossSection, &InclusiveDiffractiveCrossSections::unuranPDF, 4);
	+
	+ 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[4];
	+ for (int i=0; i<1; i++) {
	+ mUnuran->SampleMulti(xrandom);
	+ }
	+
	+ mEventCounter = 0;
	+ mTriesCounter = 0;
	+ mIsInitialized = true;
	+ cout << "Sartre for InclusiveDiffraction is initialized." << endl << endl;
	+
	+ return true;
	+}
	+
	+double SartreInclusiveDiffraction::cross_section_wrapper(double* xx, double* par){
	+ double arg[4]={xx[0], xx[1], xx[2], par[0]};
	+ double result=mCrossSection->unuranPDF(arg);
	+ return -result;
	+}
	+
	+bool SartreInclusiveDiffraction::init(const string& str) // overloaded version of init()
	+{
	+ if (str.empty())
	+ return init();
	+ else
	+ return init(str.c_str());
	+}
	+
	+vector<pair<double,double> > SartreInclusiveDiffraction::kinematicLimits()
	+{
	+ vector<pair<double,double> > 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* SartreInclusiveDiffraction::generateEvent()
	+{
	+ if (!mIsInitialized) {
	+ cout << "SartreInclusiveDiffraction::generateEvent(): Error, Sartre is not initialized yet." << endl;
	+ cout << " Call init() before trying to generate events." << endl;
	+ return 0;
	+ }
	+
	+ double xrandom[4];
	+
	+ //
	+ // Generate one event
	+ //
	+ while (true) {
	+ mTriesCounter++;
	+ delete mCurrentEvent;
	+ mCurrentEvent = new Event;
	+ //
	+ // xrandom[i]
	+ // i=0: beta
	+ // i=1: Q2
	+ // i=2: W2
	+ // i=3: z
	+ //
	+ mUnuran->SampleMulti(xrandom);
	+ xrandom[1] = exp(xrandom[1]); // log(Q2) -> Q2
	+ double MX2=xrandom[1]*(1-xrandom[0])/xrandom[0];
	+ bool isValidEvent;
	+ isValidEvent = Kinematics::valid( mS, xrandom[0], xrandom[1], xrandom[2], xrandom[3], sqrt(MX2), true, (mSettings->verboseLevel() > 1));
	+ if (!isValidEvent) {
	+ if (mSettings->verboseLevel() > 2)
	+ cout << "SartreInclusiveDiffraction::generateEvent(): event rejected, not within kinematic limits" << endl;
	+ continue;
	+ }
	+ //
	+ // Fill beam particles in Event structure
	+ // Kinematics for eA is reported as 'per nucleon'
	+ //
	+ mCurrentEvent->eventNumber = mEventCounter;
	+ mCurrentEvent->beta = xrandom[0]; // beta
	+ mCurrentEvent->Q2 = xrandom[1]; // Q2
	+ mCurrentEvent->x = Kinematics::x(xrandom[1], xrandom[2]); // x
	+ mCurrentEvent->xpom=mCurrentEvent->x/xrandom[0]; //xpom=x/beta
	+ mCurrentEvent->y = Kinematics::y(xrandom[1], mCurrentEvent->x, mS); // y
	+ mCurrentEvent->s = mS; // s
	+ mCurrentEvent->W = sqrt(xrandom[2]);
	+ mCurrentEvent->z = xrandom[3];
	+ mCurrentEvent->MX = sqrt(MX2);
	+ mCurrentEvent->polarization = mCrossSection->polarizationOfLastCall();
	+ // mCurrentEvent->diffractiveMode = mCrossSection->diffractiveModeOfLastCall();
	+ mCurrentEvent->quarkSpecies = mCrossSection->quarkSpeciesOfLastCall();
	+ mCurrentEvent->quarkMass=mCrossSection->quarkMassOfLastCall();
	+ mCurrentEvent->crossSection=mCrossSection->crossSectionOfLastCall();
	+ double dsigmadb2=0;
	+ double z=xrandom[3];
	+ double eps2=z(1-z)(mCurrentEvent->Q2+MX2);
	+ double r_average=1./sqrt(eps2);
	+ if(mA==1)
	+ dsigmadb2=mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom);
	+ else{
	+ dsigmadb2=mCrossSection->dipoleModel()->coherentDsigmadb2(r_average*hbarc, 0, mCurrentEvent->xpom);
	+ if(mCurrentEvent->xpom<1e-2){
	+ cout<<"==============="<<endl;
	+ PR(dsigmadb2);
	+ PR(mCrossSection->dipoleModel()->dsigmadb2ep(r_average*hbarc, 0, mCurrentEvent->xpom));
	+ PR(mCurrentEvent->xpom);
	+ PR(r_average*hbarc);
	+ cout<<"==============="<<endl;
	+ }
	+ }
	+ mCurrentEvent->Omega=-2.*log(1-dsigmadb2/2.);
	+ 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(mA, mCurrentEvent);
	+ }
	+ else {
	+ ok = mFinalStateGenerator.generate(mA, mCurrentEvent);
	+ }
	+ if (!ok) {
	+ if (mSettings->verboseLevel() > 1) cout << "SartreInclusiveDiffraction::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<BreakupProduct>& products = mNucleus->breakupProducts();
	+ for (int i=0; i<nFragments; i++) {
	+ Particle fragment;
	+ fragment.index = mCurrentEvent->particles.size();
	+ fragment.pdgId = products[i].pdgId;
	+ fragment.status = 1;
	+ fragment.p = products[i].p*(1/static_cast<double>(products[i].A));
	+ fragment.parents.push_back(indexOfScatteredHadron);
	+ mCurrentEvent->particles.push_back(fragment);
	+ }
	+ }
	+ return mCurrentEvent;
	+}
	+
	+double SartreInclusiveDiffraction::totalCrossSection()
	+{
	+ if (mTotalCrossSection == 0) {
	+ //
	+ // Limits of integration in beta, Q2, W2, z
	+ //
	+ double xmin[4];
	+ double xmax[4];
	+ copy(mLowerLimit, mLowerLimit+4, xmin);
	+ copy(mUpperLimit, mUpperLimit+4, 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 SartreInclusiveDiffraction::totalCrossSection(double lower[4], double upper[4]) // beta, Q2, W, z
	+{
	+ lower[2] = lower[2]; upper[2] = upper[2]; // W -> W2
	+ double result = calculateTotalCrossSection(lower, upper);
	+ return result;
	+}
	+
	+double SartreInclusiveDiffraction::integrationVolume(){
	+ //
	+ // Limits of integration in beta, Q2, W2, z
	+ //
	+ double xmin[4];
	+ double xmax[4];
	+ copy(mLowerLimit, mLowerLimit+4, xmin);
	+ copy(mUpperLimit, mUpperLimit+4, xmax);
	+ //
	+ // At this point mLowerLimit[1] and mUpperLimit[1]
	+ // are in log(Q2)
	+ //
	+ xmin[1] = exp(xmin[1]); // log Q2 limit
	+ xmax[1] = exp(xmax[1]); // log Q2 limit
	+
	+ double volume = (xmax[0]-xmin[0]) * (xmax[1]-xmin[1]) * (xmax[2]-xmin[2]) * (xmax[3]-xmin[3]);
	+ return volume; //GeV4
	+}
	+
	+EventGeneratorSettings* SartreInclusiveDiffraction::runSettings()
	+{
	+ return EventGeneratorSettings::instance();
	+}
	+
	+double SartreInclusiveDiffraction::calculateTotalCrossSection(double lower[4], double upper[4])
	+{
	+ double result = 0;
	+
	+ if (!mIsInitialized) {
	+ cout << "SartreInclusiveDiffraction::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 = 1e-1;//5e-4;
	+ unsigned int numberofcalls = 10000;//1000000;
	+ ROOT::Math::Functor wfL((*mCrossSection), 4);
	+ ROOT::Math::IntegratorMultiDim ig(ROOT::Math::IntegrationMultiDim::kADAPTIVE,
	+ 0, precision, numberofcalls);
	+
	+ 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<float>::epsilon()) {
	+ cout << "SartreInclusiveDiffraction::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);
	+ }
	+
	+ return result;
	+}
	+
	+const FrangibleNucleus* SartreInclusiveDiffraction::nucleus() const {return mNucleus;}
	+
	+void SartreInclusiveDiffraction::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 SartreInclusiveDiffraction::runTime() const
	+{
	+ time_t now = time(0);
	+ return now-mStartTime;
	+}
	+
	+//==============================================================================
	+//
	+// 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;
	+}

	Property changes on: trunk/src/SartreInclusiveDiffraction.cpp
	___________________________________________________________________
	Added: svn:executable
	## -0,0 +1 ##
	+*
	\ No newline at end of property

File Metadata

Mime Type: text/x-diff
Expires: Tue, Nov 19, 6:01 PM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 3805478
Default Alt Text: (50 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions