diff --git a/offlineForBarbara.py b/offlineForBarbara.py index 4415784..53c464b 100644 --- a/offlineForBarbara.py +++ b/offlineForBarbara.py @@ -261,25 +261,27 @@ MaxDaughtersRedChi2, MinDaughtersNdf = None, None def addOfflineToTree(elenaTree): - global DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal - global DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 - global MaxDaughtersRedChi2, MinDaughtersNdf, HNLw, NuWeight - elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') # - elenaTree, DOCA = tools.AddVect(elenaTree, 'DOCA', 'float') # - elenaTree, vtxx = tools.AddVect(elenaTree, 'vtxx', 'float') # - elenaTree, vtxy = tools.AddVect(elenaTree, 'vtxy', 'float') # - elenaTree, vtxz = tools.AddVect(elenaTree, 'vtxz', 'float') # - elenaTree, IP0 = tools.AddVect(elenaTree, 'IP0', 'float') # - #elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int') # - elenaTree, DaughtersAlwaysIn = tools.AddVect(elenaTree, 'DaughtersAlwaysIn', 'int') # - elenaTree, BadTruthVtx = tools.AddVect(elenaTree, 'BadTruthVtx', 'int') # - elenaTree, Has1Muon1 = tools.AddVect(elenaTree, 'Has1Muon1', 'int') # - elenaTree, Has1Muon2 = tools.AddVect(elenaTree, 'Has1Muon2', 'int') # - elenaTree, Has2Muon1 = tools.AddVect(elenaTree, 'Has2Muon1', 'int') # - elenaTree, Has2Muon2 = tools.AddVect(elenaTree, 'Has2Muon2', 'int') # - elenaTree, HasEcal = tools.AddVect(elenaTree, 'HasEcal', 'int') # - elenaTree, MaxDaughtersRedChi2 = tools.AddVect(elenaTree, 'MaxDaughtersRedChi2', 'float') # - elenaTree, MinDaughtersNdf = tools.AddVect(elenaTree, 'MinDaughtersNdf', 'int') # + global DaughtersFitConverged, DOCA, vtxx, vtxy, vtxz, IP0, HasEcal + global DaughtersAlwaysIn, BadTruthVtx, Has1Muon1, Has1Muon2, Has2Muon1, Has2Muon2 + global MaxDaughtersRedChi2, MinDaughtersNdf, HNLw, NuWeight + elenaTree, DaughtersFitConverged = tools.AddVect(elenaTree, 'DaughtersFitConverged', 'int') # + elenaTree, DOCA = tools.AddVect(elenaTree, 'DOCA', 'float') # + elenaTree, vtxx = tools.AddVect(elenaTree, 'vtxx', 'float') # + elenaTree, vtxy = tools.AddVect(elenaTree, 'vtxy', 'float') # + elenaTree, vtxz = tools.AddVect(elenaTree, 'vtxz', 'float') # + elenaTree, IP0 = tools.AddVect(elenaTree, 'IP0', 'float') # + #elenaTree, NParticles = tools.AddVar(elenaTree, 'NParticles', 'int') # + elenaTree, DaughtersAlwaysIn = tools.AddVect(elenaTree, 'DaughtersAlwaysIn', 'int') # + elenaTree, BadTruthVtx = tools.AddVect(elenaTree, 'BadTruthVtx', 'int') # + elenaTree, Has1Muon1 = tools.AddVect(elenaTree, 'Has1Muon1', 'int') # + elenaTree, Has1Muon2 = tools.AddVect(elenaTree, 'Has1Muon2', 'int') # + elenaTree, Has2Muon1 = tools.AddVect(elenaTree, 'Has2Muon1', 'int') # + elenaTree, Has2Muon2 = tools.AddVect(elenaTree, 'Has2Muon2', 'int') # + elenaTree, HasEcal = tools.AddVect(elenaTree, 'HasEcal', 'int') # + elenaTree, MaxDaughtersRedChi2 = tools.AddVect(elenaTree, 'MaxDaughtersRedChi2', 'float') # + elenaTree, MinDaughtersNdf = tools.AddVect(elenaTree, 'MinDaughtersNdf', 'int') # + # Add liquid scintillator segmentation + tools.makeLSsegments(nodes, elenaTree) nodes = None def loadNodes(fGeo): @@ -418,12 +420,15 @@ def pushOfflineByEvent(tree, datatype, verbose): - # True HNL decay vertex (only for signal normalisation) - signalNormalisationZ(tree, datatype, verbose) - ## Number of particles - #nParticles(tree) - # Empties arrays filled by particle - prepareFillingsByParticle() + # True HNL decay vertex (only for signal normalisation) + signalNormalisationZ(tree, datatype, verbose) + ## Number of particles + #nParticles(tree) + # Empties arrays filled by particle + prepareFillingsByParticle() + # Liquid scintillator segments + global nodes + tools.hitSegments(tree, nodes) def pushOfflineByParticle(tree, particle): hasEcalAndMuons(tree, particle) diff --git a/tools.py b/tools.py index 80c66b6..d250aba 100755 --- a/tools.py +++ b/tools.py @@ -269,20 +269,25 @@ zstarts = [] zends = [] lsNames = [] -nZ, Nphi, nBins = None, None, None +Nz, Nphi, nBins = None, None, None Nphi = 16 dphi = 2.*math.pi/Nphi phistart = dphi/2. -nin, zout = None, None +zin, zout = None, None zstepsize = None # To hold the energy loss listOfSegments = [] -abig = 250. -b = 500. -asmall = 250. +abig = None +asmall = None +b = None def makeLSsegments(nodes, tree): global zstarts, zends, lsNames, Nz, Nphi, nBins, zin, zout, zstepsize + global abig, asmall, b + abig = 250.#nodes['T2LiSc_1_1']['x']['dim'] + asmall = 186.#nodes['T1LiSc_1_1']['x']['dim'] + b = 500.#nodes['T2LiSc_1_1']['y']['dim'] + print abig, asmall, b for node in nodes.keys(): if 'LiSc' in node: zstarts.append(nodes[node]['z']['pos']-nodes[node]['z']['dim']) @@ -294,39 +299,51 @@ # It cannot be so bad... zstarts.sort() zends.sort() - Nz = len(zstarts) + #Nz = len(zstarts) zin = min(zstarts) - zout = max(zends) + zout =max(zends) + print 'zstart: ', zin, 'zend: ', zout zstepsize = 100. - #Nz = int(math.ceil((zi-zo)/zstepsize)) + Nz = int(math.ceil((zout-zin)/zstepsize)) + print Nphi, Nz, nBins nBins = Nphi*Nz - global listOfHitSegments, numberOfHitLSSegments, listOfSegments + global listOfHitLSSegments, numberOfHitLSSegments, listOfSegments for i in xrange(nBins): listOfSegments.append(0.) tree, listOfHitLSSegments = AddVect(tree, 'listOfHitLSSegments', 'int') tree, numberOfHitLSSegments = AddVar(tree, 'numberOfHitLSSegments', 'int') def hitSegments(event, nodes): - global listOfHitSegments, numberOfHitLSSegments, listOfSegments - global abig, asmall, b, phistart, dphi, Nphi, nBins + global listOfHitLSSegments, numberOfHitLSSegments, listOfSegments + global abig, asmall, b, phistart, dphi, Nphi, nBins, zin, zout endsmall = nodes['T1Endplate_1']['z']['pos']+nodes['T1Endplate_1']['z']['dim'] + print 'endsmall: ', endsmall + print 'nBins: ', nBins, 'len(list): ',len(listOfSegments) for hit in event.vetoPoint: z = hit.GetZ() + if (z<(nodes['VetoTimeDet_1']['z']['pos']+nodes['VetoTimeDet_1']['z']['dim']) + or z > zout): + continue x = hit.GetX() y = hit.GetY() a = abig if z < endsmall: a = asmall phi = ROOT.TMath.ATan2(a*y,b*x)+2*math.pi - zbin = int(math.floor( (z-zin)/zstepsize )) - phibin = int((math.floor( (phi+phistart)/dphi ))%Nphi) + zbin = ROOT.TMath.FloorNint((z-zin)/zstepsize) + phibin = ROOT.TMath.FloorNint((phi+phistart)/dphi)%Nphi iBin = zbin*Nphi+phibin - listOfSegments[iBin] += hit.GetEnergyLoss() + try: + listOfSegments[iBin] += hit.GetEnergyLoss() + except: + print 'hitSegments: could not find bin ', iBin + print x, y, z + sys.exit() nHit = 0 for iBin in xrange(nBins): if listOfSegments[iBin] > 0.045: nHit +=1 - Push(listOfHitSegments, iBin) + Push(listOfHitLSSegments, int(iBin)) Push(numberOfHitLSSegments, nHit)