Page MenuHomeHEPForge
Files F8309079

No OneTemporary
Actions

View Options

Index: trunk/examples/sartreInclusiveDiffractionRuncard.txt
===================================================================
--- trunk/examples/sartreInclusiveDiffractionRuncard.txt (revision 0)
+++ trunk/examples/sartreInclusiveDiffractionRuncard.txt (revision 511)
@@ -0,0 +1,108 @@
+//==============================================================================
+// sartreRuncard.txt
+//
+// 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 <http://www.gnu.org/licenses/>.
+//
+// Author: Thomas Ullrich
+// Last update:
+// $Date: 2019-06-10 23:53:14 +0530 (Mon, 10 Jun 2019) $
+// $Author: ullrich $
+//==============================================================================
+//
+// Example Runcard for Sartre Event Generator.
+//
+// Comments start with a # or //
+//
+// Name and value are separated by a "=": name = value
+// (alternatively ":" can be used as well)
+//==============================================================================
+#
+# Define beams
+#
+eBeamEnergy = 30
+hBeamEnergy = 920
+A = 1
+
+#
+# Number of events and printout frequency
+#
+numberOfEvents = 10
+timesToShow = 10
+
+#
+# Set verbosity
+#
+verbose = false
+verboseLevel = 2
+#
+# Rootfile
+#
+rootfile = example.root
+
+#
+# Dipole model
+#
+# dipoleModel: can be bSat or bNonSat
+# dipoleModelParameterSet:
+# KMW: older set but most complete
+# HMPZ: new parameters with still limited set
+# of tables. For eA use KMW.
+# STU: new parameters with still limited set
+# (uses large dipole damping)
+#
+dipoleModel = bSat
+dipoleModelParameterSet = STU
+
+#
+# User variable
+# If you use the provided sartreMain.cpp this is
+# used for vector meson decays
+# PDG: pi+ = 211, K+ = 321, e- = 11, mu- = 13
+# with BR=100%
+#
+userInt = 11
+
+#
+# Kinematic range min > max means no limits (given by table range)
+#
+Q2min = 10
+Q2max = -11
+Wmin = 10.
+Wmax = -10.049876
+betamin = 0.4
+betamax = 1
+
+#
+# Misc
+#
+enableNuclearBreakup = false
+maxNuclearExcitationEnergy = 3.0
+
+#
+# Random generator seed (if not given current time is used)
+#
+#seed = 2011987
+
+#
+# User parameters
+#
+# userDouble = 0.
+# userString = "Hello World!"
+# userInt = 0
+
+#
+# Expert flags
+#
+# applyPhotonFlux = true
Property changes on: trunk/examples/sartreInclusiveDiffractionRuncard.txt
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Index: trunk/examples/cmake/modules/FindPythia8.cmake
===================================================================
--- trunk/examples/cmake/modules/FindPythia8.cmake (revision 0)
+++ trunk/examples/cmake/modules/FindPythia8.cmake (revision 511)
@@ -0,0 +1,94 @@
+ ################################################################################
+ # Copyright (C) 2014 GSI Helmholtzzentrum fuer Schwerionenforschung GmbH #
+ # #
+ # This software is distributed under the terms of the #
+ # GNU Lesser General Public Licence (LGPL) version 3, #
+ # copied verbatim in the file "LICENSE" #
+ ################################################################################
+# - Try to find PYTHIA8 instalation
+# Once done this will define
+# PYTHIA8_INCLUDE_DIR where to locate Pythia.h file
+# PYTHIA8_LIBRARY where to find the libpythia8 library
+# PYTHIA8_<lib>_LIBRARY Addicional libraries
+# PYTHIA8_LIBRARIES (not cached) the libraries to link against to use Pythia8
+# PYTHIA8_FOUND if false, you cannot build anything that requires Pythia8
+# PYTHIA8_VERSION version of Pythia8 if found
+################################################################################
+
+set(_pythia8dirs
+ ${PYTHIA8_DIR}
+ $ENV{PYTHIA8_DIR}
+ ${SIMPATH}
+ ${SIMPATH}/generators
+ /usr
+ /opt/pythia8)
+
+find_path(PYTHIA8_INCLUDE_DIR
+ NAMES Pythia8/Pythia.h pythia8/Pythia.h pythia8/pythia.h Pythia8/pythia.h
+ HINTS ${_pythia8dirs}
+ PATH_SUFFIXES include include/Pythia8 include/pythia8
+ DOC "Specify the directory containing Pythia.h.")
+
+if(PYTHIA8_INCLUDE_DIR MATCHES "include$")
+ get_filename_component(Pythia8_PREFIX ${PYTHIA8_INCLUDE_DIR}/.. ABSOLUTE)
+elseif()
+ get_filename_component(Pythia8_PREFIX ${PYTHIA8_INCLUDE_DIR}/../.. ABSOLUTE)
+endif()
+
+find_library(PYTHIA8_LIBRARY
+ NAMES pythia8 Pythia8
+ HINTS ${_pythia8dirs}
+ PATH_SUFFIXES lib
+ DOC "Specify the Pythia8 library here.")
+
+find_path(PYTHIA8_LIB_DIR NAMES libpythia8.so libpythia8.dylib PATHS
+ $ENV{PYTHIA8_DIR}
+ ${PYTHIA8_DIR}
+ ${SIMPATH}
+ ${SIMPATH}/generators
+ PATH_SUFFIXES lib
+ )
+
+find_library(PYTHIA8_hepmcinterface_LIBRARY
+ NAMES hepmcinterface pythia8tohepmc
+ HINTS ${_pythia8dirs}
+ PATH_SUFFIXES lib)
+
+find_library(PYTHIA8_lhapdfdummy_LIBRARY
+ NAMES lhapdfdummy
+ HINTS ${_pythia8dirs}
+ PATH_SUFFIXES lib)
+
+
+find_path(PYTHIA8DATA
+ NAMES MainProgramSettings.xml
+ HINTS
+ ${Pythia8_PREFIX}/share/Pythia8/xmldoc
+ ${Pythia8_PREFIX}/share/pythia8/xmldoc)
+
+foreach(_lib PYTHIA8_LIBRARY PYTHIA8_hepmcinterface_LIBRARY PYTHIA8_lhapdfdummy_LIBRARY)
+ if(${_lib})
+ set(PYTHIA8_LIBRARIES ${PYTHIA8_LIBRARIES} ${${_lib}})
+ endif()
+endforeach()
+set(PYTHIA8_INCLUDE_DIRS ${PYTHIA8_INCLUDE_DIR} ${PYTHIA8_INCLUDE_DIR}/Pythia8 )
+
+find_file(Pythia8_VERSION_FILE
+ NAMES Pythia.h pythia.h
+ HINTS ${PYTHIA8_INCLUDE_DIR} ${PYTHIA8_INCLUDE_DIR}/Pythia8 ${PYTHIA8_INCLUDE_DIR}/pythia8)
+
+if(Pythia8_VERSION_FILE)
+ file(READ "${Pythia8_VERSION_FILE}" __version_raw)
+ string(REGEX MATCH "#define PYTHIA_VERSION ([0-9\.]+)"
+ __version_raw "${__version_raw}"
+ )
+ set(Pythia8_VERSION "${CMAKE_MATCH_1}")
+ unset(__version_raw)
+endif()
+
+# handle the QUIETLY and REQUIRED arguments and set PYTHIA8_FOUND to TRUE if
+# all listed variables are TRUE
+
+include(FindPackageHandleStandardArgs)
+find_package_handle_standard_args(Pythia8 DEFAULT_MSG PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY PYTHIA8DATA)
+mark_as_advanced(PYTHIA8_INCLUDE_DIR PYTHIA8_LIBRARY PYTHIA8_hepmcinterface_LIBRARY PYTHIA8_lhapdfdummy_LIBRARY)
Index: trunk/examples/sartreInclusiveDiffractionMain.cpp
===================================================================
--- trunk/examples/sartreInclusiveDiffractionMain.cpp (revision 0)
+++ trunk/examples/sartreInclusiveDiffractionMain.cpp (revision 511)
@@ -0,0 +1,512 @@
+//==============================================================================
+// sartreMain.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 <http://www.gnu.org/licenses/>.
+//
+// Author: Thomas Ullrich
+// Last update:
+// $Date: 2019-03-09 00:43:19 +0530 (Sat, 09 Mar 2019) $
+// $Author: ullrich $
+//==============================================================================
+//
+// Example main program. Use to get started.
+//
+//==============================================================================
+#include <iostream>
+#include <cmath>
+#include "TTree.h"
+#include "TFile.h"
+#include "SartreInclusiveDiffraction.h"
+#include "DipoleModelParameters.h"
+#include "TGenPhaseSpace.h"
+#include "Settings.h"
+#include "TH1D.h"
+#include "TROOT.h"
+#include "TLorentzVector.h"
+#include "TClonesArray.h"
+#include "InclusiveDiffractiveCrossSections.h"
+#include "DglapEvolution.h"
+#include "Kinematics.h"
+
+#define PR(x) cout << #x << " = " << (x) << endl;
+
+using namespace std;
+InclusiveDiffractiveCrossSections* myICS=0;
+
+double uiCrossSection(double* var, double* par){
+ double beta=par[0];
+ double Q2=par[1];
+ double W2=par[2];
+ double z=var[0];
+
+ return myICS->dsigdbetadQ2dW2dz_total_checked(beta, Q2, W2, z);
+}
+
+double xp_sigmar(double Q2, double beta, double xpom, double s){
+ double x=beta*xpom;
+ double W2=Kinematics::W2(Q2, x);
+ double y=Kinematics::y(Q2, x, s);
+ TF1 fsigma("fsigma", uiCrossSection, 0, 1, 3);
+ fsigma.SetParameter(0, beta);
+ fsigma.SetParameter(1, Q2);
+ fsigma.SetParameter(2, W2);
+ ROOT::Math::WrappedTF1 wfsigma(fsigma);
+ ROOT::Math::GaussIntegrator giSigma;
+ giSigma.SetFunction(wfsigma);
+ giSigma.SetAbsTolerance(0.);
+ giSigma.SetRelTolerance(1e-5);
+
+ double cs=giSigma.Integral(0, 1); //nb/GeV4
+ double prefactor= 4*M_PI*alpha_em*alpha_em* xpom*xpom /(Q2*Q2*Q2);
+ prefactor*=1-y+y*y/2; //GeV-6
+ double sigma_r=cs/prefactor; //nb*GeV2
+ sigma_r /= 1e7; //fm2*GeV2
+ sigma_r /= hbarc2; //no unit;
+ return xpom*sigma_r;
+}
+
+void reducedCrossSections(double s){
+
+ myICS = new InclusiveDiffractiveCrossSections();
+
+ double Q2val[80]={2.5, 2.5, 5.1, 8.8, 2.5, 2.5, 2.5, 5.1, 5.1, 8.8, 8.8, 15.3, 26.5, 46, 2.5, 2.5, 2.5, 2.5, 5.1, 5.1, 5.1, 8.8, 8.8, 8.8, 15.3, 15.3, 26.5, 26.5, 46, 46, 2.5, 2.5, 2.5, 2.5, 5.1, 5.1, 5.1, 5.1, 8.8, 8.8, 8.8, 8.8, 15.3, 15.3, 15.3, 26.5, 26.5, 26.5, 46, 46, 46, 80, 80, 200, 2.5, 2.5, 2.5, 2.5, 5.1, 5.1, 5.1, 5.1, 8.8, 8.8, 8.8, 8.8, 15.3, 15.3, 15.3, 15.3, 26.5, 26.5, 46, 46, 46, 46, 80, 80, 200, 200};
+ double betaval[80]={0.1780, 0.5620, 0.5620, 0.5620, 0.0562, 0.1780, 0.5620, 0.1780, 0.5620, 0.1780, 0.5620, 0.5620, 0.5620, 0.8160, 0.0178, 0.0562, 0.1780, 0.5620, 0.0562, 0.1780, 0.5620, 0.0562, 0.1780, 0.5620, 0.1780, 0.5620, 0.1780, 0.5620, 0.5620, 0.8160, 0.0056, 0.0178, 0.0562, 0.1780, 0.0178, 0.0562, 0.1780, 0.5620, 0.0178, 0.0562, 0.1780, 0.5620, 0.0562, 0.1780, 0.5620, 0.0562, 0.1780, 0.5620, 0.1780, 0.5620, 0.8160, 0.1780, 0.5620, 0.5620, 0.0056, 0.0178, 0.0562, 0.1780, 0.0056, 0.0178, 0.0562, 0.1780, 0.0178, 0.0562, 0.1780, 0.5620, 0.0178, 0.0562, 0.1780, 0.5620, 0.0562, 0.1780, 0.0562, 0.1780, 0.5620, 0.8160, 0.1780, 0.5620, 0.1780, 0.5620};
+ double xpomval[80]={0.0003, 0.0003, 0.0003, 0.0003, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0009, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0025, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0085, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160, 0.0160};
+ double xsigma_rval[80]={0.0156, 0.0214, 0.0275, 0.0250, 0.0114, 0.0102, 0.0172, 0.0121, 0.0187, 0.0177, 0.0207, 0.0180, 0.0241, 0.0145, 0.0099, 0.0074, 0.0068, 0.0110, 0.0107, 0.0118, 0.0153, 0.0125, 0.0129, 0.0166, 0.0136, 0.0173, 0.0138, 0.0189, 0.0196, 0.0131, 0.0101, 0.0076, 0.0064, 0.0074, 0.0120, 0.0088, 0.0095, 0.0137, 0.0128, 0.0106, 0.0104, 0.0142, 0.0134, 0.0111, 0.0162, 0.0157, 0.0126, 0.0140, 0.0121, 0.0135, 0.0110, 0.0206, 0.0116, 0.0109, 0.0093, 0.0073, 0.0068, 0.0116, 0.0135, 0.0111, 0.0088, 0.0075, 0.0124, 0.0108, 0.0090, 0.0102, 0.0176, 0.0122, 0.0098, 0.0151, 0.0150, 0.0113, 0.0163, 0.0133, 0.0124, 0.0092, 0.0133, 0.0090, 0.0145, 0.0093};
+ double deltaval[80]={12, 19, 16, 13, 17, 8.8, 31, 18, 12, 9.7, 10, 11, 25, 22, 18, 12, 12, 36, 6.8, 7.4, 8.7, 7.1, 6.1, 10, 6.2, 7.2, 7.6, 12, 8.6, 20, 25, 12, 13, 15, 10, 6.7, 6.7, 27, 7.8, 6.3, 5.9, 13, 9.7, 6.2, 6.8, 12, 7.5, 7.6, 8.9, 8.5, 19, 16, 13, 21, 15, 13, 16, 20, 7.4, 6.5, 6.6, 19, 6.0, 6.4, 7.1, 22, 7.4, 6.8, 6.8, 22, 7.2, 7.9, 11, 8.1, 8.9, 28, 14, 16, 21, 24};
+
+ double beta=0.00562;//betaval[i];
+ double xpom=0.0085;//xpomval[i];
+ double Q2=2.5;
+ cout<<"sartre_0_0=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0085_beta_0_0178=[[2.5,5.1,8.8,],[0.0076,0.012,0.0128],
+
+ beta=.0178;
+ xpom=0.0085;
+ double Q2_1[3]={2.5,5.1,8.8};
+ cout<<"sartre_0_1=[";
+ for(int i=0; i<3; i++){
+ cout<<xp_sigmar(Q2_1[i], beta, xpom, s);
+ if(i<2) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0085_beta_0_0562=[[2.5,5.1,8.8,15.3,26.5],
+ beta=.0562;
+ xpom=0.0085;
+ double Q2_2[5]={2.5,5.1,8.8,15.3,26.5};
+ cout<<"sartre_0_2=[";
+ for(int i=0; i<5; i++){
+ cout<<xp_sigmar(Q2_2[i], beta, xpom, s);
+ if(i<4) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0085_beta_0_178=[[2.5,5.1,8.8,15.3,26.5,46,80]
+ xpom=0.0085;
+ beta=0.178;
+ double Q2_3[7]={2.5,5.1,8.8,15.3,26.5,46,80};
+ cout<<"sartre_0_3=[";
+ for(int i=0; i<7; i++){
+ cout<<xp_sigmar(Q2_3[i], beta, xpom, s);
+ if(i<6) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0085_beta_0_562=[[5.1,8.8,15.3,26.5,46.0,80.0,200.0]
+ xpom=0.0085;
+ beta=.562;
+ double Q2_4[7]={5.1,8.8,15.3,26.5,46.0,80.0,200.0};
+ cout<<"sartre_0_4=[";
+ for(int i=0; i<7; i++){
+ cout<<xp_sigmar(Q2_4[i], beta, xpom, s);
+ if(i<6) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0085_beta_0_816=[[46.0]
+ xpom=0.0085;
+ beta=0.816;
+ Q2=46.;
+ cout<<"sartre_0_5=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0025_beta_0_0178=[[2.5],
+ xpom=0.0025;
+ beta=0.0178;
+ Q2=2.5;
+ cout<<"sartre_1_1=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0025_beta_0_0562=[[2.5,5.1,8.8]
+ xpom=0.0025;
+ beta=0.0562;
+ double Q2_5[3]={2.5,5.1,8.8};
+ cout<<"sartre_1_2=[";
+ for(int i=0; i<3; i++){
+ cout<<xp_sigmar(Q2_5[i], beta, xpom, s);
+ if(i<2) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0025_beta_0_178=[[2.5,5.1,8.8,15.3,26.5]
+ xpom=0.0025;
+ beta=0.178;
+ double Q2_6[5]={2.5,5.1,8.8,15.3,26.5};
+ cout<<"sartre_1_3=[";
+ for(int i=0; i<5; i++){
+ cout<<xp_sigmar(Q2_6[i], beta, xpom, s);
+ if(i<4) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0025_beta_0_562=[[2.5,5.1,8.8,15.3,26.5,46.0]
+ xpom=0.0025;
+ beta=0.562;
+ double Q2_7[6]={2.5,5.1,8.8,15.3,26.5,46.0};
+ cout<<"sartre_1_4=[";
+ for(int i=0; i<6; i++){
+ cout<<xp_sigmar(Q2_7[i], beta, xpom, s);
+ if(i<5) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0025_beta_0_816=[[46],
+ xpom=0.0025;
+ beta=0.816;
+ Q2=46.;
+ cout<<"sartre_1_5=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0009_beta_0_0562=[[2.5]
+ xpom=0.0009;
+ beta=0.0562;
+ Q2=2.5;
+ cout<<"sartre_2_2=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0009_beta_0_178=[[2.5,5.1,8.8]
+ xpom=0.0009;
+ beta=0.178;
+ double Q2_8[3]={2.5,5.1,8.8};
+ cout<<"sartre_2_3=[";
+ for(int i=0; i<3; i++){
+ cout<<xp_sigmar(Q2_8[i], beta, xpom, s);
+ if(i<2) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0009_beta_0_562=[[2.5,5.1,8.8,15.3,26.5]
+ xpom=0.0009;
+ beta=0.562;
+ double Q2_9[5]={2.5,5.1,8.8,15.3,26.5};
+ cout<<"sartre_2_4=[";
+ for(int i=0; i<5; i++){
+ cout<<xp_sigmar(Q2_9[i], beta, xpom, s);
+ if(i<4) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+ // xp_0_0009_beta_0_816=[[46.0]
+ xpom=0.0009;
+ beta=0.816;
+ Q2=46.;
+ cout<<"sartre_2_5=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+ // xp_0_0003_beta_0_178=[[2.5],[0.0156]
+ xpom=0.0003;
+ beta=0.178;
+ Q2=2.5;
+ cout<<"sartre_3_3=["<<xp_sigmar(Q2, beta, xpom, s)<<"]"<<endl;
+
+// xp_0_0003_beta_0_562=[[2.5,5.1,8.8]
+ xpom=0.0003;
+ beta=0.562;
+ double Q2_10[3]={2.5,5.1,8.8};
+ cout<<"sartre_3_4=[";
+ for(int i=0; i<3; i++){
+ cout<<xp_sigmar(Q2_10[i], beta, xpom, s);
+ if(i<2) cout<<", ";
+ }
+ cout<<"]"<<endl;
+
+}
+void checksForCrossSections(){
+ cout<<"#######################################"<<endl;
+ cout<<"### Doing some sanity checks ###"<<endl;
+ cout<<"#######################################"<<endl;
+ TH1D cs_beta_T("cs_beta_T", "cs_beta_T", 100, 0, 1);
+ TH1D cs_beta_L("cs_beta_L", "cs_beta_L", 100, 0, 1);
+ InclusiveDiffractiveCrossSections myCS;
+ double Q2=5;
+ double xpom=1e-3;
+ double z=0.1;
+ cout<<"Filling the cross section histograms binned in beta:"<<endl;
+ for(int i=1; i<=100; i++){
+ cout<<"Filling bin# "<<i<<endl;
+ double beta=cs_beta_T.GetBinCenter(i);
+ double MX2=(1-beta)/beta*Q2;
+ double W2=Kinematics::W2(Q2, xpom, MX2);
+ myCS.setQuarkIndex(0);
+ double csT=myCS.dsigmadbetadz_T(beta, Q2, W2, z);
+ double csL=myCS.dsigmadbetadz_L(beta, Q2, W2, z);
+ PR(csT);
+ PR(csL);
+ cs_beta_T.SetBinContent(i, csT);
+ cs_beta_L.SetBinContent(i, csL);
+ }
+ TFile *hfile = new TFile("CStestsBeta.root","RECREATE");
+ cs_beta_T.Write();
+ cs_beta_L.Write();
+ hfile->Write();
+ delete hfile;
+}
+
+struct rootSartreEvent {
+ double t;
+ double beta;
+ double Q2;
+ double x;
+ double s;
+ double y;
+ double W;
+ double z;
+ double MX;
+ double xpom;
+ int iEvent;
+ int pol; // 0=transverse or 1=longitudinal
+ int dmode; // 0=coherent, 1=Incoherent
+// double crossSection;
+};
+
+rootSartreEvent myRootSartreEvent;
+void randomlyReverseBeams(Event* ); // for UPC mode
+
+int main(int argc, char *argv[])
+{
+ //
+ // Check command line arguments
+ //
+ if (! (argc == 2 || argc == 3) ) {
+ cout << "Usage: " << argv[0] << " runcard [rootfile]" << endl;
+ return 2;
+ }
+ string runcard = argv[1];
+
+
+ //To run some sanity checks on the cross sections:
+ // checksForCrossSections();
+ // exit(1);
+
+ //
+ // Create the generator and initialize it.
+ // Once initialized you cannot (should not) change
+ // the settings w/o re-initialing sartre.
+ //
+ SartreInclusiveDiffraction sartre;
+ //
+ // DGLAP Evolution can be speed up by using lookup tables
+ //
+ DglapEvolution &dglap = DglapEvolution::instance();
+ dglap.generateLookupTable(1000, 1000);
+ dglap.useLookupTable(true);
+
+ bool ok = sartre.init(runcard);
+ if (!ok) {
+ cerr << "Initialization of sartre failed." << endl;
+ return 1;
+ }
+
+ EventGeneratorSettings* settings = sartre.runSettings();
+
+ //
+ // ROOT file
+ // Use the one from command line unless not given
+ // in which case the one in the runcard is used.
+ //
+ string rootfile;
+ if (argc == 3) {
+ rootfile = argv[2];
+ settings->setRootfile(argv[2]);
+ }
+ else
+ rootfile = settings->rootfile();
+
+ //
+ // Print out all Sartre settings
+ //
+ settings->list();
+
+ TFile *hfile = 0;
+ if (rootfile.size()) {
+ hfile = new TFile(rootfile.c_str(),"RECREATE");
+ cout << "ROOT file is '" << rootfile.c_str() << "'." << endl;
+ }
+// double s=4*settings->electronBeamEnergy()*settings->hadronBeamEnergy();
+// reducedCrossSections(s);
+// exit(1);
+//
+ //
+ // Setup ROOT tree
+ //
+ TLorentzVector *eIn = new TLorentzVector;
+ TLorentzVector *pIn = new TLorentzVector;
+ TLorentzVector *X = new TLorentzVector;
+ TLorentzVector *eOut = new TLorentzVector;
+ TLorentzVector *pOut = new TLorentzVector;
+ TLorentzVector *gamma = new TLorentzVector;
+ TLorentzVector *quark = new TLorentzVector;
+ TLorentzVector *antiquark = new TLorentzVector;
+
+ TClonesArray protons("TLorentzVector");
+ TClonesArray neutrons("TLorentzVector");
+ TClonesArray remnants("TLorentzVector");
+
+ TTree tree("tree","sartre");
+ tree.Branch("event", &myRootSartreEvent.t,
+ "t/D:beta/D:Q2/D:x/D:s/D:y/D:W/D:z/D:MX/D:xpom/D:iEvent/I:pol/I:dmode/I");
+ tree.Branch("eIn", "TLorentzVector", &eIn, 32000, 0);
+ tree.Branch("pIn", "TLorentzVector", &pIn, 32000, 0);
+ tree.Branch("X", "TLorentzVector", &X, 32000, 0);
+ tree.Branch("eOut", "TLorentzVector", &eOut, 32000, 0);
+ tree.Branch("pOut", "TLorentzVector", &pOut, 32000, 0);
+ tree.Branch("gamma","TLorentzVector", &gamma, 32000, 0);
+ tree.Branch("quark", "TLorentzVector", &quark, 32000, 0);
+ tree.Branch("antiquark", "TLorentzVector", &antiquark, 32000, 0);
+
+ if(settings->enableNuclearBreakup()){
+ tree.Branch("protons", &protons);
+ tree.Branch("neutrons", &neutrons);
+ tree.Branch("nuclearRemnants", &remnants);
+ }
+
+ int nPrint;
+ if (settings->timesToShow())
+ nPrint = settings->numberOfEvents()/settings->timesToShow();
+ else
+ nPrint = 0;
+
+ unsigned long maxEvents = settings->numberOfEvents();
+
+ cout << "Generating " << maxEvents << " events." << endl << endl;
+
+ //
+ // Event generation
+ //
+ for (unsigned long iEvent = 0; iEvent < maxEvents; iEvent++) {
+ //
+ // Generate one event
+ //
+ Event *event = sartre.generateEvent();
+ if (nPrint && (iEvent+1)%nPrint == 0 && iEvent != 0) {
+ cout << "processed " << iEvent+1 << " events" << endl;
+ }
+
+ //
+ // If Sartre is run in UPC mode, half of the events needs to be
+ // rotated around and axis perpendicular to z:
+ //
+ if(settings->UPC() and settings->A()==settings->UPCA()){
+ randomlyReverseBeams(event);
+ }
+
+ //
+ // Print out (here only for the first few events)
+ //
+ if (iEvent < 10) event->list();
+
+ //
+ // Add up the cross-section
+ //
+
+ //
+ // Fill ROOT tree
+ //
+ myRootSartreEvent.iEvent = event->eventNumber;
+ myRootSartreEvent.t = event->t;
+ myRootSartreEvent.beta = event->beta;
+ myRootSartreEvent.Q2 = event->Q2;
+ myRootSartreEvent.x = event->x;
+ myRootSartreEvent.s = event->s;
+ myRootSartreEvent.y = event->y;
+ myRootSartreEvent.W = event->W;
+ myRootSartreEvent.z = event->z;
+ myRootSartreEvent.MX = event->MX;
+ myRootSartreEvent.xpom = event->xpom;
+ myRootSartreEvent.pol = event->polarization == transverse ? 0 : 1;
+ myRootSartreEvent.dmode = event->diffractiveMode == coherent ? 0 : 1;
+ eIn = &event->particles[0].p;
+ pIn = &event->particles[1].p;
+ eOut = &event->particles[2].p;
+ pOut = &event->particles[6].p;
+ X = &event->particles[4].p;
+ gamma = &event->particles[3].p;
+ quark = &event->particles[7].p;
+ antiquark= &event->particles[8].p;
+
+ //If the event is incoherent, and nuclear breakup is enabled, fill the remnants to the tree
+ if(settings->enableNuclearBreakup() and event->diffractiveMode == incoherent){
+ protons.Clear();
+ neutrons.Clear();
+ remnants.Clear();
+ for(unsigned int iParticle=7; iParticle < event->particles.size(); iParticle++){
+ if(event->particles[iParticle].status == 1) { // Final-state particle
+ const Particle& particle = event->particles[iParticle];
+ switch (abs(particle.pdgId)) {
+ case 2212: // (Anti-)proton
+ new(protons[protons.GetEntries()]) TLorentzVector(particle.p);
+ break;
+ case 2112: // (Anti-)neutron
+ new(neutrons[neutrons.GetEntries()]) TLorentzVector(particle.p);
+ break;
+ default: // Any other remnant
+ new(remnants[remnants.GetEntries()]) TLorentzVector(particle.p);
+ break;
+ } // switch
+ } // if
+ } // for
+ } // if
+ tree.Fill();
+ }
+ cout << "All events processed\n" << endl;
+
+ //
+ // That's it, finish up
+ //
+ cout<<"Calculating the total cross section on given interval."<<endl;
+ cout<<"This may take a while, use your radical freedom to do something else..."<<endl;
+ double totalCS=0;// #TTsartre.totalCrossSection();
+ TH1D* histoForCSandNumberOfEvents = new TH1D("histoForCSandNumberOfEvents", "Cross-section and Number of Events", 2, 0., 1.);
+ histoForCSandNumberOfEvents->SetBinContent(1, totalCS);
+ histoForCSandNumberOfEvents->SetBinContent(2, maxEvents);
+
+ double runTime = sartre.runTime();
+ hfile->Write();
+ cout << rootfile.c_str() << " written." << endl;
+ cout << "Total cross-section: " << totalCS << " nb" << endl;
+ sartre.listStatus();
+ cout << "CPU Time/event: " << runTime/maxEvents << " seconds/event" << endl;
+ return 0;
+}
+
+// UPC only
+void randomlyReverseBeams(Event* myEvent) {
+
+ TRandom3 *random = EventGeneratorSettings::randomGenerator();
+
+ if(random->Uniform(1) > 0.5){
+ for(unsigned int i=0; i<myEvent->particles.size(); i++)
+ myEvent->particles.at(i).p.RotateX(M_PI);
+ }
+}
+
+
Property changes on: trunk/examples/sartreInclusiveDiffractionMain.cpp
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property

File Metadata

Mime Type
text/x-diff
Expires
Sat, Dec 21, 2:14 PM (16 h, 15 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4023079
Default Alt Text
(26 KB)

Event Timeline