diff --git a/distribsForHNLandBG.py b/distribsForHNLandBG.py index 9c63993..4433df3 100644 --- a/distribsForHNLandBG.py +++ b/distribsForHNLandBG.py @@ -3,7 +3,7 @@ import funcsByBarbara as tools r.gROOT.SetBatch(r.kTRUE) -r.gROOT.ProcessLine(".L lhcbStyle.C+") +r.gROOT.ProcessLine(".X lhcbStyle.C") plots = {} def addPlot(key, x_axis, tree, string, cuts=None): @@ -19,22 +19,56 @@ plots[key].Draw() c1.Print('accPlots/%s.pdf'%key) c1.Print('accPlots/%s.png'%key) + key += "-WEIGHTED" + if 'bg' in key: + if not cuts: + tree.Draw(string+'>>h1', "NuWeight") + else: + cuts = "NuWeight*("+cuts+")" + 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) + else: + if not cuts: + tree.Draw(string+'>>h1',"HNLw") + else: + cuts = "HNLw*("+cuts+")" + 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') +#pimu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/old_samples/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root') +pimu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/MUPI/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-2985556.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') +#mumunu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/old_samples/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008.root') +mumunu_geo = tools.searchForNodes2_xyz_dict('../Geraldine/MUMUNU/geofile_full.10.0.Pythia8-TGeant4-1.0GeV-008-TOTAL.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') +#bg_geo = tools.searchForNodes2_xyz_dict('../Geraldine/KLONG/geofile_full.10.0.Genie-TGeant4.root') +bg_geo = tools.searchForNodes2_xyz_dict('../Geraldine/KLONG/geofile_full.10.0.Genie-TGeant4_9files.root') studies = { - 'pimu': {'data':pimu, 'geo':pimu_geo}, - 'mumunu': {'data':mumunu, 'geo':mumunu_geo}, - 'bg': {'data':bg, 'geo':bg_geo} + 'pimu': {'data':pimu, 'geo':pimu_geo, 'ntot':1000}, + 'mumunu': {'data':mumunu, 'geo':mumunu_geo, 'ntot':7000}, + 'bg': {'data':bg, 'geo':bg_geo, 'ntot':99999} } # Cuts: @@ -67,11 +101,15 @@ 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 ip_cut(event, maxip): + if event.IP0 < maxip: + return True + return False + def nothing_in_veto5(event): if len(event.veto5_x) > 0: return False @@ -82,25 +120,57 @@ return False return True +def weight(event, study): + w = event.HNLw / studies[s]['ntot'] + if 'bg' in study: + w = event.NuWeight + return w + +results = {} for s in studies: print 'Reading %s...'%s t = studies[s]['data'] ntot = t.GetEntriesFast() - nAfterCuts = [0]*6 + nAfterCuts = [0]*10 + nAfterCuts_weighted = [0.]*10 for event in t: + w = weight(event,s) if fit_converged(event): nAfterCuts[0] += 1 + nAfterCuts_weighted[0] += w 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 + nAfterCuts_weighted[1] +=w + if chi2_cut(event,15): + nAfterCuts[2]+=1 + nAfterCuts_weighted[2]+=w + if doca_cut(event, 50.): + nAfterCuts[2+1] += 1 + nAfterCuts_weighted[2+1] += w + if doca_cut(event, 30.): + nAfterCuts[3+1]+=1 + nAfterCuts_weighted[3+1]+=w + if z_cut(event, s): + nAfterCuts[4+1] += 1 + nAfterCuts_weighted[4+1] += w + if ip_cut(event,500.): + nAfterCuts[5+1]+=1 + nAfterCuts_weighted[5+1]+=w + if ip_cut(event,250.): + nAfterCuts[6+1]+=1 + nAfterCuts_weighted[6+1]+=w + if nothing_in_veto5(event): + nAfterCuts[7+1] += 1 + nAfterCuts_weighted[7+1] += w + if nothing_in_liquidscint(event): + nAfterCuts[8+1] += 1 + nAfterCuts_weighted[8+1] += w print '%s \t survived particles: '%s, nAfterCuts, ' of %s total'%ntot + print '%s \t weighted survived particles: '%s, nAfterCuts_weighted, ' of %s total'%ntot + results[s] = {} + results[s]['ntot']=studies[s]['ntot'] + results[s]['numbers']=nAfterCuts + results[s]['numbersWithWeights']=nAfterCuts_weighted 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' ) @@ -110,7 +180,7 @@ 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" + cut += " && DOCA < 25" 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) @@ -122,3 +192,7 @@ addPlot(s+'-chi2andDOCAandZcuts-dPt', 'Tracks Pt', t, 'DaughtersPt' , cut) addPlot(s+'-chi2andDOCAandZcuts-VtxZ', 'Reco vertex z position [cm]', t, 'vtxz' , cut) + +for s in ['pimu', 'mumunu', 'bg']: + for list in results[s]: + print s,list,results[s][list] diff --git a/funcsByBarbara.py b/funcsByBarbara.py index 26f6561..88471ff 100644 --- a/funcsByBarbara.py +++ b/funcsByBarbara.py @@ -1,7 +1,29 @@ from elena import * from array import array +from collections import OrderedDict import ROOT +# Custom sorting +def sortNodesBy(nodes, coord='z', what='pos', printdict=True): + d = OrderedDict(sorted(nodes.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 + +def retrieveMCParticleInfo(sTree, key): + mcPartKey = sTree.fitTrack2MC[key] + mcPart = sTree.MCTrack[mcPartKey] + if not mcPart: + print "Error in retrieveMCParticleInfo" + sys.exit() + mother = sTree.MCTrack[mcPart.GetMotherId()] + return int(mcPart.GetPdgCode()), int(mother.GetPdgCode()), mcPart.GetStartX(), mcPart.GetStartY(), mcPart.GetStartZ() + + def AddVect(t,name,vectType): vect = ROOT.vector(vectType)() t.Branch( name, vect ) @@ -67,3 +89,85 @@ #return flag, flag_eff, nstat, nstat_eff, flag_Ethr return flag, flag_Ethr + +#ff = ROOT.TFile("histoForWeights.root","read") +#h_GioHans = ff.Get("h_Gio") + +def calcWeight(x,y,z, E, nodes, entries, weightHist, datatype): + if datatype == 'sig': + return -999 + params = {'OPERA': {'material':'Fe','lenght':60.,}, + 'window1': {'material':'Fe','lenght':nodes["lidT1O_1"]['z']['dim']}, + 'window2': {'material':'Al','lenght':nodes["lidT1I_1"]['z']['dim']}, + 'volumeOut': {'material':'Fe','lenght':30.,}, + 'volumeIn': {'material':'Al','lenght':8.}} + materialPars = {'Fe':{'A':55.847, #g/mol # fuori + 'rho': 7874 * 1.e+3/1.e+6,#g/cm3 + }, + 'Al':{'A':26.98, #g/mol # dentro + 'rho': 2700 * 1.e+3/1.e+6 #g/cm3 + }} + perc = {'OPERA':0.33, + 'window1':0.015, + 'window2':0.015, + 'volumeOut':0.32, + 'volumeIn':0.32} + intElName = interactionElement(x,y,z, nodes) + NA = 6.022e+23 + material = params[intElName]['material'] + crossSec = 8.8e-39*E #1.5e-38*fabs(E) ##3./2*(1.2+0.35)*1.e-38*fabs(E) + flux = 2.e+18#5.e+16 + L = params[intElName]['lenght'] + rho = materialPars[material]['rho'] + n_element = perc[intElName]*entries + binN = weightHist.GetXaxis().FindBin(E) + weight = crossSec * NA * rho * L * flux / n_element * weightHist.GetBinContent(binN) + return weight + +def interactionElement(x,y,z, nodes): + # window1 + if nodes["lidT1O_1"]['z']['pos']-nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim']: + return "window1" + # window2 + if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT1I_1"]['z']['pos']+nodes["lidT1I_1"]['z']['dim']: + return "window2" + # vacuum tank borders + if nodes["lidT1O_1"]['z']['pos']+nodes["lidT1O_1"]['z']['dim'] < z < nodes["lidT6O_1"]['z']['pos']+nodes["lidT6O_1"]['z']['dim']: # include le lid x la parte fuori ma non x la parte dentro + if inInnerLid(x,y,z, nodes): + return "volumeIn" + else: + return "volumeOut" + # opera + try: + minz, maxz = nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'], nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim'] + except KeyError: + try: + minz, maxz = nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'], nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim'] + except KeyError: + # Pazienza esaurita, ora lo hanno chiamato "volMagneticSpectrometer_1" ed e' lungo piu' di 9m!!! + try: + minz, maxz = nodes["volMagneticSpectrometer_1"]['z']['pos']-nodes["volMagneticSpectrometer_1"]['z']['dim'], nodes["volMagneticSpectrometer_1"]['z']['pos']+nodes["volMagneticSpectrometer_1"]['z']['dim'] + except KeyError: + print "Cannot find OPERA" + sys.exit() + #if nodes["volIron_12"]['z']['pos']-nodes["volIron_12"]['z']['dim'] < z < nodes["volIron_24"]['z']['pos']+nodes["volIron_24"]['z']['dim']: + #if nodes["volLayer_11"]['z']['pos']-nodes["volLayer_11"]['z']['dim'] < z < nodes["volLayer_0"]['z']['pos']+nodes["volLayer_0"]['z']['dim']: + if minz < z < maxz: + return "OPERA" + + +def inInnerLid(x,y,z, nodes): + if nodes["lidT1I_1"]['z']['pos']-nodes["lidT1I_1"]['z']['dim'] < z < nodes["lidT6I_1"]['z']['pos']+nodes["lidT6I_1"]['z']['dim']: + if z < nodes["Veto_5"]['z']['pos']-nodes["Veto_5"]['z']['dim']: + a,b = nodes['T1I_1']['x']['dim'], nodes['T1I_1']['y']['dim'] # smaller part of the tank (first 5 m) + else: + a,b = nodes['T2I_1']['x']['dim'], nodes['T2I_1']['y']['dim'] # larger part of the tank + if inEllipse(x,y,a,b): + return True + return False + + +def inEllipse(x,y,a,b): + if ( (x*x)/(a*a) + (y*y)/(b*b) ) < 1.: + return True + return False \ No newline at end of file diff --git a/histoForWeights.root b/histoForWeights.root new file mode 100644 index 0000000..a6d91b1 --- /dev/null +++ b/histoForWeights.root Binary files differ