diff --git a/ntuple_nuVeto_small.py b/ntuple_nuVeto_small.py index a6c2c26..51f6645 100644 --- a/ntuple_nuVeto_small.py +++ b/ntuple_nuVeto_small.py @@ -470,127 +470,127 @@ startPartRange = [1] primaryMum = 0 # Look at the primary particle - for ip in startPartRange: - skipEvent=False - # Were there any hit in the tracking stations? - TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, None, None, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] ) - # If no, skip this event - if TrackSyst[0]==0: - skipEvent=True - #continue - # exit from loop on particles + ip = startPartRange[0] + skipEvent=False + # Were there any hit in the tracking stations? + TrackSyst[0],TrackSyst_eff[0],dummy,dummy,TrackSyst_Ethr[0], dummy = wasFired(None, strawPoints, trackStationsPos, None, None, checkOn=False, pointVects=None)#[strawPoint_x, strawPoint_y, strawPoint_z] ) + # If no, skip this event + if TrackSyst[0]==0: + skipEvent=True + continue + # exit from loop on particles + #break + # Select the primary MC particle + part = particles[ip] + pdgPart = pdg.GetParticle(part.GetPdgCode()) + if sampleType=="nuBg": + assert("nu" in pdgPart.GetName()) + ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother + #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena + ## (and de-indented the following block) + ## Starting the counter of how many particles were produced by the interaction of this specific nu + for recoP in recoParts: + nRecoed[0] += 1 + offline.pushOfflineByParticle(t, recoP) + if nRecoed[0]>0: + recoed[0]=1 + else: + skipEvent=True + continue + #break + nNu[0]+=1 + if part.GetMotherId()==primaryMum: + assert(primaryDone==False) + primaryDone=True + isPrimary[0] = int(True) + assert(len(idPart_fromNu)==0) + assert(len(idNeutrPart_fromNu)==0) + else: + continue + # Where did this guy interact? + tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName() + #nuInteractionNode.push_back(tmpName) + tmpName = convertion(tmpName) + if not tmpName in dictNodeNames: + tmpName = "others" + nuIntNumSimpl[0] = dictNodeNames[tmpName] + # Store interaction point coordinates + startZ_nu[0] = part.GetStartZ() + startY_nu[0] = part.GetStartY() + startX_nu[0] = part.GetStartX() + # Primary particle information + nuE[0] = part.GetEnergy() + # Looking for particles produced by the neutrino interaction + interacted = False + partKidsId = [] + NC[0] = 0 + # Add information for offline selection + # (mainly signal normalisation and zeroing of particle-wise arrays) + offline.pushOfflineByEvent(t, vetoScint, sampleType, verbose, threshold) + # Loop on MCTracks + for ip2 in xrange(0,len(particles)): + part2 = particles[ip2] + # exit if we have reached the empty part of the array + if not (type(part2)==type(ShipMCTrack())): break - # Select the primary MC particle - part = particles[ip] - pdgPart = pdg.GetParticle(part.GetPdgCode()) - if sampleType=="nuBg": - assert("nu" in pdgPart.GetName()) - # Add information for offline selection - # (mainly signal normalisation and zeroing of particle-wise arrays) - offline.pushOfflineByEvent(t, sampleType, verbose, threshold) - ## Looking for a neutrino: it should have the correct pdg code and it should not have a mother - #if (("nu" in pdgPart.GetName())):# and part.GetMotherId()==-1): # commented out by elena - ## (and de-indented the following block) - ## Starting the counter of how many particles were produced by the interaction of this specific nu - for recoP in recoParts: - nRecoed[0] += 1 - offline.pushOfflineByParticle(t, recoP) - if nRecoed[0]>0: - recoed[0]=1 - else: - skipEvent=True - #continue - break - nNu[0]+=1 - if part.GetMotherId()==primaryMum: - assert(primaryDone==False) - primaryDone=True - isPrimary[0] = int(True) - assert(len(idPart_fromNu)==0) - assert(len(idNeutrPart_fromNu)==0) - else: - continue - # Where did this guy interact? - tmpName = fGeo.FindNode(part.GetStartX(),part.GetStartY(),part.GetStartZ()).GetVolume().GetName() - #nuInteractionNode.push_back(tmpName) - tmpName = convertion(tmpName) - if not tmpName in dictNodeNames: - tmpName = "others" - nuIntNumSimpl[0] = dictNodeNames[tmpName] - # Store interaction point coordinates - startZ_nu[0] = part.GetStartZ() - startY_nu[0] = part.GetStartY() - startX_nu[0] = part.GetStartX() - # Primary particle information - nuE[0] = part.GetEnergy() - # Looking for particles produced by the neutrino interaction - interacted = False - partKidsId = [] - NC[0] = 0 - # Loop on MCTracks - for ip2 in xrange(0,len(particles)): - part2 = particles[ip2] - # exit if we have reached the empty part of the array - if not (type(part2)==type(ShipMCTrack())): - break - if part2.GetMotherId()==ip: - interacted = True - part2Id = part2.GetPdgCode() - partKidsId.append(part2Id) - # Was this a neutral current interaction? - if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0): - NC[0] = 1 - nPart_fromNu[0]+=1 - idPart_fromNu.push_back(part2Id) - if pdg.GetParticle(part2.GetPdgCode()) == None: - part2Name = str(part2Id) - else: - part2Name = getPartName(part2Id) - #idStrPart_fromNu.push_back(part2Name) - # Counting how many particles produced by the nu-interaction are charged or neutral - if (pdg.GetParticle(part2.GetPdgCode())) == None or (int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0)): - nNeutrPart_fromNu[0]+=1 - idNeutrPart_fromNu.push_back(part2Id) - #idStrNeutrPart_fromNu.push_back(part2Name) - else: - nChargedPart_fromNu[0]+=1 - idChargedPart_fromNu.push_back(part2Id) - #idStrChargedPart_fromNu.push_back(part2Name) - # How many are produced by the interaction in the last passive element of the "opera-mu system"? - # Finding out the "interaction element code" - part2Z = part2.GetStartZ() - part2X = part2.GetStartX() - part2Y = part2.GetStartY() - somewhere = False - nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName() - if nodeName in lastPassive_nodeName: - somewhere = True - interactionElement[0] = 0 - intElName = 'OPERA' - # To know if it was in the full OPERA-system (excluded last passive) - if tmpName in OPERA_nodeName and somewhere==False: - somewhere = True - interactionElement[0] = 3 - intElName = 'OPERA' - # To know if it was between the two windows - if tmpName in entrance_nodeName: - assert(somewhere==False) - somewhere = True - interactionElement[0] = 1 - # Vacuum tank outer window - if tmpName in volumeIn_nodeName and somewhere==False: - interactionElement[0] = 5 - somewhere = True - # Vacuum tank inner window - elif nodeName in volumeOut_nodeName and somewhere==False: - interactionElement[0] = 6 - somewhere = True - # Vacuum tank ribs - if "Rib" in tmpName: - interactionElement[0] = 7 - # Somewhere else - if somewhere==False: - interactionElement[0] = -1 + if part2.GetMotherId()==ip: + interacted = True + part2Id = part2.GetPdgCode() + partKidsId.append(part2Id) + # Was this a neutral current interaction? + if not (pdg.GetParticle(part2.GetPdgCode()) == None) and ("nu" in pdg.GetParticle(part2.GetPdgCode()).GetName()) and (part2.GetMotherId()==0): + NC[0] = 1 + nPart_fromNu[0]+=1 + idPart_fromNu.push_back(part2Id) + if pdg.GetParticle(part2.GetPdgCode()) == None: + part2Name = str(part2Id) + else: + part2Name = getPartName(part2Id) + #idStrPart_fromNu.push_back(part2Name) + # Counting how many particles produced by the nu-interaction are charged or neutral + if (pdg.GetParticle(part2.GetPdgCode())) == None or (int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0)): + nNeutrPart_fromNu[0]+=1 + idNeutrPart_fromNu.push_back(part2Id) + #idStrNeutrPart_fromNu.push_back(part2Name) + else: + nChargedPart_fromNu[0]+=1 + idChargedPart_fromNu.push_back(part2Id) + #idStrChargedPart_fromNu.push_back(part2Name) + # How many are produced by the interaction in the last passive element of the "opera-mu system"? + # Finding out the "interaction element code" + part2Z = part2.GetStartZ() + part2X = part2.GetStartX() + part2Y = part2.GetStartY() + somewhere = False + nodeName = fGeo.FindNode(part2X,part2Y,part2Z).GetName() + if nodeName in lastPassive_nodeName: + somewhere = True + interactionElement[0] = 0 + intElName = 'OPERA' + # To know if it was in the full OPERA-system (excluded last passive) + if tmpName in OPERA_nodeName and somewhere==False: + somewhere = True + interactionElement[0] = 3 + intElName = 'OPERA' + # To know if it was between the two windows + if tmpName in entrance_nodeName: + assert(somewhere==False) + somewhere = True + interactionElement[0] = 1 + # Vacuum tank outer window + if tmpName in volumeIn_nodeName and somewhere==False: + interactionElement[0] = 5 + somewhere = True + # Vacuum tank inner window + elif nodeName in volumeOut_nodeName and somewhere==False: + interactionElement[0] = 6 + somewhere = True + # Vacuum tank ribs + if "Rib" in tmpName: + interactionElement[0] = 7 + # Somewhere else + if somewhere==False: + interactionElement[0] = -1 # End loop on MCTracks # Should be changed to avoid entering in the loop if there isn't anything in the tracking if skipEvent: @@ -617,4 +617,4 @@ nf.Close() f.Close() print "Output wrote to %s"%outputFileName -print "Number of events with vertex outside Veto_5+500 - Tr4_4: %s"%offline.num_bad_z +print "Number of events with vertex outside Veto_5-500 - Tr4_4: %s"%offline.num_bad_z diff --git a/offlineForBarbara.py b/offlineForBarbara.py index b395ffd..d332e0a 100644 --- a/offlineForBarbara.py +++ b/offlineForBarbara.py @@ -419,7 +419,7 @@ tools.PutToZero(DaughtersFitConverged) -def pushOfflineByEvent(tree, datatype, verbose, threshold): +def pushOfflineByEvent(tree, vetoPoints, datatype, verbose, threshold): # True HNL decay vertex (only for signal normalisation) signalNormalisationZ(tree, datatype, verbose) ## Number of particles @@ -428,7 +428,7 @@ prepareFillingsByParticle() # Liquid scintillator segments global nodes - tools.hitSegments(tree, nodes, threshold) + tools.hitSegments(vetoPoints, nodes, threshold) def pushOfflineByParticle(tree, particle): hasEcalAndMuons(tree, particle) diff --git a/tools.py b/tools.py index c918adc..b1b25c0 100755 --- a/tools.py +++ b/tools.py @@ -276,12 +276,12 @@ zin, zout = None, None zstepsize = None # To hold the energy loss -listOfSegments = [] +EnergyDeposits = [] abig = None asmall = None b = None -def makeLSsegments(nodes, tree): +def makeLSsegments(nodes, tree, savelist=False): global zstarts, zends, lsNames, Nz, Nphi, nBins, zin, zout, zstepsize global abig, asmall, b abig = 250.#nodes['T2LiSc_1_1']['x']['dim'] @@ -307,19 +307,24 @@ Nz = int(math.ceil((zout-zin)/zstepsize)) #print Nphi, Nz, nBins nBins = Nphi*Nz - global listOfHitLSSegments, numberOfHitLSSegments, listOfSegments + # Init energy deposit to 0. for i in xrange(nBins): - listOfSegments.append(0.) - tree, listOfHitLSSegments = AddVect(tree, 'listOfHitLSSegments', 'int') + EnergyDeposits.append(0.) + global listOfHitLSSegments, numberOfHitLSSegments, EnergyDeposits + if savelist: + tree, listOfHitLSSegments = AddVect(tree, 'listOfHitLSSegments', 'int') tree, numberOfHitLSSegments = AddVar(tree, 'numberOfHitLSSegments', 'int') -def hitSegments(event, nodes, threshold): - global listOfHitLSSegments, numberOfHitLSSegments, listOfSegments +def hitSegments(vetoPoints, nodes, threshold, savelist=False): + global listOfHitLSSegments, numberOfHitLSSegments, EnergyDeposits global abig, asmall, b, phistart, dphi, Nphi, nBins, zin, zout + # Init energy deposit to 0. + for i in xrange(nBins): + EnergyDeposits[i]=0. 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: + #print 'nBins: ', nBins, 'len(list): ',len(EnergyDeposits) + for hit in vetoPoints: z = hit.GetZ() if (z<(nodes['VetoTimeDet_1']['z']['pos']+nodes['VetoTimeDet_1']['z']['dim']) or z > zout): @@ -334,16 +339,17 @@ phibin = ROOT.TMath.FloorNint((phi+phistart)/dphi)%Nphi iBin = zbin*Nphi+phibin try: - listOfSegments[iBin] += hit.GetEnergyLoss() + EnergyDeposits[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] > threshold: + if EnergyDeposits[iBin] > threshold: nHit +=1 - Push(listOfHitLSSegments, int(iBin)) + if savelist: + Push(listOfHitLSSegments, int(iBin)) Push(numberOfHitLSSegments, nHit)