diff --git a/ntuple_nuVeto_small.py b/ntuple_nuVeto_small.py index 2035142..c5dc16e 100644 --- a/ntuple_nuVeto_small.py +++ b/ntuple_nuVeto_small.py @@ -5,10 +5,12 @@ fileName, fileNameGeo, outputFileName, sampleType = None, None, None, None maxevents = 0 verbose = True +# Energy deposit threshold for the liquid scintillator +threshold = 0.015 # Parse commandline inputs try: - opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:", []) + opts, args = getopt.getopt(sys.argv[1:], "n:t:f:o:v:", ['--Ethr=']) except getopt.GetoptError: # print help information and exit: print '\tUSAGE: -f inputfile.root -o outputfile.root -t sampleType (sig, nuBg or cosmics) -n maxevents -v verbose (0 or 1)' @@ -24,6 +26,8 @@ sampleType = str(a) if o in ("-v"): verbose = bool(int(a)) + if o in ("--Ethr="): + threshold = float(a) fileNameGeo = fileName.replace('ship', 'geofile_full') @@ -35,7 +39,7 @@ print '\tFATAL! %s not defined!'%namestr(item)[0] sys.exit() -for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose]: +for item in [fileName, fileNameGeo, outputFileName, sampleType, maxevents, verbose, threshold]: print '\tINFO: %s set to %s'%(namestr(item)[0], item) # Load SHiP (LHCb) style @@ -467,13 +471,15 @@ 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 - skipEvent=False + #continue + # exit from loop on particles + break # Select the primary MC particle part = particles[ip] pdgPart = pdg.GetParticle(part.GetPdgCode()) @@ -481,117 +487,118 @@ assert("nu" in pdgPart.GetName()) # Add information for offline selection # (mainly signal normalisation and zeroing of particle-wise arrays) - offline.pushOfflineByEvent(t, sampleType, verbose) + 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 - if True: - ## 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 - 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 - # End loop on MCTracks + ## (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 + # End loop on MCTracks # Should be changed to avoid entering in the loop if there isn't anything in the tracking if skipEvent: continue strawVetoAny[0],strawVetoAny_eff[0],dummy,dummy,strawVetoAny_Ethr[0],dummy = wasFired(None, strawPoints, strawVetoPos,None,None, checkOn=False, pointVects=None)#[strawVetoPoint_x, strawVetoPoint_y, strawVetoPoint_z]) upstreamVetoany[0], upstreamVetoany_eff[0], dummy,dummy, upstreamVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, upstreamVetoPos, None, None, checkOn=False, pointVects=None)#[upstreamVetoPoint_x, upstreamVetoPoint_y, upstreamVetoPoint_z]) RPCany[0], RPCany_eff[0],dummy,dummy, RPCany_Ethr[0],dummy = wasFired(None, rpc, RPCstationsPos, None, None, checkOn=False, pointVects=None)#[rpcPoint_x, rpcPoint_y, rpcPoint_z]) - scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=0.015) + scintVetoAny[0], scintVetoAny_eff[0], dummy,dummy, scintVetoAny_Ethr[0],dummy = wasFired(None, vetoScint, vetoWall, None, None, checkOn=False, pointVects=None, Ethr=threshold) #assert(len(idStrPart_fromNu)==len(idPart_fromNu)) #assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu)) #assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu)) diff --git a/offlineForBarbara.py b/offlineForBarbara.py index 53c464b..b395ffd 100644 --- a/offlineForBarbara.py +++ b/offlineForBarbara.py @@ -419,7 +419,7 @@ tools.PutToZero(DaughtersFitConverged) -def pushOfflineByEvent(tree, datatype, verbose): +def pushOfflineByEvent(tree, 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) + tools.hitSegments(tree, nodes, threshold) def pushOfflineByParticle(tree, particle): hasEcalAndMuons(tree, particle) diff --git a/tools.py b/tools.py index 65caf87..c918adc 100755 --- a/tools.py +++ b/tools.py @@ -313,7 +313,7 @@ tree, listOfHitLSSegments = AddVect(tree, 'listOfHitLSSegments', 'int') tree, numberOfHitLSSegments = AddVar(tree, 'numberOfHitLSSegments', 'int') -def hitSegments(event, nodes): +def hitSegments(event, nodes, threshold): 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'] @@ -341,7 +341,7 @@ sys.exit() nHit = 0 for iBin in xrange(nBins): - if listOfSegments[iBin] > 0.015: + if listOfSegments[iBin] > threshold: nHit +=1 Push(listOfHitLSSegments, int(iBin)) Push(numberOfHitLSSegments, nHit)