Differences

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

Link to this comparison view

Both sides previous revision Previous revision
Next revision 		Previous revision
fcs:fccee:tutorial [2023/04/10 14:35]
hepsim17 [Using files from HepSim] 		fcs:fccee:tutorial [2023/05/24 14:06] (current)
hepsim17 [Data streaming]
Line 69: Line 69:
 </code> </code>
  
-======  Validation ====== +=====  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. 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: 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> +<hidden Click here to show the Python/Jython code> 
-<code python>+<code>
 # Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/ # Part of =HepMC= : https://atlaswww.hep.anl.gov/asc/hepsim/
 # S.Chekanov (ANL) # S.Chekanov (ANL)
Line 105: Line 105:
 cross=0; nev=0;  # Nfiles=len(flist) cross=0; nev=0;  # Nfiles=len(flist)
 for m in range(Nfiles):             # loop over all files in this list     for m in range(Nfiles):             # loop over all files in this list    
-   file=FileMC(url+flist[m])        # get input file+   file=FileMC(flist[m])        # get input file
    header = file.getHeader()    header = file.getHeader()
    un=float(header.getMomentumUnit()) # conversion units    un=float(header.getMomentumUnit()) # conversion units
Line 123: Line 123:
            py=pa.getPy(j)/un         # pY            py=pa.getPy(j)/un         # pY
            pz=pa.getPz(j)/un         # pZ            pz=pa.getPz(j)/un         # pZ
-           e= pa.getEnergy(j)/un     # energy+           e= pa.getEnergy(j)/un        # energy
            mass= pa.getMass(j)/un    # mass            mass= pa.getMass(j)/un    # mass
            apdg=abs(pa.getPdgId(j))            apdg=abs(pa.getPdgId(j))
Line 169: Line 169:
 # h1.toTable()  # h1.toTable() 
  
-l1=HLabel("Durham jets  &sigma;=%.3e pb"%cross, 0.3, 0.75, "NDC")+l1=HLabel("Durham jets  &sigma;=%.3e pb"%cross, 0.4, 0.75, "NDC")
 l1.setFont(Font("Helvetica", Font.PLAIN, 14)) l1.setFont(Font("Helvetica", Font.PLAIN, 14))
 c1.add(l1) c1.add(l1)
-l2=HLabel("=HepSim=",0.7,0.87, "NDC") +c1.cd(1,2) # plot rapidity
-l2.setColor(Color.gray) +
-l2.setFont(Font("Helvetica", Font.PLAIN, 14)) +
-c1.add(l2) +
-# plot rapidity +
-c1.cd(1,2)+
 c1.setNameX("Rapidity") c1.setNameX("Rapidity")
 c1.setNameY("Events") c1.setNameY("Events")
Line 215: Line 210:
  
 {{:fcs:validation_fccee.png?nolink&400|}} {{: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 ===== ===== Preparation for fast simulation =====
Line 300: Line 422:
 </code> </code>
  
-This script creates the output file in the directory "out". You ca view events in this file as:+This script creates the output file in the directory "out". One can view events in this file as:
  
 <code bash> <code bash>
Line 328: Line 450:
  
 The output file will be  data/gev250ee_zh_nunubb_001.root  The output file will be  data/gev250ee_zh_nunubb_001.root 
-===== Full simulation =====+ 
 +===== 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 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. 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  ===== ===== 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). Unlike the tutorial discussed above, this example 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.+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> <code bash>
fcs/fccee/tutorial.1681137358.txt.gz · Last modified: 2023/04/10 14:35 by hepsim17
Except where otherwise noted, content on this wiki is licensed under the following license: CC Attribution-Share Alike 4.0 International


  • Show page
  • Old revisions
  • Backlinks
  • Recent Changes
  • Media Manager
  • Sitemap