Differences

This shows you the differences between two versions of the page.

@@ Line 1: / Line 1: @@
+====== 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.
+ --- //[[[email protected]|Sergei Chekanov]] 2023/04/03 21:08//