gev380_eetunes_groot.py (raw text file)

# Validation script of =HepMC= : https://atlaswww.hep.anl.gov/hepsim/
# Run it inside hs-tools or Jas4pp framework  
# Task: Calculations of event shapes in e+e-
# S.Chekanov (ANL)

from java.awt import *
from java.lang import *
from javax.swing import JFrame
from proto import FileMC 
from org.jlab.groot.data import H1F
from org.jlab.groot.graphics import EmbeddedCanvas
from jhplot.utils import FileList
from hepsim import HepSim
from java.util import ArrayList
from hep.physics.vec import BasicHepLorentzVector 
from hep.physics.jet import EventShape
from hephysics.particle import LParticle 
import sys

evshape=EventShape()


NMax=2              # max number of files per tune. Set to -1 to run over all files

# default data location 
default_www="http://mc.hep.anl.gov/asc/hepsim/events/ee/380gev/pythia8_tunes_qedoff/"

# get data either from local file system or stream them from URL
url="" 
if len(sys.argv)>1:
   if sys.argv[1].startswith("http"):
         flist=HepSim.getList(sys.argv[1])
         url=sys.argv[1]
   else: flist=FileList.get(sys.argv[1],"promc")
   if (len(sys.argv)>2): TotalEvents=int(sys.argv[2])
else:
   url=default_www; flist=HepSim.getList(url)
if len(flist)==0: print "Error: No input file!"; sys.exit(0)
else: print "Reading "+str(len(flist))+" files"

# extract files with in certain tune ranges
xmin=HepSim.getRanges(flist,"_tune","_");

# reduce range
# del xmin[-1]; del xmin[-1]; 
print "Tune ranges: ",xmin

# create file list using mass ranges 
filelist=[]
histo1=[] # oblateness: Major Thrust - Minor Thrust 
histo2=[] # Major Thrust
histo3=[] # Minor Thrust 
 

n=0
for i in xmin:
     ma=[]
     for j in flist:
        if (j.find("tune"+str(i)+"_")>-1): ma.append(j)
     filelist.append(ma)
     hh=H1F("tune="+str(i),40,0,0.4)   
     hh.setTitleX("O");
     hh.setTitleY("d sigma / d O [pb]");
     hh.setLineColor(n) # set color map 
     histo1.append( hh ) 
     
     hh=H1F("tune="+str(i),40,0,0.4)
     hh.setLineColor(n) #  color map 
     hh.setTitleX("T_{Major}");
     hh.setTitleY("d sigma / d T_{Major}  [pb]");
     histo2.append( hh ) 

     hh=H1F("tune="+str(i),40,0,0.4)
     hh.setLineColor(n) #  color map 
     hh.setTitleX("T_{Minor}");
     hh.setTitleY("d sigma / d T_{Minor}  [pb]");
     histo3.append( hh )       

     n=n+1


print "Total number of e+e- tunes =",len(filelist)
Total=len(filelist)
for ipt in range(Total):
  ptlist=filelist[ipt]
  print ipt, "Tune = "+str(xmin[ipt]) 
  cross=0; xcross=0; nev=0;  Nfiles=len(ptlist)
  if (NMax>-1):  Nfiles=NMax           # limit to NMax files 
  for m in range(Nfiles):              # loop over all files in this list    
    file=FileMC(url+ptlist[m])         # get input file
    print "Read=",ptlist[m]
    header = file.getHeader()
    un=float(header.getMomentumUnit()) # conversion units
    for i in range(file.size()):
      nev=nev+1
      if (nev%200==0):
           if (m==0): print "Event=",nev
           else: print "Event=",nev," done=",int((100.0*m)/Nfiles),"%"
      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
           p=BasicHepLorentzVector(e,px,py,pz)
           particles.add(p)          # add particle to the list 
      # see: http://java.freehep.org/lib/freehep/api/hep/physics/jet/EventShape.html
      evshape.setEvent(particles)
      oblateness=evshape.oblateness()
      majoraxis=evshape.majorAxis()
      minoraxis=evshape.minorAxis()
      trustaxis=evshape.thrustAxis()
      trust=evshape.thrust()
      histo1[ipt].fill(oblateness)
      histo2[ipt].fill(trust.y()) # Major 
      histo3[ipt].fill(trust.z()) # Minor 

    stat = file.getStat()
    cross=cross+stat.getCrossSectionAccumulated() 
    erro=stat.getCrossSectionErrorAccumulated()
    file.close()
  cross=cross/Nfiles
  lumi=nev/cross
  binwidth=histo1[ipt].getAxis().getBinWidth(1)
  histo1[ipt].divide(lumi*binwidth)   # for cross section 
  histo2[ipt].divide(lumi*binwidth)   # for cross section 
  histo3[ipt].divide(lumi*binwidth)   # for cross section 

  # print ipt, ") Cross section = %.3e pb"%cross, " Lumi=%.3e pb"%lumi

crossD=histo1[0].integral() # get total cross section 


c1 = EmbeddedCanvas()
c1.divide(1,3)

c1.cd(0)
c1.getPad(0).getAxisX().setRange(0, 0.3)
c1.getPad(0).getAxisY().setRange(0, 200)
c1.draw(histo1[0])
for i in range(1,len(histo1)): c1.draw(histo1[i],"same")


c1.cd(1)
for h in histo2: c1.draw(h)
c1.getPad(1).getAxisX().setRange(0, 0.3)
c1.getPad(1).getAxisY().setRange(0, 500)
c1.draw(histo2[0])
for i in range(1,len(histo2)): c1.draw(histo2[i],"same")


c1.cd(2)
c1.getPad(2).getAxisX().setRange(0, 0.3)
c1.getPad(2).getAxisY().setRange(0, 500)
c1.draw(histo3[0])
for i in range(1,len(histo3)): c1.draw(histo3[i],"same")


frame = JFrame("Event shapes")
frame.setSize(600,700)
frame.add(c1)
frame.setVisible(True)

name="output"
if len(sys.argv[0])>0: name=sys.argv[0].replace(".py","")
c1.save(name+".png");    print name+".png created"
print "All done"