simpleFluxPlot.py

Go to the documentation of this file.

import ROOT, sys, os



fluxpaths = ["/lbne/data/users/kirby/g4lbne_output/test2_300k_B_011.root","/lbne/data/users/ljf26/fluxfiles/g4lbne/v3r2p4/QGSP_BERT/Nominal/200kA/flux/"]
fluxfiles = []
fluxchains = []

for i in range(0,len(fluxpaths)):
    fluxfiles.append([])
    if os.path.isdir(fluxpaths[i]):
        print "The path you specified ("+fluxpaths[i]+") appears to be a directory... looking for flux files in that directory."
        
        for root,dirs,files in os.walk(fluxpaths[i]):
            for file in files:
                temppath = os.path.join(root,file)

                if file.startswith("g4lbne_") and file.endswith(".root"):
                    fluxfiles[i].append(temppath)
        print "Found ",len(fluxfiles[i])," flux files to consider"
    else:
        print "Drawing fluxes from file:",fluxpaths[i]
        fluxfiles[i].append(fluxpaths[i])


##### CREATE CHAINS OF FLUX NTUPLES #####            

    fluxchains.append(ROOT.TChain("nudata"))
    for file in fluxfiles[i]:
        fluxchains[i].Add(file)

vars_to_draw1 = ["beamX","beamY","protonX","protonY","protonPx","protonPy","protonPz","nuTarZ","hornCurrent","Ndxdz","Ndydz","Npz","Nenergy","NdxdzNear[0]","NdydzNear[0]","NenergyN[0]","NWtNear[0]","NdxdzFar[0]","NdydzFar[0]","NenergyF[0]","NWtFar[0]","Norig","Ndecay","Ntype","Vx","Vy","Vz","ppdxdz","ppdydz","pppz","ppenergy","ptype","ppvx","ppvy","ppvz","muparpx","muparpy","muparpz","Necm","Nimpwt","tvx","tvy","tvz","tpx","tpy","tpz","tptype","tgen"]
vars_to_draw2 = ["beamX","beamY","protonX","protonY","protonPx","protonPy","protonPz","nuTarZ","hornCurrent","Ndxdz","Ndydz","Npz","Nenergy","NdxdzNear[0]","NdydzNear[0]","NenergyN[0]","NWtNear[0]","NdxdzFar[0]","NdydzFar[0]","NenergyF[0]","NWtFar[0]","Norig","Ndecay","Ntype","Vx","Vy","Vz","ppdxdz","ppdydz","pppz","ppenergy","ptype","ppvx","ppvy","ppvz","muparpx","muparpy","muparpz","Necm","Nimpwt","tvx","tvy","tvz","tpx","tpy","tpz","tptype","tgen"]


canv = ROOT.TCanvas("canvas","canvas")
hists1 = []
hists2 = []
for i in range(0,len(vars_to_draw1)):
    fluxchains[0].Draw(vars_to_draw1[i]+">>"+vars_to_draw1[i]+"_1");
    hists1.append(ROOT.gDirectory.Get(vars_to_draw1[i]+"_1"))
    hists2.append(hists1[i].Clone(vars_to_draw2[i]+"_2"))
    fluxchains[1].Draw(vars_to_draw1[i]+">>"+vars_to_draw2[i]+"_2");
    
    hists2[i].Scale(3.0/190.0);

    hists1[i].SetLineColor(4)
    hists1[i].SetLineWidth(2)
    hists1[i].Draw()
    hists2[i].SetLineColor(2);
    hists2[i].Draw("SAME")
    
    canv.Print("plots/"+vars_to_draw1[i]+".eps")
    canv.Print("plots/"+vars_to_draw2[i]+".png")

os.system("mkdir ~/public_html/temp/ntuple_comps/")
os.system("cp plots/*.png ~/public_html/temp/ntuple_comps/")
os.system("cp dirindex.sh ~/public_html/temp/ntuple_comps/")
os.system("chmod u+x ~/public_html/temp/ntuple_comps/dirindex.sh")
os.system("cd plots/; ~/public_html/temp/ntuple_comps/dirindex.sh | tee ~/public_html/temp/ntuple_comps/index.html")

"""
canv = ROOT.TCanvas("NeutrinoTheta","NeutrinoTheta")
canv.cd()
H_NeutrinoTheta = ROOT.TH1F("H_NeutrinoTheta","H_NeutrinoTheta",100,0,0.5);
FluxChain.Draw("360.0/(2*TMath::Pi())*TMath::ACos(TMath::Sqrt(1/(1+pow(NdxdzNear[0],2)+pow(NdydzNear[0],2))))>>H_NeutrinoTheta","Nimpwt*NWtNear[0]")
#FluxChain.Draw("Nimpwt");

canv.Print("NeutrinoTheta.png");
canv..Print("NeutrinoTheta.eps");


H_NeutrinoDxdz = ROOT.TH1F("H_NeutrinoDxdz","H_NeutrinoDxdz",200,-0.03,0.03);
FluxChain.Draw("NdxdzNear[0]>>H_NeutrinoDxdz","Nimpwt*NWtNear[0]")

canv.Print("NeutrinoDxdz.png");
canv.Print("NeutrinoDxdz.eps");


H_NeutrinoDydz = ROOT.TH1F("H_NeutrinoDyDz","H_NeutrinoDyDz",200,-0.03,0.03);
FluxChain.Draw("NdydzNear[0]>>H_NeutrinoDxdz","Nimpwt*NWtNear[0]")

canv.Print("NeutrinoDydz.png");
canv.Print("NeutrinoDydz.eps");

canv.cd()
H_NeutrinoEnergy = ROOT.TH1F("H_NeutrinoEnergy","H_NeutrinoEnergy",100,0,20);
FluxChain.Draw("NenergyN[0]>>H_NeutrinoEnergy","Nimpwt*NWtNear[0]")
H_NeutrinoEnergy.Draw();

canv.Print("NeutrinoEnergy.png");
canv.Print("NeutrinoEnergy.eps");
"""