Differences

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

--- fcs:fccee:tutorial [2023/04/18 20:11]
hepsim17 [Full simulation and file conversions]
+++ fcs:fccee:tutorial [2023/05/24 13:25]
hepsim17 [Data streaming]
@@ Line 74: / Line 74: @@
 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/
 # S.Chekanov (ANL)
@@ Line 210: / Line 210: @@
 {{: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 |}}
 ===== Preparation for fast simulation =====