Page MenuHomeHEPForge
Files F7879090

No OneTemporary
Actions

View Options

Index: trunk/code/Makefile
===================================================================
--- trunk/code/Makefile (revision 473)
+++ trunk/code/Makefile (revision 474)
@@ -1,27 +1,71 @@
-all: jewel-2.1.0-vac jewel-2.1.0-simple jewel-2.1.1-vac jewel-2.1.1-simple
+all: jewel-2.3.0-vac jewel-2.3.0-simple jewel-2.3.0-hydro jewel-2.4.0-vac jewel-2.4.0-simple jewel-emmi-vac jewel-emmi-simple jewel-emmi-hydro jewel-230-hilmi-vac jewel-230-hilmi-simple jewel-240-hilmi-vac jewel-240-hilmi-simple
# path to LHAPDF library
-LHAPDF_PATH := /home/lhapdf/install/lib/
+LHAPDF_PATH :=/media/hdmobil/korinna/arbeit/durham/sherpa/lhapdf/lhapdf-5.8.4/lhapdf-install/lib/
FC := gfortran
-FFLAGS := -g -static -fno-align-commons
+FFLAGS := -g -static
jewel-2.1.0-vac: jewel-2.1.0.o medium-vac.o pythia6425mod.o
$(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
jewel-2.1.0-simple: jewel-2.1.0.o medium-simple.o pythia6425mod.o
$(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
jewel-2.1.1-vac: jewel-2.1.1.o medium-vac.o pythia6425mod.o
$(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
jewel-2.1.1-simple: jewel-2.1.1.o medium-simple.o pythia6425mod.o
$(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+jewel-2.2.0-vac: jewel-2.2.0.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.2.0-simple: jewel-2.2.0.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.3.0-vac: jewel-2.3.0.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.3.0-simple: jewel-2.3.0.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.3.0-hydro: jewel-2.3.0.o medium-hydro.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.4.0-vac: jewel-2.4.0.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-2.4.0-simple: jewel-2.4.0.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-emmi-vac: jewel-emmi.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-emmi-simple: jewel-emmi.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-emmi-hydro: jewel-emmi.o medium-hydro.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-230-hilmi-vac: jewel-230-hilmi.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-230-hilmi-simple: jewel-230-hilmi.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-240-hilmi-vac: jewel-240-hilmi.o medium-vac.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+jewel-240-hilmi-simple: jewel-240-hilmi.o medium-simple.o pythia6425mod.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF
+
+
+
clean:
rm -f medium-*.o
rm -f jewel*.o
rm -f pythia6425mod.o
rm -f *~
.PHONY: all
Index: trunk/code/generate-eventfile.sh
===================================================================
--- trunk/code/generate-eventfile.sh (revision 0)
+++ trunk/code/generate-eventfile.sh (revision 474)
@@ -0,0 +1,14 @@
+for ((j = 1000; j < 1010; j++))
+do
+ runcard=params.$j.dat
+ #jobfile=lhc-batch-$j.sh
+ echo $runcard
+ cp params_master.dat $runcard
+ #cp params.alice.dat $runcard
+ sed -i 's/xxxx/'$j'/g' $runcard
+ chmod u+x $runcard
+ #cp job_master.sh $jobfile
+ #sed -i 's/RUNCARD/'$runcard'/g' $jobfile
+ #chmod u+x $jobfile
+done
+
Property changes on: trunk/code/generate-eventfile.sh
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Index: trunk/code/JEWEL_BKGSUBTRACTION.cc
===================================================================
--- trunk/code/JEWEL_BKGSUBTRACTION.cc (revision 473)
+++ trunk/code/JEWEL_BKGSUBTRACTION.cc (revision 474)
@@ -1,481 +1,481 @@
//!-----------------------------------
/*
Authors: Raghav Kunnawalkam Elayavalli and Korinna Christine Zapp
July 5th 2017,
RIVET Analysis for Jet - Background subtraction in JEWEL
To be run with jewel-v2.2.0, with including recoils,
dummy particles and the scattering centers (HepMC:StatusCode 3)
*/
//!-----------------------------------
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "HepMC/GenParticle.h"
#include "HepMC/GenEvent.h"
#include "HepMC/GenCrossSection.h"
namespace Rivet {
class JEWEL_BKGSUBTRACTION : public Analysis {
public:
//! Constructor
JEWEL_BKGSUBTRACTION(string name = "JEWEL_BKGSUBTRACTION")
: Analysis(name)
{
setNeedsCrossSection(true);
}
//! Cells which make up the grid
struct candidate{
//! these are the boundaries
double etaMin;
double phiMin;
double etaMax;
double phiMax;
//! Pseudojets will have the objects inside
PseudoJets objects;
vector <int> jetID;
//! Four momeentum to get useful informations like eta, phi, mass, pT etc...
FourMomentum candMom;
FourMomentum bkgMom;
double sumnegpT;
double sumpT;
};
std::pair<int,int> FindCandidateBin(double eta, double phi){
int etaloc = -9;
int philoc = -9;
for(unsigned x = 0; x<_etabins.size()-1; ++x){
for(unsigned y = 0; y<_phibins.size()-1; ++y){
if(eta >= _etabins[x] && eta < _etabins[x+1] &&
phi >= _phibins[y] && phi < _phibins[y+1]){
etaloc = x;
philoc = y;
}
}
}
if(etaloc == -9 || philoc == -9)
return std::pair<int,int>(-1,-1);
else
return std::pair<int,int>(etaloc,philoc);
}
//! 4-Momentum Subtraction Procedure
PseudoJets do4MomSub(const Jets inputCollection,
vector <HepMC::GenParticle*> scatteringCenters,
bool doSubtraction){
PseudoJets jets_bkgsub;
//! loop over the input jet collection
foreach (const PseudoJet jet, inputCollection){
FourMomentum jt, bkg, jt_sub;
//! constituents map
foreach (const fastjet::PseudoJet &c, jet.constituents()){
jt += FourMomentum(c.E(), c.px(), c.py(), c.pz());
}
if(doSubtraction){
for(unsigned ip = 0; ip<scatteringCenters.size(); ++ip){
bool isinJet = false;
double _delRSmall = 100;
double _delRposition = 0;
int counter = 0;
foreach (const fastjet::PseudoJet &c, jet.constituents()){
const double _delRTruth = deltaR(scatteringCenters[ip]->momentum().eta(),
scatteringCenters[ip]->momentum().phi(),
c.eta(), c.phi_std());
if(_delRTruth < 1e-5){
_delRSmall = _delRTruth;
_delRposition = counter;
isinJet = true;
}
counter++;
}// jet constituents loop
// sum the background contribution
if(isinJet){
bkg += FourMomentum(scatteringCenters[ip]->momentum().e(),
scatteringCenters[ip]->momentum().px(),
scatteringCenters[ip]->momentum().py(),
scatteringCenters[ip]->momentum().pz());
}
}//! scattering center loop
jt_sub = jt - bkg;
}else
jt_sub = jt;
if(jt_sub.pT() <= _pTCut) continue;
jets_bkgsub.push_back(PseudoJet(jt_sub.px(), jt_sub.py(), jt_sub.pz(), jt_sub.E()));
}
if(jets_bkgsub.size()==0){
if(verbose) std::cout<<"no 4MomSub jets in the event with given kinematic cuts"<<std::endl;
}
return jets_bkgsub;
}
//! Grid Based Subtraction - 1 Procedure (Jet Clustering before Discretization)
PseudoJets doGridSub1(const Jets inputCollection,
vector <HepMC::GenParticle*> scatteringCenters,
bool doSubtraction){
//! initialize the grid
vector <vector<candidate> > grid;
for(int x = 0; x<_Nbounds_eta; ++x){
vector <candidate> xgrid;
for(int y = 0; y<_Nbounds_phi; ++y){
candidate test;
test.etaMin = _etabins[x];
test.phiMin = _phibins[y];
test.etaMax = _etabins[x+1];
test.phiMax = _phibins[y+1];
test.sumnegpT = 0.0;
xgrid.push_back(test);
}
grid.push_back(xgrid);
xgrid.clear();
}
//! Now loop over all the jets and fill their constituents into the grid
int jetcounter = 0;
foreach ( const PseudoJet jet, inputCollection ) {
jetcounter++;
//! loop over the candidates of the jet and include them in the grid
foreach ( const PseudoJet c, jet.constituents() ){
double ceta = c.eta();
double cphi = c.phi_std();
std::pair<int,int> candpos = FindCandidateBin(ceta, cphi);
if(candpos.first!=-1 && candpos.second!=-1){
grid[candpos.first][candpos.second].objects.push_back(c);
grid[candpos.first][candpos.second].candMom+=FourMomentum(c.E(), c.px(), c.py(), c.pz());
grid[candpos.first][candpos.second].jetID.push_back(jetcounter);
grid[candpos.first][candpos.second].sumpT+=c.pt();
}
}
}// jet loop
//! Now loop over the scattering centers that are near jets and add that as BKG
if(doSubtraction){
for(unsigned ip = 0; ip<scatteringCenters.size(); ++ip){
double px,py,pz,E;
px = scatteringCenters[ip]->momentum().px();
py = scatteringCenters[ip]->momentum().py();
pz = scatteringCenters[ip]->momentum().pz();
E = scatteringCenters[ip]->momentum().e();
PseudoJet part(px,py,pz,E);
bool iJ = false;
foreach(const PseudoJet j, inputCollection){
double delR = deltaR(j.eta(), j.phi_std(),
part.eta(), part.phi_std());
if(delR < _jetR){
iJ = true;
break;
}
}
if(iJ){
std::pair<int,int> candpos = FindCandidateBin(part.eta(), part.phi_std());
if(candpos.first!=-1 && candpos.second!=-1){
grid[candpos.first][candpos.second].objects.push_back(part);
grid[candpos.first][candpos.second].bkgMom+=FourMomentum(E, px, py, pz);
grid[candpos.first][candpos.second].sumpT-=part.pt();
}
}
}// scattering center loop
}// dosub
//! build up the pseudojets to do the clustering
PseudoJets pJet_sub;
for(int x = 0; x<_Nbounds_eta; ++x){
for(int y = 0; y<_Nbounds_phi; ++y){
FourMomentum jet_sub = grid[x][y].candMom;
if(doSubtraction)
jet_sub -= grid[x][y].bkgMom;
double px, py, pz, E;
if(doSubtraction && (grid[x][y].candMom.pT() - grid[x][y].bkgMom.pT() < 0.0)){
grid[x][y].sumnegpT += grid[x][y].bkgMom.pT() - grid[x][y].candMom.pT();
px = 0.0; py = 0.0; pz = 0.0; E = 0.0;
}else {
px = jet_sub.px();
py = jet_sub.py();
pz = jet_sub.pz();
E = jet_sub.E();
}
PseudoJet part_sub(px,py,pz,E);
if(part_sub.pt() != 0 && part_sub.E() != 0)
pJet_sub.push_back(part_sub);
}
}//!grid loop
fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, _jetR);
fastjet::ClusterSequence cs_sub(pJet_sub, jet_def);
PseudoJets jets_bkgsub = sorted_by_pt(cs_sub.inclusive_jets(_pTCut));
if(jets_bkgsub.size()==0){
if(verbose) std::cout<<"no GridSub1 jets in the event with given kinematic cuts"<<std::endl;
}
return jets_bkgsub;
}
//! Grid Based Subtraction - 2 Procedure (Discretization before jet clustering)
PseudoJets doGridSub2(const Event& event,
const Jets inputCollection,
vector <HepMC::GenParticle*> scatteringCenters,
bool doSubtraction){
//! initialize the grid
vector <vector<candidate> > grid;
for(int x = 0; x<_Nbounds_eta; ++x){
vector <candidate> xgrid;
for(int y = 0; y<_Nbounds_phi; ++y){
candidate test;
test.etaMin = _etabins[x];
test.phiMin = _phibins[y];
test.etaMax = _etabins[x+1];
test.phiMax = _phibins[y+1];
test.sumnegpT = 0.0;
xgrid.push_back(test);
}
grid.push_back(xgrid);
xgrid.clear();
}
//! Now loop over all the final state particles and fill their constituents into the grid
const ParticleVector& FS = applyProjection<FinalState>(event, "FS").particlesByPt();
foreach ( const Particle& p, FS) {
double peta = p.eta();
double pphi = p.phi(MINUSPI_PLUSPI);
double ppT = p.pT();
std::pair<int,int> candpos = FindCandidateBin(peta, pphi);
if(candpos.first!=-1 && candpos.second!=-1){
grid[candpos.first][candpos.second].objects.push_back(PseudoJet(p.px(), p.py(), p.pz(), p.E()));
grid[candpos.first][candpos.second].candMom+=FourMomentum(p.E(), p.px(), p.py(), p.pz());
grid[candpos.first][candpos.second].sumpT+=ppT;
}
}
//! This part of GridSub2 is same as GridSub1
//! Now loop over the scattering centers that are near jets and add that as BKG
if(doSubtraction){
for(unsigned ip = 0; ip<scatteringCenters.size(); ++ip){
double px,py,pz,E;
px = scatteringCenters[ip]->momentum().px();
py = scatteringCenters[ip]->momentum().py();
pz = scatteringCenters[ip]->momentum().pz();
E = scatteringCenters[ip]->momentum().e();
PseudoJet part(px,py,pz,E);
bool iJ = false;
foreach(const PseudoJet j, inputCollection){
double delR = deltaR(j.eta(), j.phi_std(),
part.eta(), part.phi_std());
if(delR < _jetR){
iJ = true;
break;
}
}
if(iJ){
std::pair<int,int> candpos = FindCandidateBin(part.eta(), part.phi_std());
if(candpos.first!=-1 && candpos.second!=-1){
grid[candpos.first][candpos.second].objects.push_back(part);
grid[candpos.first][candpos.second].bkgMom+=FourMomentum(E, px, py, pz);
grid[candpos.first][candpos.second].sumpT-=part.pt();
}
}
}// scattering center loop
}// dosub
//! build up the pseudojets to do the clustering
PseudoJets pJet_sub;
for(int x = 0; x<_Nbounds_eta; ++x){
for(int y = 0; y<_Nbounds_phi; ++y){
FourMomentum jet_sub = grid[x][y].candMom;
if(doSubtraction)
jet_sub -= grid[x][y].bkgMom;
double px, py, pz, E;
if(doSubtraction && (grid[x][y].candMom.pT() - grid[x][y].bkgMom.pT() < 0.0)){
grid[x][y].sumnegpT += grid[x][y].bkgMom.pT() - grid[x][y].candMom.pT();
px = 0.0; py = 0.0; pz = 0.0; E = 0.0;
}else {
px = jet_sub.px();
py = jet_sub.py();
pz = jet_sub.pz();
E = jet_sub.E();
}
PseudoJet part_sub(px,py,pz,E);
if(part_sub.pt() != 0 && part_sub.E() != 0)
pJet_sub.push_back(part_sub);
}
}//!grid loop
fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, _jetR);
fastjet::ClusterSequence cs_sub(pJet_sub, jet_def);
PseudoJets jets_bkgsub = sorted_by_pt(cs_sub.inclusive_jets(_pTCut));
if(jets_bkgsub.size()==0){
if(verbose) std::cout<<"no GridSub2 jets in the event with given kinematic cuts"<<std::endl;
}
return jets_bkgsub;
}
//! Book histograms and initialise projections before the run
void init() {
//! verbosity flag
verbose = false;
//! Jet Radius
_jetR=0.4;
//! Jet Kinematic cuts
_pTCut=100.0;
//! Grid initialization
_etaMax=2.4;
_phiMax=M_PI;
_delRMin=0.05;
//! cell boundaries
_Nbounds_eta = (int)2*_etaMax/(_delRMin);
_Nbounds_phi = (int)2*_phiMax/(_delRMin);
for(int i = 0; i<=_Nbounds_eta; ++i){
double etax = -1*_etaMax + i*_delRMin;
_etabins.push_back(etax);
}
for(int i = 0; i<=_Nbounds_phi; ++i){
double phix = -1*_phiMax + i*_delRMin;
_phibins.push_back(phix);
}
if(verbose)
std::cout<<"Grid we are using is ("<<_Nbounds_eta<<" x "<<_Nbounds_phi
<<") eta: [-"<<_etaMax<<", "<<_etaMax
<<"] phi: [-"<<_phiMax<<", "<<_phiMax<<"]"<<std::endl;
//! Initialize Histograms
//! Example jet pT histograms for 3 jet collections
_h_JetpT_noSub = bookHisto1D("JetpT_noSub", 30, 100, 400);
_h_JetpT_4MomSub = bookHisto1D("JetpT_4MomSub", 30, 100, 400);
_h_JetpT_GridSub1 = bookHisto1D("JetpT_GridSub1", 30, 100, 400);
_h_JetpT_GridSub2 = bookHisto1D("JetpT_GridSub2", 30, 100, 400);
//! final state and jet projection
FinalState fs(-3.0, 3.0, 0.*GeV);
addProjection(fs, "FS");
addProjection(FastJets(fs, FastJets::ANTIKT, _jetR), "Jets");
}
//! Perform the per-event analysis
void analyze(const Event& event) {
const double weight = event.weight();
//! **************************************
//! Jet Collection from all particles in the event
//! Used for w/o recoils and in vacuum
if(verbose) std::cout<<"Jet Collection built without subtraction"<<std::endl;
const FastJets& AJets = applyProjection<FastJets>(event, "Jets");
Cut cuts = Cuts::etaIn(-2.0, 2.0) & (Cuts::pT > _pTCut*GeV);
const Jets jets_noSub = AJets.jetsByPt(cuts);
foreach ( const PseudoJet jet, jets_noSub ) {
_h_JetpT_noSub->fill(jet.pt(), weight);
}
//! Checking if additional scattering centers are available from HepMC file
//! get the scattered particles from the gen event:
vector <HepMC::GenParticle*> pscat;
foreach (const HepMC::GenParticle* p, particles(event.genEvent())) {
if(p->status() == 3)
pscat.push_back((HepMC::GenParticle*)p);
}
bool doSubtraction = true;
if(pscat.size() == 0)
doSubtraction = false;
//! **************************************
//! BKG-Sub Jet Collection from all particles w/ recoils inclding the dummy and Scattering centers
if(verbose) std::cout<<"Jet Collection w/ 4MomSub Subtraction"<<std::endl;
PseudoJets jets_4MomSub = do4MomSub(jets_noSub, pscat, doSubtraction);
foreach ( const PseudoJet jet, jets_4MomSub ) {
_h_JetpT_4MomSub->fill(jet.pt(), weight);
}
//! **************************************
//! BKG-Sub Jet Collection from all particles w/ recoils inclding the Scattering centers
- if(verbose) std::cout<<"Jet Collection w/ GridSub1 Subtraction"<<std::endl;
- PseudoJets jets_GridSub1 = doGridSub1(jets_noSub, pscat, doSubtraction);
- foreach ( const PseudoJet jet, jets_GridSub1 ) {
- _h_JetpT_GridSub1->fill(jet.pt(), weight);
- }
+ //if(verbose) std::cout<<"Jet Collection w/ GridSub1 Subtraction"<<std::endl;
+ //PseudoJets jets_GridSub1 = doGridSub1(jets_noSub, pscat, doSubtraction);
+ //foreach ( const PseudoJet jet, jets_GridSub1 ) {
+ //_h_JetpT_GridSub1->fill(jet.pt(), weight);
+ //}
//! **************************************
//! BKG-Sub Jet Collection from all particles w/ recoils inclding the Scattering centers
- if(verbose) std::cout<<"Jet Collection w/ GridSub2 Subtraction"<<std::endl;
- PseudoJets jets_GridSub2 = doGridSub2(event, jets_noSub, pscat, doSubtraction);
- foreach ( const PseudoJet jet, jets_GridSub2 ) {
- _h_JetpT_GridSub2->fill(jet.pt(), weight);
- }
+ //if(verbose) std::cout<<"Jet Collection w/ GridSub2 Subtraction"<<std::endl;
+ //PseudoJets jets_GridSub2 = doGridSub2(event, jets_noSub, pscat, doSubtraction);
+ //foreach ( const PseudoJet jet, jets_GridSub2 ) {
+ //_h_JetpT_GridSub2->fill(jet.pt(), weight);
+ //}
}
//! Normalise histograms etc., after the run
void finalize() {
scale(_h_JetpT_noSub, crossSection()/picobarn/sumOfWeights());
scale(_h_JetpT_4MomSub, crossSection()/picobarn/sumOfWeights());
- scale(_h_JetpT_GridSub1, crossSection()/picobarn/sumOfWeights());
- scale(_h_JetpT_GridSub2, crossSection()/picobarn/sumOfWeights());
+ //scale(_h_JetpT_GridSub1, crossSection()/picobarn/sumOfWeights());
+ //scale(_h_JetpT_GridSub2, crossSection()/picobarn/sumOfWeights());
}
protected:
double _pTCut;
double _jetR;
bool verbose;
vector<double> _etabins;
vector<double> _phibins;
int _Nbounds_eta;
int _Nbounds_phi;
double _delRMin;
double _etaMax;
double _phiMax;
private:
Histo1DPtr _h_JetpT_noSub;
Histo1DPtr _h_JetpT_4MomSub;
Histo1DPtr _h_JetpT_GridSub1;
Histo1DPtr _h_JetpT_GridSub2;
};
//! The hook for the plugin system
DECLARE_RIVET_PLUGIN(JEWEL_BKGSUBTRACTION);
}
Index: trunk/code/runit-lhc5.sh
===================================================================
--- trunk/code/runit-lhc5.sh (revision 0)
+++ trunk/code/runit-lhc5.sh (revision 474)
@@ -0,0 +1,37 @@
+#!/bin/bash
+
+source /media/hdmobil/korinna/arbeit/rivet2/local/rivetenv.sh
+export RIVET_ANALYSIS_PATH=$RIVET_ANALYSIS_PATH:/media/hdmobil/korinna/arbeit/MC/jewel/trunc/code/
+#export RIVET_REF_PATH=$RIVET_REF_PATH:/home/jewel/trunk/code
+export LD_LIBRARY_PATH=/media/hdmobil/korinna/arbeit/durham/sherpa/lhapdf/lhapdf-5.8.4/lhapdf-install/lib/:$LD_LIBRARY_PATH
+export LHAPATH=/media/hdmobil/korinna/arbeit/durham/sherpa/lhapdf/lhapdf-5.8.4/lhapdf-install/share/lhapdf/PDFsets
+
+
+fifo=fifo-vac5.hepmc
+runcard=params.vac5.dat
+#outfile=data/jewel-2.3.0_vac_hl_5TeV_1M.yoda
+#outfile=data/jewel-2.1.1_vac_cms-zg.yoda
+outfile=data/times300.yoda
+
+rm $fifo
+mkfifo $fifo
+sed -i 's/FIFO/'$fifo'/g' $runcard
+#./jewel-2.3.0-vac $runcard &
+./jewel-emmi-vac $runcard &
+
+#rivet -a MC_XS -a MC_SOFTDROP_REC -a MC_SOFTDROP_NOREC_GRID -a MC_SOFTDROP_NOREC_PART -H $outfile $fifo
+#rivet -a MC_XS -a MC_SOFTDROP_REC -a MC_SOFTDROP_NOREC_GRID -H $outfile $fifo
+#rivet -a MC_XS -a MC_SOFTDROP_UNSUB_GRID -a MC_SOFTDROP_UNSUB_PART -H $outfile $fifo
+#rivet -a MC_XS -a MC_SOFTDROP_4MOMSUB -a MC_SOFTDROP_NOREC -H $outfile $fifo
+#rivet --pwd -a MC_XS -a MC_SOFTDROP_4MOMSUB -a MC_SOFTDROP_NOREC -a MC_SOFTDROP_REC -a MC_SOFTDROP_NOREC_GRID -H $outfile $fifo
+#rivet -a MC_XS -a MC_ETOT -H $outfile $fifo
+#rivet --pwd -a MC_XS -a TEST_JETMASS -H $outfile $fifo
+#rivet --pwd -a MC_XS -a CMS_SOFTDROP_4MOMSUB -H $outfile $fifo
+#rivet --pwd -a MC_XS -a MC_SOFTDROP_4MOMSUB -H $outfile $fifo
+#rivet --pwd -a MC_XS -a CMS_SOFTDROP_4MOMSUB -a CMS_GROOMEDMASS_4MOMSUB -a MC_JETMASS_TINY -a ATLAS_2018_I1673184 -H $outfile $fifo
+rivet --pwd -a MC_XS -a MC_TIMES -H $outfile $fifo
+
+
+sed -i 's/'$fifo'/FIFO/g' $runcard
+
+rm $fifo
Property changes on: trunk/code/runit-lhc5.sh
___________________________________________________________________
Added: svn:executable
## -0,0 +1 ##
+*
\ No newline at end of property
Index: trunk/code/analyse-lhc2-med.sh
===================================================================
--- trunk/code/analyse-lhc2-med.sh (revision 473)
+++ trunk/code/analyse-lhc2-med.sh (revision 474)
@@ -1,22 +1,26 @@
#!/bin/bash
#source /home/rivet/install/rivetenv.sh
#source /home/rivet2/local/rivetenv.sh
#export RIVET_ANALYSIS_PATH=$PWD
#export RIVET_ANALYSIS_PATH=$RIVET_ANALYSIS_PATH:$PWD
#export RIVET_REF_PATH=$PWD
-#for file in eventfiles/jewel-2.1.0_360MeV_rapfix_geo_md09_0-0_wo_20??.hepmc
-#for file in eventfiles/jewel-2.0.0_360MeV_nops_geo-0_md09_0-0_wo_10??.hepmc
+#for file in eventfiles/jewel-2.3.0_360MeV_jj_0-10_wrec_wors_100?.hepmc
+#for file in eventfiles/jewel-2.1.1_360MeV_rapfix_md09_0-10_w_100[4-9].hepmc
#for file in eventfiles/jewel-2.1.1_360MeV_massfix_geo_md09_0-10_wo_10??.hepmc
-#for file in eventfiles/jewel-2.1.0_360MeV_rapfix_geo_md09_0-10_w_pytest1_10??.hepm
-for file in eventfiles/jewel-2.1.1_360MeV_sc-soft_md09_0-10_wo_10*.hepmc
-#for file in eventfiles/jewel-2.1.0_360MeV_rapfix_geo_md09_0-10_wo_10??.hepmc
+#for file in eventfiles/jewel-2.1.1_360MeV_sc-soft_md09_0-10_wo_10*.hepmc
+#for file in eventfiles/jewel-2.3.0_2TeV_jj_km1_rm0m_0-10_wrec-0_wors_10[1-9]?.hepmc
+#for file in eventfiles/jewel-2.3.0_2TeV_q0_jj_km1_rm0r_0-10_wrec-0_wors_100?.hepmc
+for file in eventfiles/jewel-240-hilmi_2TeV_jj_km1_rm0m_0-10_wrec_100?.hepmc
do
outfile=${file/.hepmc/.yoda}
echo $outfile
- rivet --pwd -a MC_XS -a ALICE_2012_I1127497 -a ALICE_2012_I1210881 -a ALICE_2014_I1263194 -a ALICE_2015_I1343112 -a ALICE_2017_I1512107_4MOMSUB -a ALICE_CHPART -a ALICE_GIRTH_4MOMSUB -a ALICE_JETRAA2 -a ALICE_JETRAA -a ATLAS_2011_S8888116 -a ATLAS_2012_I1126965 -a ATLAS_2013_I1228693 -a ATLAS_2013_I1240088 -a ATLAS_2014_I1300152_4MOMSUB -a ATLAS_CONF_2011_075 -a ATLAS_CONF_2014_025 -a ATLAS_CONF_2014_028 -a CMS_2011_I889010 -a CMS_2012_I1088823 -a CMS_2012_I1090064 -a CMS_2012_I1116250 -a CMS_2013_I1256590_4MOMSUB -a CMS_CHPTSPEC -a CMS_HIN_12_004 -a HADRONSPECTRA -a JETSPECTRUM -a MC_JETS -H $outfile $file
+ rivet --pwd -a MC_XS -a MC_JETS_ISRAD -H $outfile $file
+ #rivet --pwd -a MC_XS -a CMS_2014_I1299142 -a ATLAS_2017_I1511869 -H $outfile $file
+ #rivet --pwd -a MC_XS -a ALICE_2017_I1512107_4MOMSUB -a ATLAS_2014_I1300152_4MOMSUB -a CMS_2013_I1256590_4MOMSUB -a CMS_2014_I1299142_4MOMSUB -a JEWEL_BKGSUBTRACTION -a JETSPECTRUM -H $outfile $file
+ #rivet --pwd -a MC_XS -a ALICE_2012_I1127497 -a ALICE_2012_I1210881 -a ALICE_2014_I1263194 -a ALICE_2015_I1343112 -a ALICE_2017_I1512107_4MOMSUB -a ALICE_CHPART -a ALICE_GIRTH_4MOMSUB -a ALICE_JETRAA2 -a ALICE_JETRAA -a ATLAS_2011_S8888116 -a ATLAS_2012_I1126965 -a ATLAS_2013_I1228693 -a ATLAS_2013_I1240088 -a ATLAS_2014_I1300152_4MOMSUB -a ATLAS_CONF_2011_075 -a ATLAS_CONF_2014_025 -a ATLAS_CONF_2014_028 -a CMS_2011_I889010 -a CMS_2012_I1088823 -a CMS_2012_I1090064 -a CMS_2012_I1116250 -a CMS_2013_I1256590_4MOMSUB -a CMS_CHPTSPEC -a CMS_HIN_12_004 -a HADRONSPECTRA -a JETSPECTRUM -a MC_JETS -H $outfile $file
#rivet --pwd -a MC_XS -H $outfile $file
#rivet -a MC_XS -a ALICE_2017_I1512107_4MOMSUB -a ALICE_2017_I1512107_GRIDSUB1 -a ALICE_2017_I1512107_GRIDSUB2 -a JEWEL_MethodComparisons_4MomSub -a JEWEL_MethodComparisons_GridSub -H $outfile $file
#rivet -a MC_XS -a JEWEL_MethodComparisons_4MomSub -a JEWEL_MethodComparisons_GridSub -H $outfile $file
done
Index: trunk/code/medium-simple.f
===================================================================
--- trunk/code/medium-simple.f (revision 473)
+++ trunk/code/medium-simple.f (revision 474)
@@ -1,774 +1,776 @@
SUBROUTINE MEDINIT(FILE,id,etam,mass)
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--longitudinal boost of momentum distribution
common/boostmed/boost
logical boost
C--factor to vary Debye mass
COMMON/MDFAC/MDFACTOR,MDSCALEFAC
DOUBLE PRECISION MDFACTOR,MDSCALEFAC
C--nuclear thickness function
COMMON /THICKFNC/ RMAX,TA(100,2)
DOUBLE PRECISION RMAX,TA
C--geometrical cross section
COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
DOUBLE PRECISION IMPMAX,CROSS
C--identifier of log file
common/logfile/logfid
integer logfid
DATA RAU/10./
DATA D3/0.9d0/
DATA ZETA3/1.2d0/
C--local variables
INTEGER I,LUN,POS,IOS,id,mass
double precision etam
CHARACTER*100 BUFFER,LABEL,tempbuf
CHARACTER*80 FILE
character firstchar
logical fileexist
etamax2 = etam
logfid = id
IOS=0
LUN=77
C--default settings
TAUI=0.6d0
TI=0.36d0
TC=0.17d0
WOODSSAXON=.TRUE.
CENTRMIN=0.d0
CENTRMAX=10.d0
NF=3
A=mass
N0=0.17d0
D=0.54d0
SIGMANN=6.2
MDFACTOR=0.45d0
MDSCALEFAC=0.9d0
boost = .true.
C--read settings from file
write(logfid,*)
inquire(file=FILE,exist=fileexist)
if(fileexist)then
write(logfid,*)'Reading medium parameters from ',FILE
OPEN(unit=LUN,file=FILE,status='old',err=10)
do 20 i=1,1000
READ(LUN, '(A)', iostat=ios) BUFFER
if (ios.ne.0) goto 30
firstchar = buffer(1:1)
if (firstchar.eq.'#') goto 20
POS=SCAN(BUFFER,' ')
LABEL=BUFFER(1:POS)
BUFFER=BUFFER(POS+1:)
IF (LABEL=="TAUI")THEN
READ(BUFFER,*,IOSTAT=IOS) TAUI
ELSE IF (LABEL=="TI") THEN
READ(BUFFER,*,IOSTAT=IOS) TI
ELSE IF (LABEL=="TC") THEN
READ(BUFFER,*,IOSTAT=IOS) TC
ELSE IF (LABEL=="WOODSSAXON") THEN
READ(BUFFER,*,IOSTAT=IOS) WOODSSAXON
ELSE IF (LABEL=="CENTRMIN") THEN
READ(BUFFER,*,IOSTAT=IOS) CENTRMIN
ELSE IF (LABEL=="CENTRMAX") THEN
READ(BUFFER,*,IOSTAT=IOS) CENTRMAX
ELSE IF (LABEL=="NF") THEN
READ(BUFFER,*,IOSTAT=IOS) NF
ELSE IF (LABEL=="N0") THEN
READ(BUFFER,*,IOSTAT=IOS) N0
ELSE IF (LABEL=="D") THEN
READ(BUFFER,*,IOSTAT=IOS) D
ELSE IF (LABEL=="SIGMANN") THEN
READ(BUFFER,*,IOSTAT=IOS) SIGMANN
ELSE IF (LABEL=="MDFACTOR") THEN
READ(BUFFER,*,IOSTAT=IOS) MDFACTOR
ELSE IF (LABEL=="MDSCALEFAC") THEN
READ(BUFFER,*,IOSTAT=IOS) MDSCALEFAC
else
write(logfid,*)'unknown label ',label
endif
20 continue
30 close(LUN,status='keep')
write(logfid,*)'...done'
goto 40
10 write(logfid,*)'Could not open medium parameter file, '//
& 'will run with default settings.'
else
write(logfid,*)'No medium parameter file found, '//
& 'will run with default settings.'
endif
40 write(logfid,*)'using parameters:'
write(logfid,*)'TAUI =',TAUI
write(logfid,*)'TI =',TI
write(logfid,*)'TC =',TC
write(logfid,*)'WOODSSAXON =',WOODSSAXON
write(logfid,*)'CENTRMIN =',CENTRMIN
write(logfid,*)'CENTRMAX =',CENTRMAX
write(logfid,*)'NF =',NF
write(logfid,*)'A =',A
write(logfid,*)'N0 =',N0
write(logfid,*)'D =',D
write(logfid,*)'SIGMANN =',SIGMANN
write(logfid,*)'MDFACTOR =',MDFACTOR
write(logfid,*)'MDSCALEFAC =',MDSCALEFAC
write(logfid,*)
write(logfid,*)
write(logfid,*)
C--calculate T_A(x,y)
CALL CALCTA
C--calculate geometrical cross section
CALL CALCXSECTION
END
SUBROUTINE MEDNEXTEVT
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--geometrical cross section
COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
DOUBLE PRECISION IMPMAX,CROSS
C--local variables
integer i,j
DOUBLE PRECISION PYR,R,b1,b2,gettemp
C--pick an impact parameter
r=(pyr(0)*(centrmax-centrmin)+centrmin)/100.
i=0
do 130 j=1,200
if ((r-cross(j,3)/cross(200,3)).ge.0.) then
i=i+1
else
goto 132
endif
130 continue
132 continue
b1 = (i-1)*0.1d0
b2 = i*0.1d0
breal = (b2*(cross(i,3)/cross(200,3)-r)
& +b1*(r-cross(i+1,3)/cross(200,3)))/
& (cross(i,3)/cross(200,3)-cross(i+1,3)/cross(200,3))
centr = r;
END
double precision function getcentrality()
implicit none
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
getcentrality=centr
end
SUBROUTINE PICKVTX(X,Y)
IMPLICIT NONE
DOUBLE PRECISION X,Y
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
C--local variables
DOUBLE PRECISION X1,X2,Y1,Y2,Z,XVAL,YVAL,ZVAL,NTHICK,PYR
X1=BREAL/2.-RAU
X2=RAU-BREAL/2.
Y1=-SQRT(4*RAU**2-BREAL**2)/2.
Y2=SQRT(4*RAU**2-BREAL**2)/2.
C--small system
! x1=-1.
! x2=1.
! y1=-1.
! y2=1.
C--end small system
131 XVAL=PYR(0)*(X2-X1)+X1
YVAL=PYR(0)*(Y2-Y1)+Y1
IF((NTHICK(XVAL-BREAL/2.,YVAL).EQ.0.d0).OR.
& NTHICK(XVAL+BREAL/2.,YVAL).EQ.0.d0) GOTO 131
C--small system
! if (sqrt(xval**2+yval**2).gt.1.d0) goto 131
C--end small system
ZVAL=PYR(0)*NTHICK(-BREAL/2.,0d0)*NTHICK(BREAL/2.,0d0)
Z=NTHICK(XVAL-BREAL/2.,YVAL)*NTHICK(XVAL+BREAL/2.,YVAL)
IF(ZVAL.GT.Z) GOTO 131
X=XVAL
Y=YVAL
END
SUBROUTINE SETB(BVAL)
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
DOUBLE PRECISION BVAL
BREAL=BVAL
END
SUBROUTINE GETSCATTERER(X,Y,Z,T,TYPE,PX,PY,PZ,E,MS)
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
C--internal medium parameters
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--longitudinal boost of momentum distribution
common/boostmed/boost
logical boost
C--function calls
DOUBLE PRECISION GETTEMP,GETMD,GETMOM,GETMS
C--identifier of log file
common/logfile/logfid
integer logfid
C--local variables
DOUBLE PRECISION X,Y,Z,T,MS,PX,PY,PZ,E,MD,TEMP
INTEGER TYPE
DOUBLE PRECISION R,PYR,pmax,wt,tau,theta,phi,pi,p,ys,pz2,e2
DATA PI/3.141592653589793d0/
R=PYR(0)
IF(R.LT.(2.*12.*NF*D3/3.)/(2.*12.*NF*D3/3.+3.*16.*ZETA3/2.))THEN
TYPE=2
+ MS=GETMS(X,Y,Z,T)
ELSE
TYPE=21
+ MS=GETMD(X,Y,Z,T)
ENDIF
- MS=GETMS(X,Y,Z,T)
- MD=GETMD(X,Y,Z,T)
TEMP=GETTEMP(X,Y,Z,T)
tau=sqrt(t**2-z**2)
if (boost) then
ys = 0.5*log((t+z)/(t-z))
else
ys = 0.d0
endif
pmax = 10.*temp
IF(TEMP.LT.1.D-2)THEN
write(logfid,*)'asking for a scattering centre without medium:'
write(logfid,*)'at (x,y,z,t)=',X,Y,Z,T
write(logfid,*)'making one up to continue but '//
& 'something is wrong!'
TYPE=21
PX=0.d0
PY=0.d0
PZ=0.d0
MS=GETMS(0.d0,0.d0,0.d0,0.d0)
MD=GETMD(0.d0,0.d0,0.d0,0.d0)
E=SQRT(PX**2+PY**2+PZ**2+MS**2)
RETURN
ENDIF
10 p = pyr(0)**0.3333333*pmax
E2 = sqrt(p**2+ms**2)
if (type.eq.2) then
wt = (exp(ms/temp)-1.)/(exp(E2/temp)-1.)
else
wt = (exp(ms/temp)+1.)/(exp(E2/temp)+1.)
endif
if (wt.gt.1.) write(logfid,*)'Error in getscatterer: weight = ',wt
if (wt.lt.0.) write(logfid,*)'Error in getscatterer: weight = ',wt
if (pyr(0).gt.wt) goto 10
phi = pyr(0)*2.*pi
theta = -acos(2.*pyr(0)-1.)+pi
px = p*sin(theta)*cos(phi)
py = p*sin(theta)*sin(phi)
pz2 = p*cos(theta)
E = cosh(ys)*E2 + sinh(ys)*pz2
pz = sinh(ys)*E2 + cosh(ys)*pz2
END
SUBROUTINE AVSCATCEN(X,Y,Z,T,PX,PY,PZ,E,m)
IMPLICIT NONE
C--longitudinal boost of momentum distribution
common/boostmed/boost
logical boost
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--local variables
double precision x,y,z,t,px,py,pz,e,getms,m,ys
if (boost) then
ys = 0.5*log((t+z)/(t-z))
if ((z.eq.0.d0).and.(t.eq.0.d0)) ys =0.d0
if (ys.gt.etamax2) ys=etamax2
if (ys.lt.-etamax2) ys=-etamax2
else
ys = 0.d0
endif
m = getms(x,y,z,t)
e = m*cosh(ys)
px = 0.d0
py = 0.d0
pz = m*sinh(ys)
end
SUBROUTINE maxscatcen(PX,PY,PZ,E,m)
IMPLICIT NONE
C--longitudinal boost of momentum distribution
common/boostmed/boost
logical boost
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--local variables
double precision px,py,pz,e,getmsmax,m,ys
if (boost) then
ys = etamax2
else
ys = 0.d0
endif
m = getmsmax()
e = m*cosh(ys)
px = 0.d0
py = 0.d0
pz = m*sinh(ys)
end
DOUBLE PRECISION FUNCTION GETMD(X1,Y1,Z1,T1)
IMPLICIT NONE
C--factor to vary Debye mass
COMMON/MDFAC/MDFACTOR,MDSCALEFAC
DOUBLE PRECISION MDFACTOR,MDSCALEFAC
DOUBLE PRECISION X1,Y1,Z1,T1,GETTEMP
GETMD=MDSCALEFAC*3.*GETTEMP(X1,Y1,Z1,T1)
GETMD=MAX(GETMD,MDFACTOR)
END
DOUBLE PRECISION FUNCTION GETMS(X2,Y2,Z2,T2)
IMPLICIT NONE
DOUBLE PRECISION X2,Y2,Z2,T2,GETMD
GETMS=GETMD(X2,Y2,Z2,T2)/SQRT(2.)
END
DOUBLE PRECISION FUNCTION GETNEFF(X3,Y3,Z3,T3)
IMPLICIT NONE
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C-- local variables
DOUBLE PRECISION X3,Y3,Z3,T3,PI,GETTEMP,tau,cosheta
DATA PI/3.141592653589793d0/
tau = sqrt(t3**2-z3**2)
cosheta = t3/tau
GETNEFF=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
& *GETTEMP(X3,Y3,Z3,T3)**3/PI**2
getneff = getneff/cosheta
C--small system
! getneff = getneff/3.
C--end small system
END
DOUBLE PRECISION FUNCTION GETTEMP(X4,Y4,Z4,T4)
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--local variables
DOUBLE PRECISION X4,Y4,Z4,T4,TAU,NPART,EPS0,EPSIN,TEMPIN,PI,
&NTHICK,ys
DATA PI/3.141592653589793d0/
GETTEMP=0.D0
C--small system
! if (sqrt(x4**2+y4**2).gt.1.) return
C--end small system
-
+
IF(ABS(Z4).GT.T4)RETURN
TAU=SQRT(T4**2-Z4**2)
C--check for overlap region
IF((NTHICK(X4-BREAL/2.,Y4).EQ.0.d0).OR.
&NTHICK(X4+BREAL/2.,Y4).EQ.0.d0) RETURN
ys = 0.5*log((t4+z4)/(t4-z4))
if (abs(ys).gt.etamax2) return
C--determine initial temperature at transverse position
IF(WOODSSAXON)THEN
EPS0=(16.*8.+7.*2.*6.*NF)*PI**2*TI**4/240.
EPSIN=EPS0*NPART(X4-BREAL/2.,Y4,X4+BREAL/2.,Y4)
& *PI*RAU**2/(2.*A)
+! EPSIN=EPS0*NPART(X4-BREAL/2.,Y4,X4+BREAL/2.,Y4)/
+! & NPART(0.d0,0.d0,0.d0,0.d0)
TEMPIN=(EPSIN*240./(PI**2*(16.*8.+7.*2.*6.*NF)))**0.25
ELSE
TEMPIN=TI
ENDIF
C--calculate temperature if before initial time
IF(TAU.LE.TAUI)THEN
GETTEMP=TEMPIN*TAU/TAUI
ELSE
C--evolve temperature
GETTEMP=TEMPIN*(TAUI/TAU)**0.3333
ENDIF
IF(GETTEMP.LT.TC) GETTEMP=0.d0
END
DOUBLE PRECISION FUNCTION GETTEMPMAX()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--function call
DOUBLE PRECISION GETTEMP
GETTEMPMAX=GETTEMP(0.D0,0.D0,0.D0,TAUI)
END
DOUBLE PRECISION FUNCTION GETMDMAX()
IMPLICIT NONE
C--factor to vary Debye mass
COMMON/MDFAC/MDFACTOR,MDSCALEFAC
DOUBLE PRECISION MDFACTOR,MDSCALEFAC
DOUBLE PRECISION GETTEMPMAX
GETMDMAX=MDSCALEFAC*3.*GETTEMPMAX()
GETMDMAX=MAX(GETMDMAX,MDFACTOR)
END
DOUBLE PRECISION FUNCTION GETMDMIN()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--factor to vary Debye mass
COMMON/MDFAC/MDFACTOR,MDSCALEFAC
DOUBLE PRECISION MDFACTOR,MDSCALEFAC
DOUBLE PRECISION GETTEMPMAX
GETMDMIN=MDSCALEFAC*3.*TC
GETMDMIN=MAX(GETMDMIN,MDFACTOR)
END
DOUBLE PRECISION FUNCTION GETMSMAX()
IMPLICIT NONE
DOUBLE PRECISION GETMDMAX,SQRT
GETMSMAX=GETMDMAX()/SQRT(2.D0)
END
DOUBLE PRECISION FUNCTION GETNATMDMIN()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--factor to vary Debye mass
COMMON/MDFAC/MDFACTOR,MDSCALEFAC
DOUBLE PRECISION MDFACTOR,MDSCALEFAC,PI
DATA PI/3.141592653589793d0/
C--local variables
DOUBLE PRECISION T,GETMDMIN
T=GETMDMIN()/(MDSCALEFAC*3.)
GETNATMDMIN=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
& *T**3/PI**2
END
DOUBLE PRECISION FUNCTION GETLTIMEMAX()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C--function call
DOUBLE PRECISION GETTEMPMAX
GETLTIMEMAX=TAUI*(GETTEMPMAX()/TC)**3*cosh(etamax2)
END
DOUBLE PRECISION FUNCTION GETNEFFMAX()
IMPLICIT NONE
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--max rapidity
common/rapmax2/etamax2
double precision etamax2
C-- local variables
DOUBLE PRECISION PI,GETTEMPMAX
DATA PI/3.141592653589793d0/
GETNEFFMAX=(2.*6.*NF*D3*2./3. + 16.*ZETA3*3./2.)
& *GETTEMPMAX()**3/PI**2
END
DOUBLE PRECISION FUNCTION NPART(XX1,YY1,XX2,YY2)
IMPLICIT NONE
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--local variables
DOUBLE PRECISION XX1,YY1,XX2,YY2,NTHICK
NPART = NTHICK(XX1,YY1)*(1.-EXP(-SIGMANN*NTHICK(XX2,YY2))) +
& NTHICK(XX2,YY2)*(1.-EXP(-SIGMANN*NTHICK(XX1,YY1)))
END
DOUBLE PRECISION FUNCTION NTHICK(X1,Y1)
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--identifier of log file
common/logfile/logfid
integer logfid
C--nuclear thickness function
COMMON /THICKFNC/ RMAX,TA(100,2)
DOUBLE PRECISION RMAX,TA
INTEGER LINE,LMIN,LMAX,I
DOUBLE PRECISION X1,Y1,XA(4),YA(4),Y,DY,R,C,B,DELTA
R=SQRT(X1**2+Y1**2)
IF(R.GT.TA(100,1))THEN
NTHICK=0.
ELSE
LINE=INT(R*99.d0/TA(100,1)+1)
LMIN=MAX(LINE,1)
LMIN=MIN(LMIN,99)
IF((R.LT.TA(LMIN,1)).OR.(R.GT.TA(LMIN+1,1)))
& write(logfid,*)LINE,LMIN,R,TA(LMIN,1),TA(LMIN+1,1)
XA(1)=TA(LMIN,1)
XA(2)=TA(LMIN+1,1)
YA(1)=TA(LMIN,2)
YA(2)=TA(LMIN+1,2)
C=(YA(2)-YA(1))/(XA(2)-XA(1))
B=YA(1)-C*XA(1)
NTHICK=C*R+B
ENDIF
END
SUBROUTINE CALCTA()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C-- nuclear thickness function
COMMON /THICKFNC/ RMAX,TA(100,2)
DOUBLE PRECISION RMAX,TA
C--variables for integration
COMMON/INTEG/B,R
DOUBLE PRECISION B,R
C--local variables
INTEGER NSTEPS,I
DOUBLE PRECISION EPS,HFIRST,Y
NSTEPS=100
EPS=1.E-4
HFIRST=0.1D0
R=1.12*A**(0.33333)-0.86*A**(-0.33333)
RMAX=2.*R
DO 10 I=1,NSTEPS
C--set transverse position
B=(I-1)*2.D0*R/NSTEPS
Y=0.D0
C--integrate along longitudinal line
CALL ODEINT(Y,-2*R,2*R,EPS,HFIRST,0.d0,101)
TA(I,1)=B
TA(I,2)=Y
10 CONTINUE
END
SUBROUTINE CALCXSECTION()
IMPLICIT NONE
C--medium parameters
COMMON/MEDPARAM/CENTRMIN,CENTRMAX,BREAL,CENTR,RAU,NF
INTEGER NF
DOUBLE PRECISION CENTRMIN,CENTRMAX,BREAL,CENTR,RAU
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C-- geometrical cross section
COMMON /CROSSSEC/ IMPMAX,CROSS(200,3)
DOUBLE PRECISION IMPMAX,CROSS
C--local variables
INTEGER IX,IY,IB
DOUBLE PRECISION B,P,PROD,X,Y,NTHICK,NPART,pprev
pprev=0.
DO 30 IB=1,200
B=0.1d0*IB
PROD=1.d0
DO 10 IX=1,100
DO 20 IY=1,100
X=-20.d0+IX*0.4d0
Y=-20.d0+IY*0.4d0
PROD=PROD*
&EXP(-NTHICK(X+B/2.D0,Y)*SIGMANN)**(0.16d0*NTHICK(X-B/2.D0,Y))
20 CONTINUE
10 CONTINUE
P=(1.D0-PROD)*8.8D0/14.D0*B
CROSS(IB,1)=B
CROSS(IB,2)=P
if (ib.eq.1) then
cross(ib,3)=0.
else
cross(ib,3)=cross(ib-1,3)+(p+pprev)/2.*0.1
endif
pprev=p
30 CONTINUE
IMPMAX=19.95
END
DOUBLE PRECISION FUNCTION MEDDERIV(XVAL,W)
IMPLICIT NONE
DOUBLE PRECISION XVAL
INTEGER W
C--medium parameters
COMMON/MEDPARAMINT/TAUI,TI,TC,D3,ZETA3,D,
&N0,SIGMANN,A,WOODSSAXON
DOUBLE PRECISION TAUI,TI,TC,ALPHA,BETA,GAMMA,D3,ZETA3,D,N0,
&SIGMANN
INTEGER A
LOGICAL WOODSSAXON
C--variables for integration
COMMON/INTEG/B,R
DOUBLE PRECISION B,R
IF (W.EQ.1) THEN
C--XVAL corresponds to z-coordinate
MEDDERIV=N0/(1+EXP((SQRT(B**2+XVAL**2)-R)/D))
ELSE
MEDDERIV=0.D0
ENDIF
END
Index: trunk/code/Makefile-lhapdf6
===================================================================
--- trunk/code/Makefile-lhapdf6 (revision 0)
+++ trunk/code/Makefile-lhapdf6 (revision 474)
@@ -0,0 +1,21 @@
+all: jewel-2.2.0-vac jewel-2.2.0-simple
+
+# path to LHAPDF library
+LHAPDF_PATH := /media/hdmobil/korinna/arbeit/lhapdf6/lib
+
+FC := gfortran
+FFLAGS := -g -static
+
+jewel-2.2.0-vac: jewel-2.2.0.o medium-vac.o pythia6425mod.o meix.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF -lstdc++
+
+jewel-2.2.0-simple: jewel-2.2.0.o medium-simple.o pythia6425mod.o meix.o
+ $(FC) -o $@ -L$(LHAPDF_PATH) $^ -lLHAPDF -lstdc++
+
+clean:
+ rm -f medium-*.o
+ rm -f jewel*.o
+ rm -f pythia6425mod.o meix.o
+ rm -f *~
+
+.PHONY: all
Index: trunk/code/runit-lhc2.sh
===================================================================
--- trunk/code/runit-lhc2.sh (revision 473)
+++ trunk/code/runit-lhc2.sh (revision 474)
@@ -1,22 +1,28 @@
#!/bin/bash
-#source /media/hdmobil/korinna/arbeit/rivet2/local/rivetenv.sh
+source /media/hdmobil/korinna/arbeit/rivet2/local/rivetenv.sh
#export RIVET_ANALYSIS_PATH=$RIVET_ANALYSIS_PATH:/home/jewel/trunk/code
#export RIVET_REF_PATH=$RIVET_REF_PATH:/home/jewel/trunk/code
-export LD_LIBRARY_PATH=/home/lhapdf/install/lib/:$LD_LIBRARY_PATH
-export LHAPATH=/home/lhapdf/install/share/lhapdf/PDFsets
+export LD_LIBRARY_PATH=/media/hdmobil/korinna/arbeit/durham/sherpa/lhapdf/lhapdf-5.8.4/lhapdf-install/lib/:$LD_LIBRARY_PATH
+export LHAPATH=/media/hdmobil/korinna/arbeit/durham/sherpa/lhapdf/lhapdf-5.8.4/lhapdf-install/share/lhapdf/PDFsets
-fifo=fifo-vac2.hepmc
+fifo=fifo-vac.hepmc
runcard=params.vac.dat
-#outfile=test.yoda
-outfile=data/jewel-2.1.1_vac-276TeV_q15.aida
+outfile=data/jetspec.yoda
rm $fifo
mkfifo $fifo
sed -i 's/FIFO/'$fifo'/g' $runcard
-./jewel-2.1.1-vac $runcard &
- rivet --pwd -a MC_XS -a ALICE_2012_I1127497 -a ALICE_2012_I1210881 -a ALICE_2014_I1263194 -a ALICE_2015_I1343112 -a ALICE_2017_I1512107_4MOMSUB -a ALICE_CHPART -a ALICE_GIRTH_4MOMSUB -a ALICE_JETRAA2 -a ALICE_JETRAA -a ATLAS_2011_S8888116 -a ATLAS_2012_I1126965 -a ATLAS_2013_I1228693 -a ATLAS_2013_I1240088 -a ATLAS_2014_I1300152_4MOMSUB -a ATLAS_CONF_2011_075 -a ATLAS_CONF_2014_025 -a ATLAS_CONF_2014_028 -a CMS_2011_I889010 -a CMS_2012_I1088823 -a CMS_2012_I1090064 -a CMS_2012_I1116250 -a CMS_2013_I1256590_4MOMSUB -a CMS_CHPTSPEC -a CMS_HIN_12_004 -a HADRONSPECTRA -a JETSPECTRUM -a MC_JETS -H $outfile $fifo
+./jewel-2.4.0-vac $runcard &
+ rivet --pwd -a MC_XS -aMC_JETSPEC -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a CMS_2014_I1299142 -a ATLAS_2017_I1511869 -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a CMS_GROOMEDMASS_4MOMSUB -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a MC_JETMASS_TINY -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a ALICE_2017_I1512107_CONSTSUB -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a ALICE_2017_I1512107_4MOMSUB -a CMS_2013_I1256590_4MOMSUB -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a ALICE_2017_I1512107_4MOMSUB -a ATLAS_2014_I1300152_4MOMSUB -a CMS_2013_I1256590_4MOMSUB -a CMS_2014_I1299142_4MOMSUB -a JEWEL_BKGSUBTRACTION -a JETSPECTRUM -H $outfile $fifo
+ #rivet --pwd -a MC_XS -a ALICE_2012_I1127497 -a ALICE_2012_I1210881 -a ALICE_2014_I1263194 -a ALICE_2015_I1343112 -a ALICE_2017_I1512107_4MOMSUB -a ALICE_CHPART -a ALICE_GIRTH_4MOMSUB -a ALICE_JETRAA2 -a ALICE_JETRAA -a ATLAS_2011_S8888116 -a ATLAS_2012_I1126965 -a ATLAS_2013_I1228693 -a ATLAS_2013_I1240088 -a ATLAS_2014_I1300152_4MOMSUB -a ATLAS_CONF_2011_075 -a ATLAS_CONF_2014_025 -a ATLAS_CONF_2014_028 -a CMS_2011_I889010 -a CMS_2012_I1088823 -a CMS_2012_I1090064 -a CMS_2012_I1116250 -a CMS_2013_I1256590_4MOMSUB -a CMS_CHPTSPEC -a CMS_HIN_12_004 -a HADRONSPECTRA -a JETSPECTRUM -a MC_JETS -H $outfile $fifo
sed -i 's/'$fifo'/FIFO/g' $runcard
rm $fifo

File Metadata

Mime Type
text/x-diff
Expires
Tue, Nov 19, 7:31 PM (1 d, 9 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
3805841
Default Alt Text
(50 KB)

Event Timeline