====== FCC-ee tutorial using HepSim files ======

Here is mini-tutorial explaining how to use Monte Carlo truth-level files from the HepSim for the FCC-ee studies. The setup of the computer environment is described in [[fcs:fccee:setup| this section]]. The HepSim part of this tutorial does not require CVMFS and can be done on any computer with java installed.


====== FCC-ee tutorial ======

See [[https://hep-fcc.github.io/fcc-tutorials/]] The source files are on GitHub at [[https://github.com/HEP-FCC/fcc-tutorials]]. Check this for making contributions: [[https://hep-fcc.github.io/fcc-tutorials/CONTRIBUTING.html]]



===== Using files from HepSim  =====

This part of the tutorial does not use any C++ specific libraries and can be done on any computers with Java installed.  Check java:
<code bash>
java -version
</code>

Typically, it will tell openjdk version "1.8.0_352" or higher Java version.

<code bash>
bash   #  set bash if you haven't done this before 
wget https://atlaswww.hep.anl.gov/hepsim/soft/hs-toolkit.tgz -O - | tar -xz;
source hs-toolkit/setup.sh
</code>

Let's look at a few events: Z-> Z H, where Z-> nunu, and H decays to  bbar. The CM energy is 250 GeV. The sample is described in https://atlaswww.hep.anl.gov/hepsim/info.php?item=353 

First, print all files with Higgs processes:

<code bash>
hs-find higgs
</code>

Then grab the file with H to bbar at e+e-:

<code bash>
hs-ls gev250ee_pythia8_zhiggs_nunubbar
</code>

Let's download 10 files (in 2 threads):

<code bash>
hs-get gev250ee_pythia8_zhiggs_nunubbar data 2 10
</code>

We should have 10 files in the directory "data". Let's take a look at a single file. We what to check how many events in the file, how many events etc:

<code bash>
hs-info  data/gev250ee_zh_nunubb_001.promc
</code>

Do you want to print 1st event? Do this:

<code bash>
hs-info data/gev250ee_zh_nunubb_001.promc 1
</code>

Want  to examine the log file? Do this:

<code bash>
hs-log  data/gev250ee_zh_nunubb_001.promc
</code>

Let's study each event in the GUI mode (needs X-session!). Start this GUI and click each event number using the left panel.

<code bash>
hs-view  data/gev250ee_zh_nunubb_001.promc
</code>

=====  Validation =====

Now we know a lot about this files. Let's make a simple code to run over 10,000 events (2 files) and draw some useful histograms.
Click the tab below to see the code. This code may look rather long, but it does a lot: It reads the truth-level events, reconstructs jets, and makes plots of the pT (differential cross section), Rapidity distribution and the invariant mass of the 2 jets. It also saves the image in SVG file  and an XML file for examination:

<hidden Click here to show the Python/Jython code>
<code>
# Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
# S.Chekanov (ANL)
from java.lang import *
from proto import FileMC   # import FileMC
from java.util import ArrayList
from hep.physics.vec import BasicHepLorentzVector
from jhplot.utils import FileList
from java.awt import Color,Font
from jhplot import  HPlot,P1D,HLabel,H1D
from jhplot.io import HBook
from hephysics.particle import LParticle
from hep.physics.jet import DurhamJetFinder,FixNumberOfJetsFinder
from hepsim import HepSim
import sys

Nfiles=2 # 2 files to process for checking 
flist=FileList.get("data","promc")
print "Found "+str(len(flist))+" files. Nr files to process= ",Nfiles

fjet=DurhamJetFinder(0.05) # DurhamJet with y_cut=0.05 
#fjet=FixNumberOfJetsFinder(2)

h1= H1D("PT_{jet}",50,10,200.)
#h1.doc()               # check its methods
h2= H1D("y_{jet}",50,-5,5.)
h3= H1D("Mass",100,100,150)

cross=0; nev=0;  # Nfiles=len(flist)
for m in range(Nfiles):             # loop over all files in this list    
   file=FileMC(flist[m])        # get input file
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # conversion units
   lunit=float(header.getLengthUnit())
   if m==0:
       print "ProMC v=",file.getVersion(), "E(varint)=",un,"L(varint)t=",lunit
   for i in range(file.size()):
      nev=nev+1
      if (nev%500==0): print "Event=",nev
      eve = file.read(i)
      pa = eve.getParticles()   # particle information
      ve = eve.getEvent()       # event information
      particles=ArrayList() # list of particles 
      for j in range(pa.getPxCount()):
        if (pa.getStatus(j)==1):     # stable
           px=pa.getPx(j)/un         # pX
           py=pa.getPy(j)/un         # pY
           pz=pa.getPz(j)/un         # pZ
           e= pa.getEnergy(j)/un        # energy
           mass= pa.getMass(j)/un    # mass
           apdg=abs(pa.getPdgId(j))
           if (apdg==12 or apdg==14 or apdg==16): continue # no neutrino
           p=BasicHepLorentzVector(e,px,py,pz)
           particles.add(p) # add particle to the list 
      fjet.setEvent(particles)
      jets=[] # make a new list with jets 
      for j in range(fjet.njets()):
         # print "Jet=",i," Nr of particles in jet =",fjet.nParticlesPerJet(i) 
         pjet=fjet.jet(j)   # make  HepLorentzVector
         ee=pjet.t()  # magnitude
         p3=pjet.v3() # 3-vector  
         pl=LParticle(p3.x(),p3.y(),p3.z(),ee) # convert to a class with convinient kinematics 
         h1.fill(pl.perp())
         h2.fill(pl.rapidity())
         jets.append(pl)
      if (len(jets)>1):
              jets[0].add(jets[1])
              h3.fill(jets[0].calcMass() )

   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

lumi=nev/cross;
print "Lumi=%.3e pb"%lumi
print "Total cross section (pb)=",cross
c1 = HPlot("=HepSim=",600,800,1,3)
c1.visible(True)
# c1.setAutoRange()
# c1.setLogScale(1,True)
c1.setMarginLeft(100)
c1.setNameX("p_{T} [GeV]")
c1.setNameY("d &sigma; / d p_{T} [pb/GeV]")
c1.setRange(0,200,0.0, 0.01)

Z=h1.getDividedByBinWidth()
Z.scale(1.0/lumi)
xsec1=P1D(Z)
xsec1.setErr(1)       # show errors
xsec1.setColor(Color.blue)
c1.draw(Z)
# h1.toTable() 

l1=HLabel("Durham jets  &sigma;=%.3e pb"%cross, 0.4, 0.75, "NDC")
l1.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l1)
c1.cd(1,2) # plot rapidity
c1.setNameX("Rapidity")
c1.setNameY("Events")
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(-5,5)
c1.draw(h2)

# plot  invariant mass 
c1.cd(1,3)
c1.setNameX("M_{jj}")
c1.setNameY("Events")
c1.setMarginLeft(100)
c1.setAutoRange()
c1.setRangeX(100,150)
c1.draw(h3)

# create file/image using name of the file
name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write(xsec1)
file.close()
c1.export(name+".svg");    print name+".svg created"
# xsec.toFile(name+".txt");  print name+".txt created"
# sys.exit(0)
</code>
</hidden>

Copy the code shown below to a file "validatation.py" and run it as 

<code>
hs-run validatation.py
</code>

It will create the image with jets and a Higgs peak for 2-jet mass.

{{:fcs:validation_fccee.png?nolink&400|}}


=====  Data streaming =====

One can also analyze truth-level files using data streaming, i.e. without downloading the actual input files. 
This example shows data for Higgs to bbar at 240 GeV e+e-.
Click the code shown below, save it as "example.py" and run it inside Jas4pp (https://atlaswww.hep.anl.gov/asc/jas4pp/) or using hs-tools package (as "hs-run example.py").
 
<hidden Click here to show the Python/Jython code>
<code python>
# illustration of how to make pT distribution in real time using data streaming
# S.Chekanov (ANL)
from java.awt import Color,Font
from java.lang import *
from proto import FileMC            
from jhplot import HPlot,H1D,H2D,P1D,HLabel  # graphics
from jhplot.utils import FileList
from hephysics.particle import LParticle 
from hephysics.hepsim import PromcUtil
from hephysics.jet import JetN2 
from hepsim import HepSim
import sys,math 

TotalEvents=10000

dataset="gev240ee_pythia8_higgs_bbar"
tag="" # only truth level
url=""
NMax=0
if len(sys.argv)>1:
   flist=FileList.get(sys.argv[1],"promc")
   if (len(sys.argv)>2): NMax=int(sys.argv[2])
else:
   sites=HepSim.urlRedirector(dataset)
   url=sites[0]+"/"+tag
   flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files. Nr of files= ",NMax

#h1.doc()               # check its methods
h1= H1D("pT [GeV]",100,0,100)

c1 = HPlot("=HepSim=",500,500)
c1.visible(True)
c1.setAutoRange(True)
c1.setMarginLeft(90)
c1.setLegend(0)

cross=0; nev=0;  Nfiles=len(flist)
for m in range(Nfiles):              # loop over all files in this list    
   file=FileMC(url+flist[m])         # get input file
   print "Read=",flist[m]
   header = file.getHeader()
   un=float(header.getMomentumUnit()) # conversion units
   lunit=float(header.getLengthUnit())
   if m==0:
       print "ProMC v=",file.getVersion(), "M unit=",un,"L unit=",lunit 
   if (nev>TotalEvents): print "Max Nr of events are done"; break # stop after some events
   for i in range(file.size()):
      if (nev>TotalEvents): break
      nev=nev+1
      if (nev%1000==0):
           if (Nfiles==1): print "Event=",nev
           else: 
                print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
                c1.clearData()
                #c1.setNameX("p_{T}(lead. jet) [GeV]")
                #c1.setNameY("Entries")
                c1.draw(h1)
                print "Update canvas"
      eve = file.read(i)
      pa = eve.getParticles()    # particle information
      nt=0; xsum=0; 
      darkpi={}
      #print " "
      for j in range(pa.getPxCount()):
         apdg=abs(pa.getPdgId(j))
         par=LParticle(pa.getPx(j)/un,pa.getPy(j)/un,pa.getPz(j)/un,pa.getEnergy(j)/un,pa.getMass(j)/un)
         nt=nt+1
         h1.fill(par.perp()) # calculate transverse mementa 
   stat = file.getStat()
   cross=stat.getCrossSectionAccumulated()
   erro=stat.getCrossSectionErrorAccumulated();
   file.close()

# show final cross section
c1.clearData()
c1.clear()
#c1.setLogScale(1,1)
#c1.setRangeY(0,01)
c1.setRangeX(0,100)
c1.setMarginLeft(90)
c1.setNameX("p_{T} [GeV]")
c1.setNameY("Events")

c1.draw(h1)
print h1.getStat()
#l1=HLabel("&sigma;=%.3e pb"%cross, 0.4, 0.75, "NDC")
#l1=HLabel("<m>=%.3e GeV"%h2.mean(), 0.4, 0.75, "NDC")
#l1.setFont(Font("Helvetica", Font.PLAIN, 15))
#c1.add(l1)

l2=HLabel("=HepSim=",0.6,0.85, "NDC")
l2.setColor(Color.gray)
l2.setFont(Font("Helvetica", Font.PLAIN, 14))
c1.add(l2)

# create file/image using name of the file
name="output"
from jhplot.io import HBook 
file=HBook(name+".jdat","w"); print name+".jdat created"
file.write("pT",h1)
file.close()
c1.export(name+".svg");    print name+".svg created"
#xsec.toFile(name+".txt");  print name+".txt created"
#sys.exit(0)
</code>
</hidden>

If you run this code as Jython script inside Jas4pp, you will see a histogram updated in real time:

{{ :fcs:fccee:screenshot_from_2023-05-24_08-57-19.png?400 |}}


=====  Analysis of LCIO files =====

Full example of how to analyse LCIO files is given in: https://atlaswww.hep.anl.gov/asc/wikidoc/doku.php?id=asc:jas4pp#reading_lcio_files

===== Preparation for fast simulation =====

For fast simulations we need to install the C++ ProMC library. After it is installed, the compiled Delphes program can start using the ProMC files natively. 
Let's take this package from https://atlaswww.hep.anl.gov/asc/promc/ and install it.  Make sure you have the gcc compiler , zip library (run "yum install python-pip python-setuptools zlib1g-dev zlib-devel" if you do not have them) and CERN ROOT. If you are using CERN or OSG login nodes, you can get these packages as:


Use the  key4hep setup with gcc11. It will also link the ROOT libraries:

<code bash>
source /cvmfs/sw.hsf.org/key4hep/setup.sh 
</code>

If ( /cvmfs/sw.hsf.org/key4hep/ ) is not available on your system, use the older gcc8 version:

<code bash>
source /cvmfs/sft.cern.ch/lcg/views/LCG_97/x86_64-centos7-gcc8-opt/setup.sh
</code>

Then install ProMC (currently it is an external package):

<code bash>
wget http://atlaswww.hep.anl.gov/asc/promc/download/current.php -O ProMC.tgz
tar -zvxf ProMC.tgz
cd ProMC
./build.sh                      # build all source files
./install.sh lib               # install into the "lib" directory
source lib/promc/setup.sh  # make it available
</code>

After the installation, make sure you see a bunch of executable files:

<code bash>
echo $PROMC/bin
</code>

This package will be used for fast simulation. In addition, you can perform various transformation of the event records. For example, move your data to ROOT:


<code bash>
cd $PROMC/examples/promc2root
make                          # creates converter promc2root
cd ../../../../../data/   # navigate to the directory with data
$PROMC/examples/promc2root/promc2root gev250ee_zh_nunubb_001.promc  gev250ee_zh_nunubb_001_truth.root
</code>
Note some system picks the native "zlib library". If you see an error in compilation, remove "-lz" string on the line  23 of Makefile.

The last command will create ROOT files with the tree "ProMC":

{{ :fcs:fccee:screenshot_from_2023-03-06_10-31-06.png?nolink&400 |}}

Similarly, one can do transformation oto the popular text formats:  See 

<code bash>
$PROMC/examples/promc2lhef
$PROMC/examples/promc2hepmc 
</code>

The same directory "examples" contains the reverse transformations.


===== Creating Pythia8 events =====

In this example, instead of downloading files with e+e- events, we will generate Pythia8 events at the truth level. We create 100k events at 250 GeV CM energy with all Higgs processes. Then we will force H to decay to 2 b jets.

(1) Step 1. Get the package:

<code bash>
wget https://atlaswww.hep.anl.gov/hepsim/examples/package_ee.tgz
tar -zvxf package_ee.tgz
</code>

(2) Compile the program. The steering card is located in "cards" that defines how many events and what are the processes.

<code bash>
cd package_ee
make
</code>

(3) Run the Pythioa8 events using the supplied bash script.   This creates 10 output files in the directory "out". You can view the output  and the log file as 

<code bash>
./A_RUN_gev250_higgs_bbar
</code>

This script creates the output file in the directory "out". One can view events in this file as:

<code bash>
cd out
promc_browser gev250ee_higgs_bbar_001.promc
</code>
===== Fast simulation =====

To perform fast simulation, install Delphes using the standard steps:

<code bash>
wget http://cp3.irmp.ucl.ac.be/downloads/Delphes-3.5.0.tar.gz
tar -zxf Delphes-3.5.0.tar.gz
cd Delphes-3.5.0
make
</code>
On some systems, the compilation may fail if the ZLIB  library is already present, i.e. the error will be "inflateValidate@ZLIB". In this case, remove "-lz" on the line 45 of the Makefile. 

The last line of the compilation should have "DelphesProMC" (binary file to process events).

<code bash>
cd ..
./Delphes-3.5.0/DelphesProMC  ./Delphes-3.5.0/cards/delphes_card_CircularEE.tcl \
data/gev250ee_zh_nunubb_001.root  \
data/gev250ee_zh_nunubb_001.promc
</code>

The output file will be  data/gev250ee_zh_nunubb_001.root 

===== Full simulation and file conversions =====

For full simulations, we will need to convert ProMC files to the LCIO files.  You can use the converter posted to https://github.com/hepsim-org/promc2other
Read the README of this simple tool.

The C++ version of the ProMC package can also convert files to HEPMC2 and LHEF.  Look at the directories:

<code>
cd $PROMC/promc2hepmc   # tool to convert to HEPMC2 
cd $PROMC/promc2lhe         #  tool to convert to LHEf
</code>

The key4hep requires HEPMC3 format (not HEPMC2). Therefore, the current option is to convert PROMC to HEMC2 and then convert HEPMC3. Alternatively, key4hep can take LHE files.

First, get the script:

<code>
wget http://fccsw.web.cern.ch/utils/converters/lhe2edm.sh
</code>
 
Then use it like this:

<code>
$ source /cvmfs/sw.hsf.org/key4hep/setup.sh
$ ./lhe2edm.sh inputfile.lhe output.edm [--nevt=1000] [--stable]
</code>
 

You can specify the number of events in the file in the case the guess based on grepping <event> is wrong.

If you specify ‘—stable’ only the stable particles (status == 1) are kept in the final file.



===== Some outreach  =====

If you are a student who has  never seen a typical e+e- event in real life, you may use this simple tool to display and view e+e- events. Here we will explore  typical e+e- events at the CM energy 250 GeV with Higgs  decaying to two Z bosons. Then Z-bosons decay to 4 leptons.  The events will be visualized in the SID detector from the ILC project (events for FCC-ee will not look very different, but the FCC-ee detector will have somewhat different geometry).  The program uses Jas3/Wired display used for the ILC R&D. 

Unlike the tutorial discussed above, this example only requires a personal Linux/Windows/Mac computer with Java installed. It will use the [[https://atlaswww.hep.anl.gov/asc/jas4pp/|Jas4pp]] program and files from HepSim.

<code bash>
wget https://atlaswww.hep.anl.gov/asc/jas4pp/download/jas4pp-1.7.tgz
tar -zvxf jas4pp-1.7.tgz
cd jas4pp
wget https://mc.hep.anl.gov/asc/hepsim/events/ee/250gev/pythia6_higgs_zz_4l/rfull001/pythia6_higgs_zz_4l_0001_pandora.slcio
jaspp pythia6_higgs_zz_4l_0001_pandora.slcio
</code>

You will see an editor. Use the left panel and expand "DataSets". Then go to "File->New->Wired 4 View. You will see a black screen. The press the small blue arrow from the top menu. It is called "Go".  

You can also display very  busy events with H-> bbar if you will download files from [[https://atlaswww.hep.anl.gov/hepsim/show.php?item=151&reco=rfull001 |151rfull001 link]].


{{ :fcs:fccee:250gev.png?nolink&500 |}}

You can explore various detector components, data collections etc. Use File->New->LCSim event browse. etc. You can also display other collision processes and different detectors using the HepSim repo.

 --- //[[hepsim@anl.gov|Sergei Chekanov]] 2023/04/03 21:08//