diff --git a/distribsForHNLandBG.py b/distribsForHNLandBG.py new file mode 100644 index 0000000..9c63993 --- /dev/null +++ b/distribsForHNLandBG.py @@ -0,0 +1,124 @@ +import ROOT as r +import numpy as np +import funcsByBarbara as tools + +r.gROOT.SetBatch(r.kTRUE) +r.gROOT.ProcessLine(".L lhcbStyle.C+") +plots = {} + +def addPlot(key, x_axis, tree, string, cuts=None): + if not cuts: + tree.Draw(string+'>>h1') + else: + tree.Draw(string+'>>h1', cuts) + plots[key] = r.gDirectory.Get('h1') + plots[key].GetXaxis().SetTitle(x_axis) + plots[key].SetTitle(key) + plots[key].SetName(key) + c1 = r.TCanvas(); c1.cd() + plots[key].Draw() + c1.Print('accPlots/%s.pdf'%key) + c1.Print('accPlots/%s.png'%key) + + +fpimu = r.TFile('../Geraldine/MUPI/ShipAna.root','read') +pimu = fpimu.Get('ShipAna') +pimu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root') +fmumunu = r.TFile('../Geraldine/MUMUNU/ShipAna.root','read') +mumunu = fmumunu.Get('ShipAna') +mumunu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root') +fbg = r.TFile('../Geraldine/KLONG/ShipAna.root','read') +bg = fbg.Get('ShipAna') +bg_geo = tools.searchForNodes2_xyz_dict('../Geraldine/KLONG/geofile_full.10.0.Genie-TGeant4.root') + +studies = { + 'pimu': {'data':pimu, 'geo':pimu_geo}, + 'mumunu': {'data':mumunu, 'geo':mumunu_geo}, + 'bg': {'data':bg, 'geo':bg_geo} +} + +# Cuts: +# - z_vertex in [vol[0]+5m, straw[0]] +# - has track converged +# - chi2/ndof < 10 +# - DOCA < 10 cm +# - has no veto? + +def fit_converged(event): + flags = event.DaughtersFitConverged + if np.product(flags): + return True + return False + +def chi2_cut(event, chi2max): + chi2 = event.DaughtersChi2 + n = event.DaughtersNPoints + for i in xrange(len(n)): + if chi2[i]/n[i] > chi2max: + return False + return True + +def doca_cut(event, doca_max): + if event.DOCA < doca_max: + return True + return False + +def z_cut(event, study): + z = event.vtxz + minz = studies[study]['geo']['lidT1I_1']['z']['pos'] + studies[study]['geo']['lidT1I_1']['z']['dim'] + 500. + maxz = studies[study]['geo']['lidT6I_1']['z']['pos'] - studies[study]['geo']['lidT6I_1']['z']['dim'] + print minz, maxz + if minz < z < maxz: + return True + return False + +def nothing_in_veto5(event): + if len(event.veto5_x) > 0: + return False + return True + +def nothing_in_liquidscint(event): + if len(event.liquidscint_x) > 0: + return False + return True + +for s in studies: + print 'Reading %s...'%s + t = studies[s]['data'] + ntot = t.GetEntriesFast() + nAfterCuts = [0]*6 + for event in t: + if fit_converged(event): + nAfterCuts[0] += 1 + if chi2_cut(event,25): + nAfterCuts[1] +=1 + if doca_cut(event, 50.): + nAfterCuts[2] += 1 + if z_cut(event, s): + nAfterCuts[3] += 1 + if nothing_in_veto5(event): + nAfterCuts[4] += 1 + if nothing_in_liquidscint(event): + nAfterCuts[5] += 1 + print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot + addPlot(s+'-0cuts-IP', 'Impact parameter to target [cm]', t, 'IP0') + addPlot(s+'-0cuts-Mass', 'Reco inv. mass', t, 'Mass' ) + addPlot(s+'-0cuts-dPt', 'Tracks Pt', t, 'DaughtersPt' ) + addPlot(s+'-0cuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' ) + cut = "DaughtersChi2<25*DaughtersNPoints" + addPlot(s+'-chi2cut-IP', 'Impact parameter to target [cm]', t, 'IP0', cut) + addPlot(s+'-chi2cut-Mass', 'Reco inv. mass', t, 'Mass' , cut) + addPlot(s+'-chi2cut-dPt', 'Tracks Pt', t, 'DaughtersPt' , cut) + addPlot(s+'-chi2cut-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut) + cut += " && DOCA < 50" + addPlot(s+'-chi2andDOCAcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut) + addPlot(s+'-chi2andDOCAcuts-Mass', 'Reco inv. mass', t, 'Mass' , cut) + addPlot(s+'-chi2andDOCAcuts-dPt', 'Tracks Pt', t, 'DaughtersPt' , cut) + addPlot(s+'-chi2andDOCAcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut) + cut += " && vtxz > %s && vtxz < %s"%(studies[s]['geo']['lidT1I_1']['z']['pos'] + studies[s]['geo']['lidT1I_1']['z']['dim'] + 500., + studies[s]['geo']['lidT6I_1']['z']['pos'] - studies[s]['geo']['lidT6I_1']['z']['dim']) + addPlot(s+'-chi2andDOCAandZcuts-IP', 'Impact parameter to target [cm]', t, 'IP0', cut) + addPlot(s+'-chi2andDOCAandZcuts-Mass', 'Reco inv. mass', t, 'Mass' , cut) + addPlot(s+'-chi2andDOCAandZcuts-dPt', 'Tracks Pt', t, 'DaughtersPt' , cut) + addPlot(s+'-chi2andDOCAandZcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut) + diff --git a/elena.py b/elena.py index 58743d4..2d46678 100644 --- a/elena.py +++ b/elena.py @@ -1,29 +1,31 @@ import os -os.system("cd "+os.path.expandvars("$FAIRSHIP")+"/../FairShipRun/") -os.system(". config.sh") -os.system("cd "+os.path.expandvars("$FAIRSHIP")+"/../Barbara") +#os.system("cd "+os.path.expandvars("$FAIRSHIP")+"/../FairShipRun/") +#os.system(". config.sh") +#os.system("cd "+os.path.expandvars("$FAIRSHIP")+"/../Barbara") +# +## Make debug ntuple +#if not os.path.isfile("nuNutple_debug691.root"): +# os.system("python ntuple_nuVeto.py") +# +## For my output +#if not os.path.exists("elena_plots"): +# os.system("mkdir elena_plots") +# +## Import amenities +#from analyseNtupleForVeto import * +# +## Reminder +#print "" +#print "Interaction Elements:" +#print "0\t Last passive element of Opera" +#print "1\t Entrance window of the vacuum tank" +#print "2\t Along the vacuum tank" +#print "3\t Along Opera" +#print "4\t Between windows (in the liquid scintillator)" +#print "" -# Make debug ntuple -if not os.path.isfile("nuNutple_debug691.root"): - os.system("python ntuple_nuVeto.py") - -# For my output -if not os.path.exists("elena_plots"): - os.system("mkdir elena_plots") - -# Import amenities -from analyseNtupleForVeto import * - -# Reminder -print "" -print "Interaction Elements:" -print "0\t Last passive element of Opera" -print "1\t Entrance window of the vacuum tank" -print "2\t Along the vacuum tank" -print "3\t Along Opera" -print "4\t Between windows (in the liquid scintillator)" -print "" +from lookAtGeo import * def findDimentionBoxElement(node): sh = node.GetVolume().GetShape() @@ -86,124 +88,128 @@ d[n.GetName()]['r']['dim'] = 0. return d -# Search nodes -nodesDict = searchForNodes2_xyz_dict('geofile_full.10.0.Genie-TGeant4.root') -# Sort nodes by z -from collections import OrderedDict -sortedNodes = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['pos'])) -def printByZ(): - # Print nodes by increasing z - print - print - print "Nodes by z:" - for node in sortedNodes.keys(): - print node, '\t\t', (sortedNodes[node]['z']['pos']-sortedNodes[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodes[node]['z']['dim']/100., ' m' +if __name__ == '__main__': -# Sort nodes by dimension along z (the walls of the vacuum tank should be ~60 m long) -sortedNodesByZdim = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['dim'])) - -def printByZdim(): - # Print nodes by increasing dimension along z - print - print - print "Nodes by z dimension:" - for node in sortedNodesByZdim.keys(): - print node, '\t\t', (sortedNodesByZdim[node]['z']['pos']-sortedNodesByZdim[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodesByZdim[node]['z']['dim']/100., ' m' - -# Custom sorting -def sortNodesBy(coord='z', what='pos', printdict=True): - d = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1][coord][what])) - if printdict: + # Search nodes + nodesDict = searchForNodes2_xyz_dict('geofile_full.10.0.Genie-TGeant4.root') + # Sort nodes by z + from collections import OrderedDict + sortedNodes = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['pos'])) + + def printByZ(): + # Print nodes by increasing z print print - print "Nodes by %s %s:"%(coord, what) - for node in d.keys(): - print node, '\t\t', (d[node][coord]['pos']-d[node][coord]['dim'])/100., ' m', '\t', 2*d[node][coord]['dim']/100., ' m' - return d - - - -# Load geometry (again?) -#r = loadGeometry('geofile_full.10.0.Genie-TGeant4.root') -#fGeo = r['fGeo'] -#someOtherNodes = ["Tr%s_%s"%(i,i) for i in xrange(1,5)] -#someOtherNodes += ["volLateralS1_d_1", "volLateralS2_d_1", "volLateralS1_c_1", "volLateralS2_c_1"] -#someOtherNodes += ["Target%s_1"%i for i in (0,1,2,3,5)] - -# Make plots -from ROOT import * -cSaver = [] -cSaver.append(TCanvas("entrance_windows","entrance_windows")) -t2.Draw("startY_nu:startX_nu","isPrimary==1 && (interactionElement==2 || interactionElement==4)") - -# Veto -vetoIn = nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'] -vetoOut = nodesDict['Veto_5']['z']['pos'] + nodesDict['Veto_5']['z']['dim'] -cSaver.append(TCanvas("veto","veto")) -t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(vetoIn,vetoOut)) -grVeto = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) -grVeto.SetMarkerColor(kGreen) -grVeto.SetMarkerStyle(6) - -# Decay volume (big one) -cSaver.append(TCanvas("decay_volume_big","decay_volume_big")) -volIn = nodesDict['T2Startplate_1']['z']['pos'] + nodesDict['T2Startplate_1']['z']['dim'] -volOut = nodesDict['lidT6I_1']['z']['pos'] - nodesDict['lidT6I_1']['z']['dim'] -print -print "\tBig volume: in = %s cm, out = %s cm"%(volIn, volOut) -print -t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(volIn,volOut)) -grBig = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) -grBig.SetMarkerColor(kBlue) -grBig.SetMarkerStyle(6) - -# Decay volume (small one) -cSaver.append(TCanvas("decay_volume_small","decay_volume_small")) -volInSmall = nodesDict['lidT1I_1']['z']['pos'] + nodesDict['lidT1I_1']['z']['dim'] -#volOutSmall = nodesDict['T1Endplate_1']['z']['pos'] - nodesDict['T1Endplate_1']['z']['dim'] -volOutSmall = nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'] -print -print "\tSmall volume: in = %s cm, out = %s cm"%(volInSmall, volOutSmall) -print -print ' veto ', nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'], ' endplate ', nodesDict['T1Endplate_1']['z']['pos'] - nodesDict['T1Endplate_1']['z']['dim'] -#volOutSmall = volInSmall+500. -t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(volInSmall,volOutSmall)) -grSmall = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) -grSmall.SetMarkerColor(kRed) -grSmall.SetMarkerStyle(6) -cSaver.append(TCanvas('whole_decay_volume','whole_decay_volume')) -grBig.Draw('ap') -grSmall.Draw('p same') -grVeto.Draw('p same') - -# Save plots -for c in cSaver: - c.Print('elena_plots/'+c.GetName()+'.eps') - -# Elements of the vacuum tank: -print -endprevious = None -startnew = None -for i in xrange(1,7): - print 'T'+str(i)+':' - if 'T'+str(i)+'Startplate_1' in nodesDict: - print '\tStartplate (z = %s to %s)'%(nodesDict['T'+str(i)+'Startplate_1']['z']['pos']-nodesDict['T'+str(i)+'Startplate_1']['z']['dim'], nodesDict['T'+str(i)+'Startplate_1']['z']['pos']+nodesDict['T'+str(i)+'Startplate_1']['z']['dim']) - startnew = nodesDict['T'+str(i)+'Startplate_1']['z']['pos']-nodesDict['T'+str(i)+'Startplate_1']['z']['dim'] - else: - startnew = nodesDict['T'+str(i)+'I_1']['z']['pos']-nodesDict['T'+str(i)+'I_1']['z']['dim'] - if endprevious: - print '\tDelta z (%s to %s) = '%(i-1, i) + str(startnew-endprevious) - print '\tI (z = %s to %s)'%(nodesDict['T'+str(i)+'I_1']['z']['pos']-nodesDict['T'+str(i)+'I_1']['z']['dim'], nodesDict['T'+str(i)+'I_1']['z']['pos']+nodesDict['T'+str(i)+'I_1']['z']['dim']) - print '\tO (z = %s to %s)'%(nodesDict['T'+str(i)+'O_1']['z']['pos']-nodesDict['T'+str(i)+'O_1']['z']['dim'], nodesDict['T'+str(i)+'O_1']['z']['pos']+nodesDict['T'+str(i)+'O_1']['z']['dim']) - print '\t\tx diff: ' + str(nodesDict['T'+str(i)+'O_1']['x']['dim']- nodesDict['T'+str(i)+'I_1']['x']['dim']) - print '\t\ty diff: ' + str(nodesDict['T'+str(i)+'O_1']['y']['dim']- nodesDict['T'+str(i)+'I_1']['y']['dim']) - print '\t\tr diff: ' + str(nodesDict['T'+str(i)+'O_1']['r']['dim']- nodesDict['T'+str(i)+'I_1']['r']['dim']) - if 'T'+str(i)+'Endplate_1' in nodesDict: - print '\tEndplate (z = %s to %s)'%(nodesDict['T'+str(i)+'Endplate_1']['z']['pos']-nodesDict['T'+str(i)+'Endplate_1']['z']['dim'], nodesDict['T'+str(i)+'Endplate_1']['z']['pos']+nodesDict['T'+str(i)+'Endplate_1']['z']['dim']) - endprevious = nodesDict['T'+str(i)+'Endplate_1']['z']['pos']+nodesDict['T'+str(i)+'Endplate_1']['z']['dim'] - else: - endprevious = nodesDict['T'+str(i)+'I_1']['z']['pos']+nodesDict['T'+str(i)+'I_1']['z']['dim'] - -print 'Lenght of the whole decay volume (small + big): ' + str(nodesDict['lidT6I_1']['z']['pos'] + nodesDict['lidT6I_1']['z']['dim'] - (nodesDict['T2Startplate_1']['z']['pos'] - nodesDict['T2Startplate_1']['z']['dim'])) + print "Nodes by z:" + for node in sortedNodes.keys(): + print node, '\t\t', (sortedNodes[node]['z']['pos']-sortedNodes[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodes[node]['z']['dim']/100., ' m' + + # Sort nodes by dimension along z (the walls of the vacuum tank should be ~60 m long) + sortedNodesByZdim = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1]['z']['dim'])) + + def printByZdim(): + # Print nodes by increasing dimension along z + print + print + print "Nodes by z dimension:" + for node in sortedNodesByZdim.keys(): + print node, '\t\t', (sortedNodesByZdim[node]['z']['pos']-sortedNodesByZdim[node]['z']['dim'])/100., ' m', '\t', 2*sortedNodesByZdim[node]['z']['dim']/100., ' m' + + # Custom sorting + def sortNodesBy(coord='z', what='pos', printdict=True): + d = OrderedDict(sorted(nodesDict.items(), key = lambda item: item[1][coord][what])) + if printdict: + print + print + print "Nodes by %s %s:"%(coord, what) + for node in d.keys(): + print node, '\t\t', (d[node][coord]['pos']-d[node][coord]['dim'])/100., ' m', '\t', 2*d[node][coord]['dim']/100., ' m' + return d + + + + # Load geometry (again?) + #r = loadGeometry('geofile_full.10.0.Genie-TGeant4.root') + #fGeo = r['fGeo'] + #someOtherNodes = ["Tr%s_%s"%(i,i) for i in xrange(1,5)] + #someOtherNodes += ["volLateralS1_d_1", "volLateralS2_d_1", "volLateralS1_c_1", "volLateralS2_c_1"] + #someOtherNodes += ["Target%s_1"%i for i in (0,1,2,3,5)] + + + # Make plots + from ROOT import * + cSaver = [] + cSaver.append(TCanvas("entrance_windows","entrance_windows")) + t2.Draw("startY_nu:startX_nu","isPrimary==1 && (interactionElement==2 || interactionElement==4)") + + # Veto + vetoIn = nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'] + vetoOut = nodesDict['Veto_5']['z']['pos'] + nodesDict['Veto_5']['z']['dim'] + cSaver.append(TCanvas("veto","veto")) + t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(vetoIn,vetoOut)) + grVeto = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) + grVeto.SetMarkerColor(kGreen) + grVeto.SetMarkerStyle(6) + + # Decay volume (big one) + cSaver.append(TCanvas("decay_volume_big","decay_volume_big")) + volIn = nodesDict['T2Startplate_1']['z']['pos'] + nodesDict['T2Startplate_1']['z']['dim'] + volOut = nodesDict['lidT6I_1']['z']['pos'] - nodesDict['lidT6I_1']['z']['dim'] + print + print "\tBig volume: in = %s cm, out = %s cm"%(volIn, volOut) + print + t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(volIn,volOut)) + grBig = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) + grBig.SetMarkerColor(kBlue) + grBig.SetMarkerStyle(6) + + # Decay volume (small one) + cSaver.append(TCanvas("decay_volume_small","decay_volume_small")) + volInSmall = nodesDict['lidT1I_1']['z']['pos'] + nodesDict['lidT1I_1']['z']['dim'] + #volOutSmall = nodesDict['T1Endplate_1']['z']['pos'] - nodesDict['T1Endplate_1']['z']['dim'] + volOutSmall = nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'] + print + print "\tSmall volume: in = %s cm, out = %s cm"%(volInSmall, volOutSmall) + print + print ' veto ', nodesDict['Veto_5']['z']['pos'] - nodesDict['Veto_5']['z']['dim'], ' endplate ', nodesDict['T1Endplate_1']['z']['pos'] - nodesDict['T1Endplate_1']['z']['dim'] + #volOutSmall = volInSmall+500. + t2.Draw("startY_nu:startX_nu","isPrimary==1 && interactionElement==2 && startZ_nu>%s && startZ_nu<%s"%(volInSmall,volOutSmall)) + grSmall = TGraph(t2.GetSelectedRows(), t2.GetV2(), t2.GetV1()) + grSmall.SetMarkerColor(kRed) + grSmall.SetMarkerStyle(6) + cSaver.append(TCanvas('whole_decay_volume','whole_decay_volume')) + grBig.Draw('ap') + grSmall.Draw('p same') + grVeto.Draw('p same') + + # Save plots + for c in cSaver: + c.Print('elena_plots/'+c.GetName()+'.eps') + + # Elements of the vacuum tank: + print + endprevious = None + startnew = None + for i in xrange(1,7): + print 'T'+str(i)+':' + if 'T'+str(i)+'Startplate_1' in nodesDict: + print '\tStartplate (z = %s to %s)'%(nodesDict['T'+str(i)+'Startplate_1']['z']['pos']-nodesDict['T'+str(i)+'Startplate_1']['z']['dim'], nodesDict['T'+str(i)+'Startplate_1']['z']['pos']+nodesDict['T'+str(i)+'Startplate_1']['z']['dim']) + startnew = nodesDict['T'+str(i)+'Startplate_1']['z']['pos']-nodesDict['T'+str(i)+'Startplate_1']['z']['dim'] + else: + startnew = nodesDict['T'+str(i)+'I_1']['z']['pos']-nodesDict['T'+str(i)+'I_1']['z']['dim'] + if endprevious: + print '\tDelta z (%s to %s) = '%(i-1, i) + str(startnew-endprevious) + print '\tI (z = %s to %s)'%(nodesDict['T'+str(i)+'I_1']['z']['pos']-nodesDict['T'+str(i)+'I_1']['z']['dim'], nodesDict['T'+str(i)+'I_1']['z']['pos']+nodesDict['T'+str(i)+'I_1']['z']['dim']) + print '\tO (z = %s to %s)'%(nodesDict['T'+str(i)+'O_1']['z']['pos']-nodesDict['T'+str(i)+'O_1']['z']['dim'], nodesDict['T'+str(i)+'O_1']['z']['pos']+nodesDict['T'+str(i)+'O_1']['z']['dim']) + print '\t\tx diff: ' + str(nodesDict['T'+str(i)+'O_1']['x']['dim']- nodesDict['T'+str(i)+'I_1']['x']['dim']) + print '\t\ty diff: ' + str(nodesDict['T'+str(i)+'O_1']['y']['dim']- nodesDict['T'+str(i)+'I_1']['y']['dim']) + print '\t\tr diff: ' + str(nodesDict['T'+str(i)+'O_1']['r']['dim']- nodesDict['T'+str(i)+'I_1']['r']['dim']) + if 'T'+str(i)+'Endplate_1' in nodesDict: + print '\tEndplate (z = %s to %s)'%(nodesDict['T'+str(i)+'Endplate_1']['z']['pos']-nodesDict['T'+str(i)+'Endplate_1']['z']['dim'], nodesDict['T'+str(i)+'Endplate_1']['z']['pos']+nodesDict['T'+str(i)+'Endplate_1']['z']['dim']) + endprevious = nodesDict['T'+str(i)+'Endplate_1']['z']['pos']+nodesDict['T'+str(i)+'Endplate_1']['z']['dim'] + else: + endprevious = nodesDict['T'+str(i)+'I_1']['z']['pos']+nodesDict['T'+str(i)+'I_1']['z']['dim'] + + print 'Lenght of the whole decay volume (small + big): ' + str(nodesDict['lidT6I_1']['z']['pos'] + nodesDict['lidT6I_1']['z']['dim'] - (nodesDict['T2Startplate_1']['z']['pos'] - nodesDict['T2Startplate_1']['z']['dim'])) diff --git a/elena_ShipAna.py b/elena_ShipAna.py new file mode 100644 index 0000000..e0a2f7b --- /dev/null +++ b/elena_ShipAna.py @@ -0,0 +1,489 @@ +# example for accessing smeared hits and fitted tracks +import os,sys,getopt,subprocess,gc + +sys.path.append("../../FairShipTools") +import funcsByBarbara as tools + +#os.chdir("../../FairShipRun") +#os.system("bash config.sh") +#os.chdir("../FairShip/macro") +## Tanto non funziona +#subprocess.call("cd ../../FairShipRun", shell=True) +#subprocess.call(". config.sh", shell=True) +#subprocess.call("cd ../FairShip/macro", shell=True) + +import ROOT +import rootUtils as ut +import shipunit as u +from ShipGeoConfig import ConfigRegistry + +PDG = ROOT.TDatabasePDG.Instance() +inputFile = None +dy = None +nEvents = 99999 +theHNLcode = 9900015 +signal_file = False +bg_file = False +file_type = None + +try: + opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:A:Y:i", ["nEvents=", "type="]) +except getopt.GetoptError: + # print help information and exit: + print ' enter file name' + sys.exit() +for o, a in opts: + if o in ("-f"): + inputFile = a + if o in ("-o"): + outputFile = a + if o in ("-Y"): + dy = float(a) + if o in ("-n", "--nEvents="): + nEvents = int(a) + if o in ("-t", "--type="): + file_type = str(a) + +if not file_type: + print " please select file type (sig or bg)" + sys.exit() + +if 'sig' in file_type: + signal_file = True +elif 'bg' in file_type: + bg_file = True +else: + print " please select file type (sig or bg)" + sys.exit() + +print +print "\tFile type is: "+file_type +print + +if not outputFile: + outputFile = "ShipAna.root" + +# If directory of output file does not exist, create it (depth=3) +if not os.path.exists(os.path.dirname(outputFile)): + if not os.path.exists(os.path.dirname(os.path.dirname(outputFile))): + if not os.path.exists(os.path.dirname(os.path.dirname(os.path.dirname(outputFile)))): + os.system("mkdir "+os.path.dirname(os.path.dirname(os.path.dirname(outputFile)))) + os.system("mkdir "+os.path.dirname(os.path.dirname(outputFile))) + os.system("mkdir "+os.path.dirname(outputFile)) + +if not dy: + # try to extract from input file name + tmp = inputFile.replace(os.path.dirname(inputFile),'') + tmp = tmp.split('.') + try: + dy = float( tmp[1]+'.'+tmp[2] ) + except: + dy = None +#else: +# inputFile = 'ship.'+str(dy)+'.Pythia8-TGeant4_rec.root' + +if not inputFile: + inputFile = 'ship.'+str(dy)+'.Pythia8-TGeant4_rec.root' + +f = ROOT.TFile(inputFile) +sTree = f.cbmsim + +# init geometry and mag. field +ShipGeo = ConfigRegistry.loadpy("$FAIRSHIP/geometry/geometry_config.py", Yheight = dy ) +# -----Create geometry---------------------------------------------- +import shipDet_conf +run = ROOT.FairRunSim() +modules = shipDet_conf.configure(run,ShipGeo) + +tgeom = ROOT.TGeoManager("Geometry", "Geane geometry") +geofile = inputFile.replace('ship.','geofile_full.').replace('_rec.','.') +gMan = tgeom.Import(geofile) +geoMat = ROOT.genfit.TGeoMaterialInterface() +ROOT.genfit.MaterialEffects.getInstance().init(geoMat) + +bfield = ROOT.genfit.BellField(ShipGeo.Bfield.max, ShipGeo.Bfield.z,2) +fM = ROOT.genfit.FieldManager.getInstance() +fM.init(bfield) + +ev = ut.PyListOfLeaves() +leaves = sTree.GetListOfLeaves() +names = ut.setAttributes(ev, leaves) + +# Get geometry positions +nodes = tools.searchForNodes2_xyz_dict(geofile) +scintTankW = [nodes["lidT1lisci_1"]['z']['pos']-nodes["lidT1lisci_1"]['z']['dim'], nodes["lidT1lisci_1"]['z']['pos']+nodes["lidT1lisci_1"]['z']['dim']] +scintWPos = [scintTankW] +trackStationsPos = [ [nodes["Tr%s_%s"%(i,i)]['z']['pos']-nodes["Tr%s_%s"%(i,i)]['z']['dim'], nodes["Tr%s_%s"%(i,i)]['z']['pos']+nodes["Tr%s_%s"%(i,i)]['z']['dim']] for i in xrange(1,5)] +muonStationsPos = [ [nodes["muondet%s_1"%(i)]['z']['pos']-nodes["muondet%s_1"%(i)]['z']['dim'], nodes["muondet%s_1"%(i)]['z']['pos']+nodes["muondet%s_1"%(i)]['z']['dim']] for i in xrange(0,4)] +strawVetoPos = [ [nodes["Veto_5"]['z']['pos']-nodes["Veto_5"]['z']['dim'], nodes["Veto_5"]['z']['pos']+nodes["Veto_5"]['z']['dim']] ] +# ecal points are before the start of the ecal, systematically. +ecalPos = [ [nodes['Ecal_1']['z']['pos']-2.*nodes['Ecal_1']['z']['dim'], nodes['Ecal_1']['z']['pos']+nodes['Ecal_1']['z']['dim']] ] +hcalPos = [ [nodes['Hcal_1']['z']['pos']-nodes['Hcal_1']['z']['dim'], nodes['Hcal_1']['z']['pos']+nodes['Hcal_1']['z']['dim']] ] +volume = [nodes["lidT1O_1"]['z']['pos']-nodes["lidT6O_1"]['z']['dim'],nodes["lidT6O_1"]['z']['pos']-nodes["lidT6O_1"]['z']['dim']] +vetoWall = [ [volume[0], nodes['Tr1_1']['z']['pos']-nodes['Tr1_1']['z']['dim']] ] +print '\t stored some geometry nodes' +#print 'ecal pos: ', ecalPos[0] + +h = {} +#ut.bookHist(h,'delPOverP','delP / P',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'delPtOverP','delPt / P',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'delPtOverPt','delPt / Pt',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'delPOverP2','delP / P chi2/nmeas<25',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'delPOverPz','delPz / Pz',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'delPOverP2z','delPz / Pz chi2/nmeas<25',100,0.,50.,100,-0.5,0.5) +#ut.bookHist(h,'chi2','chi2/nmeas after trackfit',100,0.,100.) +#ut.bookHist(h,'IP','Impact Parameter',100,0.,10.) +#ut.bookHist(h,'meas','number of measurements',40,-0.5,39.5) +#ut.bookHist(h,'measVSchi2','number of measurements vs chi2/meas',40,-0.5,39.5,100,0.,100.) +#ut.bookHist(h,'distu','distance to wire',100,0.,1.) +#ut.bookHist(h,'distv','distance to wire',100,0.,1.) +#ut.bookHist(h,'disty','distance to wire',100,0.,1.) +#ut.bookHist(h,'meanhits','mean number of hits / track',50,-0.5,49.5) +#ut.bookHist(h,'Pt','Reconstructed transverse momentum',100,0.,40.) +ut.bookHist(h,'Candidate-DOCA','DOCA between the two tracks',300,0.,50.) +ut.bookHist(h,'Candidate-IP0','Impact Parameter to target',200,0.,300.) +ut.bookHist(h,'Candidate-IP0/mass','Impact Parameter to target vs mass',100,0.,2.,100,0.,100.) +ut.bookHist(h,'Candidate-Mass','reconstructed Mass',100,0.,4.) +ut.bookHist(h,'Candidate-Pt','reconstructed Pt',100,0.,40.) +ut.bookHist(h,'Candidate-vtxx','X position of reconstructed vertex',150,-4000.,4000.) +ut.bookHist(h,'Candidate-vtxy','Y position of reconstructed vertex',150,-4000.,4000.) +ut.bookHist(h,'Candidate-vtxz','Z position of reconstructed vertex',150,-4000.,4000.) +ut.bookHist(h,'CandidateDaughters-Pt','reconstructed Pt of candidate daughters',100,0.,40.) +ut.bookHist(h,'CandidateDaughters-chi2','chi2/nmeas after trackfit',100,0.,100.) + +def myVertex(t1,t2,PosDir): + # closest distance between two tracks + V=0 + for i in range(3): V += PosDir[t1][1](i)*PosDir[t2][1](i) + S1=0 + for i in range(3): S1 += (PosDir[t1][0](i)-PosDir[t2][0](i))*PosDir[t1][1](i) + S2=0 + for i in range(3): S2 += (PosDir[t1][0](i)-PosDir[t2][0](i))*PosDir[t2][1](i) + l = (S2-S1*V)/(1-V*V) + x2 = PosDir[t2][0](0)+l*PosDir[t2][1](0) + y2 = PosDir[t2][0](1)+l*PosDir[t2][1](1) + z2 = PosDir[t2][0](2)+l*PosDir[t2][1](2) + x1 = PosDir[t1][0](0)+l*PosDir[t1][1](0) + y1 = PosDir[t1][0](1)+l*PosDir[t1][1](1) + z1 = PosDir[t1][0](2)+l*PosDir[t1][1](2) + dist = ROOT.TMath.Sqrt((x1-x2)**2+(y1-y2)**2+(z1-z2)**2) + return (x1+x2)/2.,(y1+y2)/2.,(z1+z2)/2.,dist + +def fitSingleGauss(x,ba=None,be=None): + name = 'myGauss_'+x + myGauss = h[x].GetListOfFunctions().FindObject(name) + if not myGauss: + if not ba : ba = h[x].GetBinCenter(1) + if not be : be = h[x].GetBinCenter(h[x].GetNbinsX()) + bw = h[x].GetBinWidth(1) + mean = h[x].GetMean() + sigma = h[x].GetRMS() + norm = h[x].GetEntries()*0.3 + myGauss = ROOT.TF1(name,'[0]*'+str(bw)+'/([2]*sqrt(2*pi))*exp(-0.5*((x-[1])/[2])**2)+[3]',4) + myGauss.SetParameter(0,norm) + myGauss.SetParameter(1,mean) + myGauss.SetParameter(2,sigma) + myGauss.SetParameter(3,1.) + myGauss.SetParName(0,'Signal') + myGauss.SetParName(1,'Mean') + myGauss.SetParName(2,'Sigma') + myGauss.SetParName(3,'bckgr') + h[x].Fit(myGauss,'','',ba,be) + + +def makePlots(): + ut.bookCanvas(h,key='strawanalysis',title='Distance to wire and mean nr of hits',nx=1200,ny=600,cx=2,cy=1) + cv = h['strawanalysis'].cd(1) + h['disty'].Draw() + h['distu'].Draw('same') + h['distv'].Draw('same') + cv = h['strawanalysis'].cd(2) + h['meanhits'].Draw() + ut.bookCanvas(h,key='fitresults',title='Fit Results',nx=1600,ny=1200,cx=2,cy=2) + cv = h['fitresults'].cd(1) + h['delPOverPz'].Draw('box') + cv = h['fitresults'].cd(2) + cv.SetLogy(1) + h['chi2'].Draw() + cv = h['fitresults'].cd(3) + h['delPOverPz_proj'] = h['delPOverPz'].ProjectionY() + ROOT.gStyle.SetOptFit(11111) + h['delPOverPz_proj'].Draw() + h['delPOverPz_proj'].Fit('gaus') + cv = h['fitresults'].cd(4) + h['delPOverP2z_proj'] = h['delPOverP2z'].ProjectionY() + h['delPOverP2z_proj'].Draw() + fitSingleGauss('delPOverP2z_proj') + h['fitresults'].Print('fitresults.gif') + ut.bookCanvas(h,key='fitresults2',title='Fit Results',nx=1600,ny=1200,cx=2,cy=2) + print 'finished with first canvas' + cv = h['fitresults2'].cd(1) + h['Doca'].Draw() + cv = h['fitresults2'].cd(2) + h['IP0'].Draw() + cv = h['fitresults2'].cd(3) + h['Mass'].Draw() + fitSingleGauss('Mass') + cv = h['fitresults2'].cd(4) + h['IP0/mass'].Draw('box') + h['fitresults2'].Print('fitresults2.gif') + print 'finished making plots' + +# start event loop +def myEventLoop(N): + nEvents = min(sTree.GetEntries(),N) + for n in range(nEvents): + rc = sTree.GetEntry(n) + wg = sTree.MCTrack[0].GetWeight() + if not wg>0.: wg=1. + ## make some straw hit analysis + #hitlist = {} + #for ahit in sTree.strawtubesPoint: + # detID = ahit.GetDetectorID() + # top = ROOT.TVector3() + # bot = ROOT.TVector3() + # modules["Strawtubes"].StrawEndPoints(detID,bot,top) + # dw = ahit.dist2Wire() + # if abs(top.y())==abs(bot.y()): h['disty'].Fill(dw) + # if abs(top.y())>abs(bot.y()): h['distu'].Fill(dw) + # if abs(top.y())25: continue + h['delPOverP2'].Fill(Ptruth,delPOverP) + h['delPOverP2z'].Fill(Ptruth,delPOverPz) + # try measure impact parameter + trackDir = fittedState.getDir() + trackPos = fittedState.getPos() + vx = ROOT.TVector3() + mcPart.GetStartVertex(vx) + t = 0 + for i in range(3): t += trackDir(i)*(vx(i)-trackPos(i)) + dist = 0 + for i in range(3): dist += (vx(i)-trackPos(i)-t*trackDir(i))**2 + dist = ROOT.TMath.Sqrt(dist) + h['IP'].Fill(dist) + # loop over particles, 2-track combinations + # From Thomas: + # An object in sTree.Particles is made out of two fitted tracks. + # Except of the mass assignment, and fake pattern recognition (only + # hits from same MCTrack are used), no other MC truth. No requirement + # that the two tracks come from the same MCTrack ! GetDaughter allows to + # get back to the underlying tracks, see ShipReco.py: + # particle = ROOT.TParticle(9900015,0,-1,-1,t1,t2,HNL,vx) + for HNL in sTree.Particles: + if signal_file and not (HNL.GetPdgCode() == theHNLcode): continue + #if bg_file and (HNL.GetPdgCode() == theHNLcode): continue + t1,t2 = HNL.GetDaughter(0),HNL.GetDaughter(1) + PosDir = {} + for tr in [t1,t2]: + xx = sTree.FitTracks[tr].getFittedState() + PosDir[tr] = [xx.getPos(),xx.getDir()] + h['CandidateDaughters-Pt'].Fill(xx.getMom().Pt()) + xv,yv,zv,doca = myVertex(t1,t2,PosDir) + h['Candidate-DOCA'].Fill(doca) + h['Candidate-vtxx'].Fill(xv) + h['Candidate-vtxy'].Fill(yv) + h['Candidate-vtxz'].Fill(zv) + #if doca>5 : continue + HNLPos = ROOT.TLorentzVector() + HNL.ProductionVertex(HNLPos) + HNLMom = ROOT.TLorentzVector() + HNL.Momentum(HNLMom) + tr = ROOT.TVector3(0,0,ShipGeo.target.z0) + t = 0 + for i in range(3): t += HNLMom(i)/HNLMom.P()*(tr(i)-HNLPos(i)) + dist = 0 + for i in range(3): dist += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2 + dist = ROOT.TMath.Sqrt(dist) + h['Candidate-IP0'].Fill(dist) + h['Candidate-IP0/mass'].Fill(HNLMom.M(),dist) + h['Candidate-Mass'].Fill(HNLMom.M()) + h['Candidate-Pt'].Fill(HNLMom.Pt()) + +elenaTree = ROOT.TTree('ShipAna','ShipAna') +elenaTree, DaughtersPt = tools.AddVect(elenaTree, 'DaughtersPt', 'float') +elenaTree, DaughtersChi2 = tools.AddVect(elenaTree, 'DaughtersChi2', 'float') +elenaTree, DaughtersNPoints = tools.AddVect(elenaTree, 'DaughtersNPoints', 'int') +elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') +elenaTree, straw_x = tools.AddVect(elenaTree, 'straw_x', 'float') +elenaTree, straw_y = tools.AddVect(elenaTree, 'straw_y', 'float') +elenaTree, straw_z = tools.AddVect(elenaTree, 'straw_z', 'float') +elenaTree, muon_x = tools.AddVect(elenaTree, 'muon_x', 'float') +elenaTree, muon_y = tools.AddVect(elenaTree, 'muon_y', 'float') +elenaTree, muon_z = tools.AddVect(elenaTree, 'muon_z', 'float') +elenaTree, ecal_x = tools.AddVect(elenaTree, 'ecal_x', 'float') +elenaTree, ecal_y = tools.AddVect(elenaTree, 'ecal_y', 'float') +elenaTree, ecal_z = tools.AddVect(elenaTree, 'ecal_z', 'float') +elenaTree, hcal_x = tools.AddVect(elenaTree, 'hcal_x', 'float') +elenaTree, hcal_y = tools.AddVect(elenaTree, 'hcal_y', 'float') +elenaTree, hcal_z = tools.AddVect(elenaTree, 'hcal_z', 'float') +elenaTree, veto5_x = tools.AddVect(elenaTree, 'veto5_x', 'float') +elenaTree, veto5_y = tools.AddVect(elenaTree, 'veto5_y', 'float') +elenaTree, veto5_z = tools.AddVect(elenaTree, 'veto5_z', 'float') +elenaTree, liquidscint_x = tools.AddVect(elenaTree, 'liquidscint_x', 'float') +elenaTree, liquidscint_y = tools.AddVect(elenaTree, 'liquidscint_y', 'float') +elenaTree, liquidscint_z = tools.AddVect(elenaTree, 'liquidscint_z', 'float') +elenaTree, DOCA = tools.AddVar(elenaTree, 'DOCA', 'float') +elenaTree, vtxx = tools.AddVar(elenaTree, 'vtxx', 'float') +elenaTree, vtxy = tools.AddVar(elenaTree, 'vtxy', 'float') +elenaTree, vtxz = tools.AddVar(elenaTree, 'vtxz', 'float') +elenaTree, IP0 = tools.AddVar(elenaTree, 'IP0', 'float') +elenaTree, Mass = tools.AddVar(elenaTree, 'Mass', 'float') +elenaTree, Pt = tools.AddVar(elenaTree, 'Pt', 'float') + +def elenaEventLoop(N): + nEvents = min(sTree.GetEntries(),N) + for n in range(nEvents): + rc = sTree.GetEntry(n) + key = -1 + # loop over particles, 2-track combinations + # From Thomas: + # An object in sTree.Particles is made out of two fitted tracks. + # Except of the mass assignment, and fake pattern recognition (only + # hits from same MCTrack are used), no other MC truth. No requirement + # that the two tracks come from the same MCTrack ! GetDaughter allows to + # get back to the underlying tracks, see ShipReco.py: + # particle = ROOT.TParticle(9900015,0,-1,-1,t1,t2,HNL,vx) + for HNL in sTree.Particles: + if signal_file and not (HNL.GetPdgCode() == theHNLcode): continue + #if bg_file and (HNL.GetPdgCode() == theHNLcode): continue + # Fill hit arrays + tools.PutToZero(straw_x); tools.PutToZero(straw_y); tools.PutToZero(straw_z) + tools.PutToZero(veto5_x); tools.PutToZero(veto5_y); tools.PutToZero(veto5_z) + tools.PutToZero(muon_x); tools.PutToZero(muon_y); tools.PutToZero(muon_z) + tools.PutToZero(ecal_x); tools.PutToZero(ecal_y); tools.PutToZero(ecal_z) + tools.PutToZero(hcal_x); tools.PutToZero(hcal_y); tools.PutToZero(hcal_z) + tools.PutToZero(liquidscint_x); tools.PutToZero(liquidscint_y); tools.PutToZero(liquidscint_z) + hasStraws, nothing = tools.wasFired(None, sTree.strawtubesPoint, trackStationsPos, pointsVects=[straw_x, straw_y, straw_z], Ethr=0.) + hasVeto5, nothing = tools.wasFired(None, sTree.strawtubesPoint, strawVetoPos, pointsVects=[veto5_x, veto5_y, veto5_z], Ethr=0.) + hasMuon, nothing = tools.wasFired(None, sTree.muonPoint, muonStationsPos, pointsVects=[muon_x, muon_y, muon_z], Ethr=0.) + hasEcal, ecalOverThr = tools.wasFired(None, sTree.EcalPoint, ecalPos, pointsVects=[ecal_x, ecal_y, ecal_z], Ethr=0.015) + hasHcal, hcalOverThr = tools.wasFired(None, sTree.HcalPoint, hcalPos, pointsVects=[hcal_x, hcal_y, hcal_z], Ethr=0.015) + hasliquidscint, liquidscintOverThr = tools.wasFired(None, sTree.vetoPoint, vetoWall, pointsVects=[liquidscint_x, liquidscint_y, liquidscint_z], Ethr=0.015) + # get "daughters" + t1,t2 = HNL.GetDaughter(0),HNL.GetDaughter(1) + PosDir = {} + tools.PutToZero(DaughtersFitConverged) + tools.PutToZero(DaughtersChi2) + tools.PutToZero(DaughtersPt) + tools.PutToZero(DaughtersNPoints) + for tr in [t1,t2]: + converged = 0 + x = sTree.FitTracks[tr] + xx = x.getFittedState() + if x.getFitStatus().isFitConverged(): + converged = 1 + tools.Push(DaughtersFitConverged, converged) + PosDir[tr] = [xx.getPos(),xx.getDir()] + h['CandidateDaughters-Pt'].Fill(xx.getMom().Pt()) + tools.Push(DaughtersPt, xx.getMom().Pt()) + h['CandidateDaughters-chi2'].Fill(x.getFitStatus().getChi2() / x.getNumPoints()) + tools.Push(DaughtersChi2, x.getFitStatus().getChi2()) + tools.Push(DaughtersNPoints, x.getNumPoints()) + xv,yv,zv,doca = myVertex(t1,t2,PosDir) + h['Candidate-DOCA'].Fill(doca) + tools.PutToZero(DOCA); tools.Push(DOCA, doca) + h['Candidate-vtxx'].Fill(xv) + tools.PutToZero(vtxx); tools.Push(vtxx, xv) + h['Candidate-vtxy'].Fill(yv) + tools.PutToZero(vtxy); tools.Push(vtxy, yv) + h['Candidate-vtxz'].Fill(zv) + tools.PutToZero(vtxz); tools.Push(vtxz, zv) + #if doca>5 : continue + HNLPos = ROOT.TLorentzVector() + HNL.ProductionVertex(HNLPos) + HNLMom = ROOT.TLorentzVector() + HNL.Momentum(HNLMom) + tr = ROOT.TVector3(0,0,ShipGeo.target.z0) + t = 0 + for i in range(3): t += HNLMom(i)/HNLMom.P()*(tr(i)-HNLPos(i)) + dist = 0 + for i in range(3): dist += (tr(i)-HNLPos(i)-t*HNLMom(i)/HNLMom.P())**2 + dist = ROOT.TMath.Sqrt(dist) + h['Candidate-IP0'].Fill(dist) + tools.PutToZero(IP0); tools.Push(IP0, dist) + h['Candidate-IP0/mass'].Fill(HNLMom.M(),dist) + h['Candidate-Mass'].Fill(HNLMom.M()) + tools.PutToZero(Mass); tools.Push(Mass, HNLMom.M()) + h['Candidate-Pt'].Fill(HNLMom.Pt()) + tools.PutToZero(Pt); tools.Push(Pt, HNLMom.Pt()) + elenaTree.SetDirectory(0) + elenaTree.Fill() + + +def HNLKinematics(): + ut.bookHist(h,'HNLmomNoW','momentum unweighted',100,0.,300.) + ut.bookHist(h,'HNLmom','momentum',100,0.,300.) + ut.bookHist(h,'HNLmom_recTracks','HNL momentum from reco tracks',100,0.,300.) + ut.bookHist(h,'HNLmomNoW_recTracks','HNL momentum unweighted from reco tracks',100,0.,300.) + for n in range(sTree.GetEntries()): + rc = sTree.GetEntry(n) + wg = sTree.MCTrack[1].GetWeight() + if not wg>0.: wg=1. + P = sTree.MCTrack[1].GetP() + h['HNLmom'].Fill(P,wg) + h['HNLmomNoW'].Fill(P) + for HNL in sTree.Particles: + t1,t2 = HNL.GetDaughter(0),HNL.GetDaughter(1) + for tr in [t1,t2]: + xx = sTree.FitTracks[tr].getFittedState() + Prec = xx.getMom().Mag() + h['HNLmom_recTracks'].Fill(Prec,wg) + h['HNLmomNoW_recTracks'].Fill(Prec) +# +def access2SmearedHits(): + key = 0 + for ahit in ev.SmearedHits.GetObject(): + print ahit[0],ahit[1],ahit[2],ahit[3],ahit[4],ahit[5],ahit[6] + # follow link to true MCHit + mchit = TrackingHits[key] + mctrack = MCTracks[mchit.GetTrackID()] + print mchit.GetZ(),mctrack.GetP(),mctrack.GetPdgCode() + key+=1 + +#myEventLoop(nEvents) +elenaEventLoop(nEvents) +#makePlots() +#HNLKinematics() +# output histograms +ut.writeHists(h,outputFile) +ofile = ROOT.TFile(outputFile, "update") +elenaTree.Write() +ofile.Close() +print "\tOutput saved to ", outputFile \ No newline at end of file diff --git a/funcsByBarbara.py b/funcsByBarbara.py new file mode 100644 index 0000000..26f6561 --- /dev/null +++ b/funcsByBarbara.py @@ -0,0 +1,69 @@ +from elena import * +from array import array +import ROOT + +def AddVect(t,name,vectType): + vect = ROOT.vector(vectType)() + t.Branch( name, vect ) + return t, vect + +def AddVar(t,name,varType): + var = array(varType[0],[-999]) + t.Branch(name,var,name+"/"+varType.upper()) + return t, var + +def PutToZero(var): + if not isinstance(var, array): + var.clear() + else: + var[0] = 0 + +def Push(leaf, value): + if not isinstance(leaf, array): + leaf.push_back(value) + else: + leaf[0] = value + + +def wasFired(indexKids, detPoints, detPos, pointsVects=None, Ethr=0): + def lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr): + for hit in detPoints: + if (indexKids is None) or (hit.GetTrackID() in indexKids): + #print RPChit.GetTrackID(), t_ID + # check if it is in one of the considered active layer + for pos in detPos: + if pos[0]<=hit.GetZ()<=pos[1]: + Eloss += hit.GetEnergyLoss() + #print pos, hit.GetZ() + if pointsVects is not None: + pointsVects[0].push_back(hit.GetX()) + pointsVects[1].push_back(hit.GetY()) + pointsVects[2].push_back(hit.GetZ()) + + flag = True + nstat+=1 + #eff_val = gRandom.Uniform(0.,1.) + #if eff_val<0.9: + # flag_eff = flag_eff or True + # nstat_eff+=1 + #else: + # flag_eff = flag_eff or False + flag_Ethr = Eloss>=Ethr + return flag, nstat, nstat_eff, Eloss, flag_Ethr + # Now in partKidTrackID I should have the trackID connected to my charged particles + #flag_eff = False + flag = False + nstat = 0 + nstat_eff = 0 + Eloss = 0 + flag,nstat,nstat_eff,Eloss,flag_Ethr = lookingForHits(detPoints,flag,nstat,nstat_eff,indexKids,Eloss,Ethr) + #if flag==False and checkOn: + #print "To be study event %s"%entry + #PrintEventPart(particles,pdg) + #PrintRPChits(detPoints,pdg) + #print + #print + #assert(False) + #return flag, flag_eff, nstat, nstat_eff, flag_Ethr + return flag, flag_Ethr + diff --git a/lhcbStyle.C b/lhcbStyle.C new file mode 100644 index 0000000..1915de3 --- /dev/null +++ b/lhcbStyle.C @@ -0,0 +1,200 @@ +// all users - please change the name of this file to lhcbStyle.C +// Commits to lhcbdocs svn of .C files are not allowed + +// Global scope variables +TStyle* lhcbStyle; // general lhcb style +TPaveText* lhcbName; // standard lhcb text for plot +TText* lhcbLabel; // style for Ttext +TLatex *lhcbLatex; //style for TLatex; + +// define names for colours +Int_t black=1; +Int_t red=2; +Int_t green=3; +Int_t blue=4; +Int_t yellow=5; +Int_t magenta=6; +Int_t cyan=7; +Int_t purple=9; + +void lhcbStyle() { + +//////////////////////////////////////////////////////////////////// +// PURPOSE: +// +// This macro defines a standard style for (black-and-white) +// "publication quality" LHCb ROOT plots. +// +// USAGE: +// +// Include the lines +// gROOT->ProcessLine(".L lhcbstyle.C"); +// lhcbStyle(); +// at the beginning of your root macro. +// +// Example usage is given in myPlot.C +// +// COMMENTS: +// +// Font: +// +// The font is chosen to be 62, this is helvetica-bold-r-normal with +// precision 2. +// +// "Landscape histograms": +// +// The style here is designed for more or less square plots. +// For longer histograms, or canvas with many pads, adjustements are needed. +// For instance, for a canvas with 1x5 histograms: +// TCanvas* c1 = new TCanvas("c1", "L0 muons", 600, 800); +// c1->Divide(1,5); +// Adaptions like the following will be needed: +// lhcbStyle->SetTickLength(0.05,"x"); +// lhcbStyle->SetTickLength(0.01,"y"); +// lhcbStyle->SetLabelSize(0.15,"x"); +// lhcbStyle->SetLabelSize(0.1,"y"); +// lhcbStyle->SetStatW(0.15); +// lhcbStyle->SetStatH(0.5); +// +// Authors: Thomas Schietinger, Andrew Powell, Chris Parkes +// Maintained by Editorial board member (currently Chris) +/////////////////////////////////////////////////////////////////// + +lhcbStyle= new TStyle("lhcbStyle","Standard LHCb plots style"); + +// use helvetica-bold-r-normal, precision 2 (rotatable) +Int_t lhcbFont = 62; +// line thickness +Double_t lhcbWidth = 3.00; + +// use plain black on white colors +lhcbStyle->SetFrameBorderMode(0); +lhcbStyle->SetCanvasBorderMode(0); +lhcbStyle->SetPadBorderMode(0); +lhcbStyle->SetPadColor(0); +lhcbStyle->SetCanvasColor(0); +lhcbStyle->SetStatColor(0); +lhcbStyle->SetPalette(1); + +// set the paper & margin sizes +lhcbStyle->SetPaperSize(20,26); +lhcbStyle->SetPadTopMargin(0.05); +lhcbStyle->SetPadRightMargin(0.05); // increase for colz plots +lhcbStyle->SetPadBottomMargin(0.16); +lhcbStyle->SetPadLeftMargin(0.14); + +// use large fonts +lhcbStyle->SetTextFont(lhcbFont); +lhcbStyle->SetTextSize(0.08); +lhcbStyle->SetLabelFont(lhcbFont,"x"); +lhcbStyle->SetLabelFont(lhcbFont,"y"); +lhcbStyle->SetLabelFont(lhcbFont,"z"); +lhcbStyle->SetLabelSize(0.05,"x"); +lhcbStyle->SetLabelSize(0.05,"y"); +lhcbStyle->SetLabelSize(0.05,"z"); +lhcbStyle->SetTitleFont(lhcbFont); +lhcbStyle->SetTitleSize(0.06,"x"); +lhcbStyle->SetTitleSize(0.06,"y"); +lhcbStyle->SetTitleSize(0.06,"z"); + +// use bold lines and markers +lhcbStyle->SetLineWidth(lhcbWidth); +lhcbStyle->SetFrameLineWidth(lhcbWidth); +lhcbStyle->SetHistLineWidth(lhcbWidth); +lhcbStyle->SetFuncWidth(lhcbWidth); +lhcbStyle->SetGridWidth(lhcbWidth); +lhcbStyle->SetLineStyleString(2,"[12 12]"); // postscript dashes +lhcbStyle->SetMarkerStyle(20); +lhcbStyle->SetMarkerSize(1.5); + +// label offsets +lhcbStyle->SetLabelOffset(0.015); + +// by default, do not display histogram decorations: +lhcbStyle->SetOptStat(0); +lhcbStyle->SetOptStat("emr"); // show only nent -e , mean - m , rms -r +// full opts at http://root.cern.ch/root/html/TStyle.html#TStyle:SetOptStat +lhcbStyle->SetStatFormat("6.3g"); // specified as c printf options +lhcbStyle->SetOptTitle(0); +lhcbStyle->SetOptFit(0); +//lhcbStyle->SetOptFit(1011); // order is probability, Chi2, errors, parameters + +// look of the statistics box: +lhcbStyle->SetStatBorderSize(0); +lhcbStyle->SetStatFont(lhcbFont); +lhcbStyle->SetStatFontSize(0.05); +lhcbStyle->SetStatX(0.9); +lhcbStyle->SetStatY(0.9); +lhcbStyle->SetStatW(0.25); +lhcbStyle->SetStatH(0.15); +// put tick marks on top and RHS of plots +lhcbStyle->SetPadTickX(1); +lhcbStyle->SetPadTickY(1); + +// histogram divisions: only 5 in x to avoid label overlaps +lhcbStyle->SetNdivisions(505,"x"); +lhcbStyle->SetNdivisions(510,"y"); + + +//define style for text +TText *lhcbLabel = new TText(); +lhcbLabel->SetTextFont(lhcbFont); +lhcbLabel->SetTextColor(1); +lhcbLabel->SetTextSize(0.04); +lhcbLabel->SetTextAlign(12); + +// define style of latex text +TLatex *lhcbLatex = new TLatex(); +lhcbLatex->SetTextFont(lhcbFont); +lhcbLatex->SetTextColor(1); +lhcbLatex->SetTextSize(0.04); +lhcbLatex->SetTextAlign(12); + +// set this style +gROOT->SetStyle("lhcbStyle"); +gROOT->ForceStyle(); +} + +void printLHCb(char* optLR="L", char* optPrelim="Final", char* optText="") +{ +////////////////////////////////////////////////////////////////////////// +// routine to print 'LHCb', 'LHCb Preliminary' on plots +// options: optLR=L (top left) / R (top right) of plots +// optPrelim= Final (LHCb), Prelim (LHCb Preliminary), Other +// optText= text printed if 'Other' specified +//////////////////////////////////////////////////////////////////// + if (optLR=="R"){ + lhcbName = new TPaveText(0.70 - lhcbStyle->GetPadRightMargin(), + 0.75 - lhcbStyle->SetPadTopMargin(0.05), + 0.95 - lhcbStyle->GetPadRightMargin(), + 0.85 - lhcbStyle->SetPadTopMargin(0.05), + "BRNDC"); + } + else if (optLR="L"){ + lhcbName = new TPaveText(lhcbStyle->GetPadLeftMargin() + 0.05, + 0.87 - lhcbStyle->GetPadTopMargin(), + lhcbStyle->GetPadLeftMargin() + 0.30, + 0.95 - lhcbStyle->GetPadTopMargin(), + "BRNDC"); + } + else{ + cout << "printLHCb: option unknown" << optLR << endl; + } + if (optPrelim=="Final"){ + lhcbName->AddText("LHCb"); + } + else if (optPrelim=="Prelim"){ + lhcbName->AddText("#splitline{LHCb}{#scale[1.0]{Preliminary}}"); + } + else if (optPrelim=="Other"){ + lhcbName->AddText(optText); + } + else{ + cout << "printLHCb: option unknown " << optPrelim << endl; + } + + lhcbName->SetFillColor(0); + lhcbName->SetTextAlign(12); + lhcbName->SetBorderSize(0); + lhcbName->Draw(); +}