Commit cb91bc4b authored by David Blyth's avatar David Blyth

Added jpsi_dvmp_delta_pt

parent 31cce795
......@@ -2,3 +2,4 @@
*.proio
*.pdf
*.png
go.sum
This repo is currently a mess. Please do not use until it is cleaned up.
module github.com/decibelcooper/eicplot
require (
github.com/pkg/profile v1.2.1
github.com/proio-org/go-proio v0.2.1
github.com/proio-org/go-proio-pb v0.0.0-20180827221749-ab81f0556d2d
go-hep.org/x/hep v0.14.0
go-hep.org/x/hep v0.15.0
gonum.org/v1/plot v0.0.0-20180905080458-5f3c436ce602
)
package main
import (
"flag"
"fmt"
"image/color"
"log"
"math"
"os"
"github.com/pkg/profile"
"github.com/proio-org/go-proio"
"github.com/proio-org/go-proio-pb/model/eic"
"go-hep.org/x/hep/hbook"
"go-hep.org/x/hep/hplot"
"gonum.org/v1/plot"
"gonum.org/v1/plot/vg"
"github.com/decibelcooper/eicplot"
)
func printUsage() {
fmt.Fprintf(os.Stderr, `Usage: `+os.Args[0]+` [options] <proio-input-file>
options:
`,
)
flag.PrintDefaults()
}
func main() {
defer profile.Start().Stop()
var (
title = flag.String("title", "", "plot title")
output = flag.String("output", "out.png", "output file")
)
flag.Usage = printUsage
flag.Parse()
if flag.NArg() < 1 {
printUsage()
log.Fatal("Invalid arguments")
}
p, _ := plot.New()
p.Title.Text = *title
p.X.Label.Text = "Transverse Momentum Transfer (GeV)"
p.X.Tick.Marker = eicplot.PreciseTicks{NSuggestedTicks: 5}
p.Y.Tick.Marker = eicplot.PreciseTicks{NSuggestedTicks: 5}
p.Y.Scale = eicplot.LogScale{}
var hists []*hbook.H1D
for _, filename := range flag.Args() {
hists = append(hists, makeHists(filename)...)
}
for i, hist := range hists {
lineColor := color.RGBA{A: 255}
switch i {
case 1:
lineColor = color.RGBA{G: 255, A: 255}
case 2:
lineColor = color.RGBA{B: 255, A: 255}
case 3:
lineColor = color.RGBA{R: 255, B: 127, G: 127, A: 255}
}
h := hplot.NewH1D(hist)
h.FillColor = nil
h.LineStyle.Color = lineColor
h.Infos.Style = hplot.HInfoNone
p.Add(h)
}
p.Save(6*vg.Inch, 4*vg.Inch, *output)
}
const jpsiMass = 3.096916
func makeHists(filename string) []*hbook.H1D {
deltaPTHist := hbook.NewH1D(50, 0, 1)
deltaPTTruthHist := hbook.NewH1D(50, 0, 1)
reader, err := proio.Open(filename)
if err != nil {
log.Fatal(err)
}
for event := range reader.ScanEvents() {
ids := event.TaggedEntries("Reconstructed")
tracks := []*eic.Track{}
for _, id := range ids {
track, ok := event.GetEntry(id).(*eic.Track)
if ok && len(track.Segment) > 0 {
tracks = append(tracks, track)
}
}
if len(tracks) == 3 {
var p [2]float64
for _, track := range tracks {
p[0] += *track.Segment[0].Poq.X
p[1] += *track.Segment[0].Poq.Y
}
deltaPT := math.Sqrt(math.Pow(p[0], 2) + math.Pow(p[1], 2))
deltaPTHist.Fill(deltaPT, 1)
}
ids = event.TaggedEntries("GenStable")
protons := []*eic.Particle{}
for _, id := range ids {
part, ok := event.GetEntry(id).(*eic.Particle)
if ok && *part.Pdg == 2212 {
protons = append(protons, part)
}
}
if len(protons) == 1 {
deltaPT := math.Sqrt(math.Pow(float64(*protons[0].P.X), 2) + math.Pow(float64(*protons[0].P.Y), 2))
deltaPTTruthHist.Fill(deltaPT, 1)
}
}
reader.Close()
return []*hbook.H1D{deltaPTTruthHist, deltaPTHist}
}
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment