IAP GITLAB

Commit 89710036 authored by Colin Baus's avatar Colin Baus

added example analyser

git-svn-id: https://devel-ik.fzk.de/svn/mc/crmc/trunk@4630 c7a5e08c-de06-0410-9364-b41cf42a0b17
parent a3b3ffba
SYSTEM=$(shell uname)
SYSTEM2=$(shell uname -m)
m32=$(shell if [ "$(SYSTEM2)" = "i686" ]; then echo "-m32"; fi)
SOURCE_DIR=src/
VER=199
LIBDIR=obj/
CFLAGS = $(m32) -ggdb3 -pipe -Wall
CXXFLAGS = $(m32) -ggdb3 -pipe -Wall
CXXFLAGS += -Wpacked -malign-double -mpreferred-stack-boundary=8
ifndef CXX
CXX = g++
endif
ifndef CC
CC = gcc
endif
ROOTCFLAGS = $(shell root-config --cflags)
ROOTLIBS = $(shell root-config --libs) -lEG
HEPCFLAGS = -I$(HEP_ROOT)/include
HEPLIBS = -L$(HEP_ROOT)/lib -lHepMC -lHepMCfio -lboost_iostreams
CFLAGS = $(ROOTCFLAGS) $(HEPCFLAGS)
LIBS = $(ROOTLIBS) $(HEPLIBS)
CXXFILES=analysis.cc
OBJS=$(FILES:%.f=$(LIBDIR)%.o)
CXXOBJS=$(CXXFILES:%.cc=$(LIBDIR)%.o)
all: check dirs bin_dir bin/analyse
.PHONY : check
check:
@if [ -z "$(HEP_ROOT)" ]; then echo "Please set HEP_ROOT to the root directory of HepMC2"; exit 1; fi
@if [ -z "$(ROOTSYS)" ]; then echo "Please set ROOTSYS to the root directory of root"; exit 1; fi
bin/analyse: $(OBJS) $(CXXOBJS)
$(CXX) $(CXXFLAGS) $(OBJS) $(CXXOBJS) -o $@ $(LIBS)
@(echo "")
@(echo "==> compilation successful!")
@(echo "==> type bin/analyse")
@(echo "")
$(OBJS) : $(LIBDIR)%.o : $(SOURCE_DIR)%.f
$(FC) -o $@ -c $<
$(CXXOBJS) : $(LIBDIR)%.o : $(SOURCE_DIR)%.cc src/analysis.h
$(CXX) $(CXXFLAGS) $(CFLAGS) -o $@ -c $<
dirs:
@if [ ! -d $(LIBDIR) ] ;then \
set -x; mkdir -p $(LIBDIR); set +x;\
fi
bin_dir:
@if [ ! -d bin ] ;then \
set -x; mkdir -p bin; set +x; \
fi
clean:
rm -rf $(LIBDIR)
rm -rf bin
rm -f *.copy *.check fort*
rm -f IO_GenEvent.dat
rm -f *~ \#*\# */*~
#;;; Emacs Setup
#;;; Local Variables:
#;;; mode: Makefile
#;;; End:
This diff is collapsed.
#include <algorithm>
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <string>
#include <sstream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#include <boost/iostreams/filter/zlib.hpp>
#include <boost/filesystem/operations.hpp>
#include <boost/filesystem/path.hpp>
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/GenParticle.h"
#include <HepMC/HeavyIon.h>
#include "HepMC/HEPEVT_Wrapper.h"
#include "HepMC/IO_GenEvent.h"
#include "HepMC/Units.h"
#include "HepMC/GenEvent.h"
#include "analysis.h"
namespace io = boost::iostreams;
using namespace std;
int main (int argc, char **argv)
{
//-------------------SET UP DATA
TFile* theOutFile;
string outFileName("new_histogram_file.root");
cout << " ! Opening output file: " << outFileName << endl;
theOutFile = new TFile(outFileName.c_str(),"RECREATE");
vector<string> filesModel1;
filesModel1.push_back("files/test.hepmc"); //add your files here
DataManager data;
data.SetFiles(filesModel1); //for more models, loop over models and call SetFiles each time
theOutFile->mkdir("model1");
//-------------------EVENT LOOP
while (data.GetNextEvent())
{
HepMC::HeavyIon* heavyIonInfo = NULL; //helpful to get cross section: heavyIonInfo->sigma_inel_NN ()
heavyIonInfo = data.evt->heavy_ion ();
//-------------------PARTICLE LOOP
HepMC::GenEvent::particle_const_iterator par = data.evt->particles_begin ();
for (; par != data.evt->particles_end (); ++par)
{
HepMC::GenParticle* p = (*par);
if (p->status () != 1) continue; //get final state particles status==2 would be decayed particles, status == 4 is beam particles
HepMC::GenVertex* parent_vertex = p->production_vertex();
const int id = p->pdg_id ();
const double pt = p->momentum ().perp ();
const double e = p->momentum ().e ();
const double eta = p->momentum ().eta ();
//for more advance paramters see HepMC documentation or load #include <TParticle.h> and fill object (see analysis.h)
//-------------------EVENT SELECTION
}//PARTICLE LOOP
//fill histograms.......
}//EVENT LOOP
std::cout << " ! Writing output file" << std::endl;
theOutFile->Write();
std::cout << " ! Closing output file" << std::endl;
theOutFile->Close();
delete theOutFile;
return 0;
}
#ifndef __analysis__h__
#define __analysis__h__
#include <string>
#include <map>
#include <vector>
#include <TFile.h>
#include <TGraphErrors.h>
#include <TGraphAsymmErrors.h>
#include <stdexcept>
#include <exception>
#include "HepMC/IO_GenEvent.h"
#include <boost/iostreams/device/file_descriptor.hpp>
#include <boost/iostreams/filtering_stream.hpp>
#include <boost/iostreams/filter/gzip.hpp>
#include <boost/iostreams/filter/zlib.hpp>
// #include <TParticle.h>
// #include <TParticlePDG.h>
// TParticle* CopyHepMC2ROOT(HepMC::GenParticle* a);
class DataManager
{
private:
std::vector<std::string> filelist;
std::vector<std::string>::iterator current_file;
boost::iostreams::filtering_istream in;
HepMC::IO_GenEvent* ascii_in;
int count;
bool GetNextFile()
{
if(++current_file == filelist.end())
return 0;
return ReadFile();
}
bool ReadFile()
{
std::string filename = *current_file;
bool use_compression = false;
if((filename).find(".gz") != std::string::npos) use_compression = true;
std::cout << "DataMangager::Opening file " << (filename) << " with" << (use_compression?" ":"out ") << "compression. " << std::endl;
in.reset();
if(use_compression)
in.push (boost::iostreams::gzip_decompressor ());
in.push (boost::iostreams::file_descriptor_source (filename));
if(ascii_in) delete ascii_in;
ascii_in = new HepMC::IO_GenEvent(in);
return true;
}
public:
HepMC::GenEvent* evt;
DataManager() : ascii_in(0), count(0), evt(0)
{
}
~DataManager()
{
if(ascii_in) delete ascii_in;
std::cout << "DataManger::Closing after reading " << count << " events." << std::endl;
}
void SetFiles(const std::vector<std::string>& files)
{
filelist = files; //copy
//open first file
current_file = filelist.begin();
ReadFile();
}
bool GetNextEvent()
{
delete evt;
evt = 0;
try
{
(*ascii_in) >> evt;
}
catch (std::exception& e)
{
evt = false;
std::cout << "Event could not be read: " << e.what() << std::endl;
}
if (!evt)
{
if (!GetNextFile()) return false;
else // new file selected
{
delete evt;
evt = NULL;
try
{
(*ascii_in) >> evt;
}
catch(std::exception& e)
{
evt = false;
std::cout << "Event could not be read again: " << e.what() << std::endl;
}
if (!evt) return false;
}
}
if(count % 1000 == 0)
std::cout << "Reading event: " << count << std::endl;
count++;
return true;
}
};
// TParticle* CopyHepMC2ROOT(HepMC::GenParticle* a)
// {
// TParticle* b = new TParticle;
// if(!b)
// return 0;
// b->SetPdgCode(a->pdg_id ());
// b->SetStatusCode(a->status ());
// b->SetMomentum(a->momentum ().px (),
// a->momentum ().py (),
// a->momentum ().pz (),
// a->momentum ().e ()
// );
// b->SetProductionVertex(a->production_vertex ()->position ().x (),
// a->production_vertex ()->position ().y (),
// a->production_vertex ()->position ().z (),
// a->production_vertex ()->position ().t ()
// );
// return b;
// }
#endif
......@@ -37,8 +37,8 @@ Reference : C. Baus, T. Pierog and R. Ulrich. To be published (2013)
EPOS LHC (-m0) : T. Pierog et al.
arXiv:1306.0121 [hep-ph]
QGSJETII-04 (-m7) : S.Ostapchenko,
Phys.Rev. D83 (2011) 014018
QGSJETII-04 (-m7) : S.Ostapchenko,
Phys.Rev. D83 (2011) 014018
* Pre LHC :
......@@ -48,7 +48,7 @@ EPOS 1.99 (-m1) : K. Werner, F.M. Liu and T. Pierog,
T. Pierog and K. Werner,
Nucl.Phys.Proc.Suppl. 196 (2009) 102-105
QGSJET01 (-m2) : N.N. Kalmykov, S. Ostapchenko, and A.I. Pavlov,
QGSJET01 (-m2) : N.N. Kalmykov, S. Ostapchenko, and A.I. Pavlov,
Nucl.Phys. B (Proc. Suppl.) 52B (1997) 17
QGSJETII-03 (-m11) : S. Ostapchenko,
......@@ -56,8 +56,8 @@ QGSJETII-03 (-m11) : S. Ostapchenko,
SIBYLL2.1 (-m6) : R. Engel, T.K. Gaisser, P. Lipari, and T. Stanev,
Proc. 26th Int. Cosmic Ray Conf., Salt Lake City (USA), 1 (1999) 415
and
E.-J. Ahn, R. Engel, T.K. Gaisser, P. Lipari, and T. Stanev,
and
E.-J. Ahn, R. Engel, T.K. Gaisser, P. Lipari, and T. Stanev,
Phys.Rev. D80 (2009) 094003
**********************************************************
......@@ -86,7 +86,7 @@ if you want to test some specific set of options.
**********************************************************
Troubleshoot
**********************************************************
Please first try in your build diretory to run:
make clean
rm CMakeCache.txt
......@@ -123,7 +123,7 @@ Options of CRMC:
-T [ --Test ] [=arg(=1)] Test mode
for projectile and target Id the following shortcuts are allowed :
for projectile and target Id the following shortcuts are allowed :
1 = PDG(2112) = proton
-1 = PDG(-2112) = anti-proton
12 = PDG(1000060120) = Carbon
......@@ -151,16 +151,16 @@ bin/crmc -T -m6
Options
**********************************************************
The details of the run can be controlled by the file "crmc.param".
The details of the run can be controlled by the file "crmc.param".
In this file the 3 first commented options are valid for EPOS 1.99 and EPOS LHC only.
- By defaults EPOS models use a simplified treatment of QGP in events where the energy density is high enough (including in pp). If you
- By defaults EPOS models use a simplified treatment of QGP in events where the energy density is high enough (including in pp). If you
uncomment the line "!switch fusion off" this will be disable and running time will be shorter (but data description will be worth since
the physics model will be incomplete). In that mode EPOS is comparable to PYTHIA model (no final state interaction).
the physics model will be incomplete). In that mode EPOS is comparable to PYTHIA model (no final state interaction).
- By default only the final particles are recorded in the output file like for other cosmic ray models.
Uncommenting "!set istmax 1" allows the user to have the full chain of mother/daughter from the beam to the final particles with EPOS models.
- By default only the final particles are recorded in the output file like for other cosmic ray models.
Uncommenting "!set istmax 1" allows the user to have the full chain of mother/daughter from the beam to the final particles with EPOS models.
The outfile is at least 2 times larger but includes the decayed particles and some special intermediate particles between the beam and
the real particles which allow the user to know where the particles were generated. The ids of such particles are the following :
......@@ -169,15 +169,15 @@ the real particles which allow the user to know where the particles were generat
92 : string. Mother of particles coming directly from string fragmentation of initial Pomerons (which do not participate to plasma formation).
93 : remnant. Mother of particles coming from the beam remnants.
Each primary particle has a mother with id 9x. For technical reasons all particles with the same id and same momentum are in fact the same
initial object (cluster, string or remnant) which was split in different local pieces to keep the correct production vertexes of the final
particles. If a set of particles with the same id (91 or 92) and the same momentum has 2 different mothers (one from the projectile and one
Each primary particle has a mother with id 9x. For technical reasons all particles with the same id and same momentum are in fact the same
initial object (cluster, string or remnant) which was split in different local pieces to keep the correct production vertexes of the final
particles. If a set of particles with the same id (91 or 92) and the same momentum has 2 different mothers (one from the projectile and one
from the target) it means that the string was produced by the interaction of this pair of nucleons or for a cluster that it was the closest
pair of nucleon which participate at the formation of the cluster (a cluster is formed by string pieces of different pairs but hepmc format
allows only one mother in that case).
- By default the cross-section is calculated by a numerical method with is valid only for h-p or h-A
(h being pion, kaon or nucleon) but not AB (nucleus-nucleus). If you uncomment "!set isigma 2" the inelastic
- By default the cross-section is calculated by a numerical method with is valid only for h-p or h-A
(h being pion, kaon or nucleon) but not AB (nucleus-nucleus). If you uncomment "!set isigma 2" the inelastic
cross-section will be always correct but it takes several minutes to compute it (see Note on simulating Heavy Ion events).
**********************************************************
......@@ -187,15 +187,15 @@ cross-section will be always correct but it takes several minutes to compute it
QGSJET01 and SIBYLL2.1 are limited to hA collisions with A<64 and h being a pion, a kaon or a nucleon.
Only EPOS (1.99 and LHC) and QGSJETII (03 and 04) models can run (h)AB collisions with A and B up to 208.
QGSJETII was never designed for HI collisions, so results should be interpreted with care (no final
state interactions).
EPOS 1.99 and LHC include a simplified treatment of final state interaction but doesn't include full hydro
state interactions).
EPOS 1.99 and LHC include a simplified treatment of final state interaction but doesn't include full hydro
simulations as already possible in EPOS 2. It gives a good overall description of HI data but it should not be
used to interpret QGP related observables (flow, jet quenching, etc ...) since the model is oversimplified in
used to interpret QGP related observables (flow, jet quenching, etc ...) since the model is oversimplified in
this distribution.
* For detailed simulations of HI collisions within EPOS 2, please contact K. Werner (werner@subatech.in2p3.fr).*
Concerning the cross-section calculation in EPOS, to avoid unnecessary long time calculation needed
for HI interactions, the default type of cross-section calculation is fast and good for hA but not reliable for
for HI interactions, the default type of cross-section calculation is fast and good for hA but not reliable for
AB collisions (OK for p-Pb, pi-C but not for C-C or Pb-Pb for instance). So if you really need to have
a correct cross-section for PbPb scattering for example, you should uncomment the line "set isigma 2" in the
crmc.param file but it will take several minutes to compute the cross-section each time you start CRMC.
......@@ -205,12 +205,22 @@ The HepMC (2.06) IO library has a default limitation of the number of produced p
Thus, events with more than that number of secondaries are either truncated or skipped. For heavy ion
collisions this easily can become a problem in crmc. Please consider changing thie limiation to a
higher value. We suggest 200000, which works for all models at LHC energies so far. We are happy for
feedback on this issue if problems are encountered.
The CRMC code automatically checks for the
limit in the HepMC library and uses this during compile time. In order to change the limit, you
need to re-compile HepMC package and change the value in the line
feedback on this issue if problems are encountered.
The CRMC code automatically checks for the
limit in the HepMC library and uses this during compile time. In order to change the limit, you
need to re-compile HepMC package and change the value in the line
#define HEPEVT_EntriesAllocation 199990
in the top of the file HepMC/HEPEVT_Wrapper.h
When you run into this limitation you will receive lots of warning messages.
When you run into this limitation you will receive lots of warning messages.
**********************************************************
Example for analysing hepmc files
******I****************************************************
In the directory ./ExampleAnalyser you can find code for a small
example that loops over the hepmc events and its particles. To
install make sure you have boost, hepmc, and root setup as for the
normal installation of crmc then run
"make" and "bin/analyse"
to run the code.
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment