from ROOT import *
from array import array
#gROOT.ProcessLine(".x mystyle.C")
pdg = TDatabasePDG.Instance()
gridNumber = 697
folder = ""
debug = True
## are the weights to be switched on?
ON = False
if debug:
#fileName = "../data/neutrino661/ship.10.0.Genie-TGeant4.root"
gridNumber = 691
fileName = "ship.10.0.Genie-TGeant4_n%s.root"%gridNumber
outputFileName = "nuNutple_debug%s.root"%gridNumber
else:
#fileName = "../data/neutrino661/ship.10.0.Genie-TGeant4.root"#"../data/all/ship.10.0.Genie-TGeant4-370k.root"
fileName = "ship.10.0.Genie-TGeant4_n%s.root"%(gridNumber)#"../data/all/ship.10.0.Genie-TGeant4-370k.root"
outputFileName = "nuNutpleStudy%s.root"%gridNumber
fileNameGeo = "ship.10.0.Genie-TGeant4.root"#"../data/neutrino681/ship.10.0.Genie-TGeant4.root"
print gridNumber
### Add what is missed:
### PDG nuclear states are 10-digit numbers
### 10LZZZAAAI e.g. deuteron is
### 1000010020
### from http://svn.cern.ch/guest/AliRoot/trunk/STEER/STEERBase/AliPDG.cxx and https://github.com/FairRootGroup/FairRoot/blob/master/eventdisplay/FairEventManager.cxx
### https://geant4.web.cern.ch/geant4/UserDocumentation/UsersGuides/ForApplicationDeveloper/html/AllResources/TrackingAndPhysics/particleList.src/ions/index.html
def PrintEventPart(particles,pdg):
print "Particles in the event:"
for (pid,part) in enumerate(particles):
print pid, pdg.GetParticle(part.GetPdgCode()).GetName(), part.GetMotherId()
print
def PrintRPChits(rpc,pdg):
print "Hits in the RPC:"
for (ix,RPCpt) in enumerate(rpc):
print ix,RPCpt.GetZ(), pdg.GetParticle(RPCpt.PdgCode()).GetName(), RPCpt.GetTrackID()
print
if not(pdg.GetParticle(1000010020)):
pdg.AddParticle("Deuteron","Deuteron", 1.875613e+00, kTRUE, 0,3,"Ion",1000010020)
pdg.AddAntiParticle("AntiDeuteron", - 1000010020)
if not(pdg.GetParticle(1000010030)):
pdg.AddParticle("Triton","Triton", 2.808921e+00, kFALSE, 3.885235e+17,3,"Ion",1000010030);
pdg.AddAntiParticle("AntiTriton", - 1000010030);
if not(pdg.GetParticle(1000020040) ):
#pdg.AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,khShGev/(12.33*kYear2Sec),6,"Ion",ionCode);
pdg.AddParticle("Alpha","Alpha",3.727379e+00,kFALSE,0,6,"Ion",1000020040);
pdg.AddAntiParticle("AntiAlpha", - 1000020040);
if not(pdg.GetParticle(1000020030) ):
pdg.AddParticle("HE3","HE3",2.808391e+00,kFALSE,0,6,"Ion",1000020030);
pdg.AddAntiParticle("AntiHE3", -1000020030);
if not(pdg.GetParticle(1000030070) ):
print "TO BE CHECKED the data insert for Li7Nucleus"
pdg.AddParticle("Li7Nucleus","Li7Nucleus",3.727379e+00/4.*7.,kFALSE,0,0,"Boson",1000030070);
if not(pdg.GetParticle(1000060120) ):
print "ERROR: random values insert for C12Nucleus"
pdg.AddParticle("C12Nucleus","C12Nucleus",0.1,kFALSE,0,0,"Isotope",1000060120);
if not(pdg.GetParticle(1000030060) ):
print "ERROR: random values insert for Li6Nucleus"
pdg.AddParticle("Li6Nucleus","Li6Nucleus",6.015121,kFALSE,0,0,"Isotope",1000030060);
if not(pdg.GetParticle(1000070140) ):
print "ERROR: random values insert for for N14"
pdg.AddParticle("N14","N14",0.1,kTRUE,0,0,"Isotope",1000070140);
if not(pdg.GetParticle(1000050100) ):
print "TO BE CHECKED the data insert for B10"
pdg.AddParticle("B10","B10",10.0129370,kTRUE,0,0,"Isotope",1000050100);
if not(pdg.GetParticle(1000020060) ):
print "TO BE CHECKED the data insert for He6"
pdg.AddParticle("He6","He6",6.0188891,kFALSE,806.7e-3,0,"Isotope",1000020060);
if not(pdg.GetParticle(1000040080) ):
print "TO BE CHECKED the data insert for Be8"
pdg.AddParticle("Be8","Be8",8.00530510,kFALSE,6.7e-17,0,"Isotope",1000040080);
if not(pdg.GetParticle(1000030080) ):
print "TO BE CHECKED the data insert for Li8"
pdg.AddParticle("Li8","Li8",8.002248736,kTRUE,178.3e-3,0,"Isotope",1000030080);
if not(pdg.GetParticle(1000040170) ):
print "ERROR: didn't find what is it particle with code 1000040170, random number inserted!"
pdg.AddParticle("None","None",0.1,kFALSE,0,0,"Isotope",1000040170);
if not(pdg.GetParticle(1000040100) ):
print "TO BE CHECKED the data insert for Be10"
pdg.AddParticle("Be10","Be10",10.0135338,kTRUE,5.004e+9,0,"Isotope",1000040100);
if not(pdg.GetParticle(1000040070) ):
print "TO BE CHECKED the data insert for Be7"
pdg.AddParticle("Be7","Be7",11.021658,kTRUE,13.81,0,"Isotope",1000040070);
if not(pdg.GetParticle(1000230470) ):
print "ERROR: didn't find what is it particle with code 1000230470, random number inserted!"
pdg.AddParticle("None2","None2",0.1,kFALSE,0,0,"Isotope",1000230470);
if not(pdg.GetParticle(1000080170) ):
print "ERROR: didn't find what is it particle with code 1000080170, random number inserted!"
pdg.AddParticle("None3","None3",0.1,kFALSE,0,0,"Isotope",1000080170);
if not(pdg.GetParticle(1000240500) ):
print "ERROR: didn't find what is it particle with code 1000240500, random number inserted!"
pdg.AddParticle("None4","None4",0.1,kFALSE,0,0,"Isotope",1000240500);
if not(pdg.GetParticle(1000210450) ):
print "ERROR: didn't find what is it particle with code 1000210450, random number inserted (Sc45Nucleus)!"
pdg.AddParticle("Sc45Nucleus","Sc45Nucleus",0.1,kFALSE,0,0,"Isotope",1000210450);
if not(pdg.GetParticle(1000040090) ):
print "TO BE CHECKED the data insert for Be9"
pdg.AddParticle("Be9","Be9",9.0121822,kFALSE,0,0,"Isotope",1000040090);
if not(pdg.GetParticle(1000080160) ):
print "TO BE CHECKED the data insert for O16"
pdg.AddParticle("O16","O16",15.99491461956,kFALSE,0,0,"Isotope",1000080160);
if not(pdg.GetParticle(1000220460) ):
print "TO BE CHECKED the data insert for Ar40"
pdg.AddParticle("Ar40","Ar40",39.9623831225,kFALSE,0,0,"Isotope",1000220460);
if not(pdg.GetParticle(1000050110) ):
print "TO BE CHECKED the data insert for B11"
pdg.AddParticle("B11","B11",11.0093054,kFALSE,0,0,"Isotope",1000050110);
KsID = 310
KLID = 130
NID = 2112
# E = neutrino energy
# ON flag to switch on or off the weights
def calcWeight(intElName, E, params, entries, ON=True):
if not ON:
return 1
materialPars = {'Fe':{'A':55.847, #g/mol
'rho': 7874 * 1.e+3/1.e+6,#g/cm3
},
'Al':{'A':26.98, #g/mol
'rho': 2700 * 1.e+3/1.e+6 #g/cm3
}}
perc = {'OPERA':0.33,
'window1':0.015,
'window2':0.015,
'volume1':0.32,
'volume2':0.32}
NA = 6.022e+23
material = params[intElName]['material']
# from PDG CC_crossSec = 1.2 for nu, and 0.35 for anti-nu, then we guess-estimate that the NN are 0.5 of CC -> 3/2.
crossSec = 3./2*(1.2+0.35)*1.e-38*fabs(E)
flux = 5.e+16
L = params[intElName]['lenght']
rho = materialPars[material]['rho']
n_element = perc[intElName]*entries
if intElName=="OPERA":
print crossSec
print L
print rho
print E
print "Weight: ", crossSec * NA * rho * flux * L/ n_element
assert(False)
return crossSec * NA * rho * L * flux / n_element
def addVect(t,name,vectType):
vect = vector(vectType)()
t.Branch( name, vect )
return t, vect
def addVar(t,name,varType):
var = array(varType,[-999])
t.Branch(name,var,name+"/"+varType.upper())
return t, var
def putToZero(var):
#var.clear()
#var.push_back(0)
var[0] = 0
def getPartName(partId):
return pdg.GetParticle(partId).GetName()
## pointVects is a list of pointVect [x,y,z]. The pointVect are the vector
## to be added to the tree containing the position of the hit in the detector
## system under consideration. They are put for the case of 100% efficiency.
def wasFired(partKidsId, detPoints, detPos, checkOn, pointVects=None):
def lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,partKidTrackID, t_ID):
for hit in detPoints:
if (hit.GetTrackID()==t_ID) or (partKidTrackID is None):
#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]:
#print pos, hit.GetZ()
if pointVects is not None:
pointVects[0].push_back(hit.GetX())
pointVects[1].push_back(hit.GetY())
pointVects[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
return flag, flag_eff, nstat, nstat_eff
def lookingForTrackID(partKidsId,detPoints):
partKidTrackID = []
for p_ID in partKidsId:
for hit in detPoints:
if p_ID==hit.PdgCode() and not (hit.GetTrackID() in partKidTrackID):
# I take the trackID to follow it
partKidTrackID.append(hit.GetTrackID())
### to be conservative I take only the first trackID of the first set of hits corrisponding to the particleID (nu_daughter)
break
return partKidTrackID
# Now in partKidTrackID I should have the trackID connected to my charged particles
flag_eff = False
flag = False
nstat = 0
nstat_eff = 0
# I should find the first hit connected to the particles from the interaction
if not (partKidsId is None):
partKidTrackID = lookingForTrackID(partKidsId,detPoints)
for t_ID in partKidTrackID:
flag,flag_eff,nstat,nstat_eff = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,partKidTrackID, t_ID)
else:
## in case we want to look at the presence of any hit, no matter what particle left it
flag,flag_eff,nstat,nstat_eff = lookingForHits(detPoints,flag,flag_eff,nstat,nstat_eff,None, None)
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
# RPC[0] = int(rpcFlag)
# RPC_eff[0] = int(RPCflag_eff)
# nRPCstations[0] = int(nRPCstat)
# nRPCstations_eff[0] = int(nRPCstat_eff)
f = TFile(fileName)
t = f.Get("cbmsim")
entries = t.GetEntries()
#t2.Draw("startZ_nu","startZ_nu>-2500.8 && startZ_nu<-2400.8")
from lookAtGeo import *
myNodes_name = ["volLayer_%s"%i for i in xrange(0,12)]
myNodes_name += ["lidT1lisci_1","lidT1I_1","lidT1O_1"]
myNodes_name += ["volScintLayer_%s"%i for i in xrange(0,11)]
myNodes_name += ["lidT6lisci_1","lidT6I_1","lidT6O_1"]
myNodes_name += ["volDriftLayer%s_1"%i for i in xrange(1,6)]
myNodes_name += ["Tr%s_%s"%(i,i) for i in xrange(1,5)]
myNodes_name += ["Veto_5"]
myGeoEl = findPositionGeoElement(fileNameGeo, myNodes_name)
lastPassiveEl = [myGeoEl["volLayer_0"]['z']-myGeoEl["volLayer_0"]['dimZ'],myGeoEl["volLayer_0"]['z']+myGeoEl["volLayer_0"]['dimZ']]
geo = loadGeometry(fileNameGeo)
ship_geo = geo['ShipGeo']
entranceWindows = [ [myGeoEl["lidT1O_1"]['z']-myGeoEl["lidT1O_1"]['dimZ'],myGeoEl["lidT1O_1"]['z']+myGeoEl["lidT1O_1"]['dimZ']],
[myGeoEl["lidT1I_1"]['z']-myGeoEl["lidT1I_1"]['dimZ'],myGeoEl["lidT1I_1"]['z']+myGeoEl["lidT1I_1"]['dimZ']]
]
volume = [myGeoEl["lidT1O_1"]['z']-myGeoEl["lidT6O_1"]['dimZ'],myGeoEl["lidT6O_1"]['z']-myGeoEl["lidT6O_1"]['dimZ']]#[myGeoEl["lidT1lisci_1"]['z']+myGeoEl["lidT1lisci_1"]['dimZ'],4100]
scintTankW = [myGeoEl["lidT1lisci_1"]['z']-myGeoEl["lidT1lisci_1"]['dimZ'],myGeoEl["lidT1lisci_1"]['z']+myGeoEl["lidT1lisci_1"]['dimZ']]
scintTankV = [myGeoEl["lidT1lisci_1"]['z']-myGeoEl["lidT1lisci_1"]['dimZ'],myGeoEl["lidT1lisci_1"]['z']+myGeoEl["lidT1lisci_1"]['dimZ']]
OPERA = [myGeoEl["volLayer_11"]['z']-myGeoEl["volLayer_11"]['dimZ'],myGeoEl["volLayer_0"]['z']+myGeoEl["volLayer_0"]['dimZ']]
OPERA_wrong = [-2705.3,-2625.3]
Tracking = [myGeoEl["Tr1_1"]['z']-myGeoEl["Tr1_1"]['dimZ'],myGeoEl["Tr4_4"]['z']+myGeoEl["Tr4_4"]['dimZ']]
print "Geometry Parameters:"
print "last passive (interactionElement == 0): ",lastPassiveEl
print "entrance windows (interactionElement == 1): ", entranceWindows
print "volume (interactionElement ==2): ", volume
print "scintTankW: ", scintTankW
print "scintTankV: ", scintTankV
print "OPERA-system: ",OPERA
print "OPERA_wrong: ", OPERA_wrong
print "Tracking: ", Tracking
#lastPassive = [myGeoEl["volLayer2_11"]['z'],myGeoEl["volLayer2_0"]['z']]
print "I should think of a correct way to distinguish volume from entrance window around the z of the two windows"
#assert(False)
## I don't want to consider the RPC after the last passive element (removed in some later version of the geo) so index just up to 11
RPCstationsPos = [[myGeoEl["volScintLayer_%s"%i]['z']-myGeoEl["volScintLayer_%s"%i]['dimZ'],myGeoEl["volScintLayer_%s"%i]['z']+myGeoEl["volScintLayer_%s"%i]['dimZ']] for i in xrange(0,11)]
HPTstationsPos = [[myGeoEl["volDriftLayer%s_1"%i]['z']-myGeoEl["volDriftLayer%s_1"%i]['dimZ'],myGeoEl["volDriftLayer%s_1"%i]['z']+myGeoEl["volDriftLayer%s_1"%i]['dimZ']] for i in xrange(1,6)]
n_RPC = len(RPCstationsPos)
RPCstationsPos += HPTstationsPos
trackStationsPos = [[myGeoEl["Tr%s_%s"%(i,i)]['z']-myGeoEl["Tr%s_%s"%(i,i)]['dimZ'],myGeoEl["Tr%s_%s"%(i,i)]['z']+myGeoEl["Tr%s_%s"%(i,i)]['dimZ']] for i in xrange(1,5)]
strawVetoPos = [[myGeoEl["Veto_5"]['z']-myGeoEl["Veto_5"]['dimZ'],myGeoEl["Veto_5"]['z']+myGeoEl["Veto_5"]['dimZ']]]
scintWPos = [scintTankW]
params = {'OPERA': {'material':'Fe',
'lenght':OPERA[1]-OPERA[0]-n_RPC*myGeoEl["volScintLayer_0"]['dimZ']*2.,},## I want to calculate only the passive part, remember that the dimZ is the half-width.
'window1': {'material':'Fe',
'lenght':entranceWindows[0][1]-entranceWindows[0][0]},
'window2': {'material':'Al',
'lenght':entranceWindows[1][1]-entranceWindows[1][0],},
'volume1': {'material':'Fe',
'lenght':volume[1]-volume[0],},
'volume2': {'material':'Al',
'lenght':volume[1]-volume[0]}}
### Ntuple elements:
nf = TFile(outputFileName,"RECREATE")
nt = TTree("t","t")
nt, event = addVar(nt, 'event','i')
nt, nNu = addVar(nt, 'nNu', 'i')
nt, isPrimary = addVar(nt, 'isPrimary','i')
nt, startZ_nu = addVar(nt, 'startZ_nu', 'f')
nt, startY_nu = addVar(nt, 'startY_nu', 'f')
nt, startX_nu = addVar(nt, 'startX_nu', 'f')
nt, nuE = addVar(nt, 'nuE', 'f')
nt, nuP = addVar(nt, 'nuP', 'f')
## 0: last passive material opera-mu-system
## 1: two windows around liquid scintillator
## 2: along vaccum tank
## 3: along all OPERA
## 4: between the two entrance windows
## 999: wrong OPERA place ## needed until we have the new production
## -1: anything else, but what???? #it should never been present
nt, interactionElement = addVar(nt, 'interactionElement','i')
### From Neutrino
nt, nPart_fromNu = addVar(nt, 'nPart_fromNu', 'i')
nt, idPart_fromNu = addVect(nt, 'idPart_fromNu', 'int')
nt, idStrPart_fromNu = addVect(nt, 'idStrPart_fromNu', 'string')
nt, nChargedPart_fromNu = addVar(nt, 'nChargedPart_fromNu', 'i')
nt, idChargedPart_fromNu = addVect(nt, 'idChargedPart_fromNu', 'int')
nt, idStrChargedPart_fromNu = addVect(nt, 'idStrChargedPart_fromNu', 'string')
nt, nNeutrPart_fromNu = addVar(nt, 'nNeutrPart_fromNu', 'i')
nt, idNeutrPart_fromNu = addVect(nt, 'idNeutrPart_fromNu', 'int')
nt, idStrNeutrPart_fromNu = addVect(nt, 'idStrNeutrPart_fromNu', 'string')
### With Ks
nt, nNeutrPart_withKs = addVar(nt,'nNeutrPart_withKs', 'i')
nt, idNeutrPart_withKs = addVect(nt, 'idNeutrPart_withKs', 'int')
nt, idStrNeutrPart_withKs = addVect(nt, 'idStrNeutrPart_withKs', 'string')
nt, nChargedPart_withKs = addVar(nt, 'nChargedPart_withKs', 'i')
nt, idChargedPart_withKs = addVect(nt, 'idChargedPart_withKs', 'int')
nt, idStrChargedPart_withKs = addVect(nt, 'idStrChargedPart_withKs', 'string')
nt, nPart_withKs = addVar(nt, 'nPart_withKs', 'i')
nt, idPart_withKs = addVect(nt, 'idPart_withKs', 'int')
nt, idStrPart_withKs = addVect(nt, 'idStrPart_withKs', 'string')
### With KL
nt, nNeutrPart_withKL = addVar(nt, 'nNeutrPart_withKL', 'i')
nt, idNeutrPart_withKL = addVect(nt,'idNeutrPart_withKL', 'int')
nt, idStrNeutrPart_withKL = addVect(nt,'idStrNeutrPart_withKL', 'string')
nt, nChargedPart_withKL = addVar(nt,'nChargedPart_withKL', 'i')
nt, idChargedPart_withKL = addVect(nt,'idChargedPart_withKL', 'int')
nt, idStrChargedPart_withKL = addVect(nt,'idStrChargedPart_withKL', 'string')
nt, nPart_withKL = addVar(nt, 'nPart_withKL', 'i')
nt, idPart_withKL = addVect(nt, 'idPart_withKL', 'int')
nt, idStrPart_withKL = addVect(nt, 'idStrPart_withKL', 'string')
### With N
nt, nNeutrPart_withN = addVar(nt,'nNeutrPart_withN', 'i')
nt, idNeutrPart_withN = addVect(nt,'idNeutrPart_withN', 'int')
nt, idStrNeutrPart_withN = addVect(nt,'idStrNeutrPart_withN', 'string')
nt, nChargedPart_withN = addVar(nt, 'nChargedPart_withN', 'i')
nt, idChargedPart_withN = addVect(nt,'idChargedPart_withN', 'int')
nt, idStrChargedPart_withN = addVect(nt,'idStrChargedPart_withN', 'string')
nt, nPart_withN = addVar(nt, 'nPart_withN', 'i')
nt, idPart_withN = addVect(nt, 'idPart_withN', 'int')
nt, idStrPart_withN = addVect(nt, 'idStrPart_withN', 'string')
#is there a Ks?
nt, hasKs = addVar(nt, 'hasKs', 'i')
nt, hasKL = addVar(nt, 'hasKL', 'i' )
nt, hasN = addVar(nt, 'hasN', 'i')
nt, rpcPoint_x = addVect(nt, 'rpcPoint_x', 'float')
nt, rpcPoint_y = addVect(nt, 'rpcPoint_y', 'float')
nt, rpcPoint_z = addVect(nt, 'rpcPoint_z', 'float')
nt, strawPoint_x = addVect(nt, 'strawPoint_x', 'float')
nt, strawPoint_y = addVect(nt, 'strawPoint_y', 'float')
nt, strawPoint_z = addVect(nt, 'strawPoint_z', 'float')
nt, strawVetoPoint_x = addVect(nt, 'strawVetoPoint_x', 'float')
nt, strawVetoPoint_y = addVect(nt, 'strawVetoPoint_y', 'float')
nt, strawVetoPoint_z = addVect(nt, 'strawVetoPoint_z', 'float')
nt, scintWPoint_x = addVect(nt, 'scintWPoint_x', 'float')
nt, scintWPoint_y = addVect(nt, 'scintWPoint_y', 'float')
nt, scintWPoint_z = addVect(nt, 'scintWPoint_z', 'float')
'''
# was the liquid scintillator window fired?
nt, scintW = addVar(nt,'scintW', 'i' )
# accounting for an eff. of each station of 90%
nt, scintW_eff = addVar(nt,'scintW_eff', 'i' )
# was the liquid scintillator wall (around the volume) fired?
nt, scintV = addVar(nt,'scintV', 'i' )
# accounting for an eff. of each station of 90%
nt, scintV_eff = addVar(nt,'scintV_eff', 'i' )
'''
nt, maxPcharged = addVar(nt, 'maxPcharged', 'f')
nt, maxPneutr = addVar(nt, 'maxPneutr', 'f')
### RPC ###
# was at least an RPC fired by something coming from the nu interaction?
nt, RPC = addVar(nt,'RPC', 'i' )
# accounting for an eff. of each station of 90%
nt, RPC_eff = addVar(nt,'RPC_eff', 'i' )
nt, nRPCstations = addVar(nt, 'nRPCstations', 'i')
nt, nRPCstations_eff = addVar(nt, 'nRPCstations_eff', 'i')
## In case something (no matter if it comes from the nu-interaction kids) fired the RPC
nt, nRPCstationsAny = addVar(nt, 'nRPCstationsAny', 'i')
nt, RPCany = addVar(nt,'RPCany', 'i' )
nt, nRPCstationsAny_eff = addVar(nt, 'nRPCstationsAny_eff', 'i')
nt, RPCany_eff = addVar(nt,'RPCany_eff', 'i' )
############
### scintW ###
# was at least an scintW fired by something coming from the nu interaction?
nt, scintW = addVar(nt,'scintW', 'i' )
# accounting for an eff. of each station of 90%
nt, scintW_eff = addVar(nt,'scintW_eff', 'i' )
nt, nScintWstations = addVar(nt, 'nScintWstations', 'i')
nt, nScintWstations_eff = addVar(nt, 'nScintWstations_eff', 'i')
## In case something (no matter if it comes from the nu-interaction kids) fired the scintW
nt, nScintWstationsAny = addVar(nt, 'nScintWstationsAny', 'i')
nt, scintWany = addVar(nt,'scintWany', 'i' )
nt, nScintWstationsAny_eff = addVar(nt, 'nScintWstationsAny_eff', 'i')
nt, scintWany_eff = addVar(nt,'scintWany_eff', 'i' )
###############
### strawVeto ###
# was at least an strawVeto fired by something coming from the nu interaction?
nt, strawVeto = addVar(nt,'strawVeto', 'i' )
# accounting for an eff. of each station of 90%
nt, strawVeto_eff = addVar(nt,'strawVeto_eff', 'i' )
nt, nStrawVetoStations = addVar(nt, 'nStrawVetoStations', 'i')
nt, nStrawVetoStations_eff = addVar(nt, 'nStrawVetoStations_eff', 'i')
## In case something (no matter if it comes from the nu-interaction kids) fired the strawVeto
nt, nStrawVetoStationsAny = addVar(nt, 'nStrawVetoStationsAny', 'i')
nt, strawVetoAny = addVar(nt,'strawVetoAny', 'i' )
nt, nStrawVetoStationsAny_eff = addVar(nt, 'nStrawVetoStationsAny_eff', 'i')
nt, strawVetoAny_eff = addVar(nt,'strawVetoAny_eff', 'i' )
###############
### TrackSyst ###
## only the case of anything fired it is considered, obviously
nt, TrackSyst = addVar(nt, "TrackSyst", 'i')
nt, TrackSyst_eff = addVar(nt, "TrackSyst_eff", 'i')
nt, nTrackStations = addVar(nt, 'nTrackStations','i')
nt, nTrackStations_eff = addVar(nt, 'nTrackStations_eff','i')
##################
# did the neutrino interact in the last passive element of the "opera-mu-system"?
nt, lastPassive = addVar(nt, 'lastPassive', 'i' )
nt, nuKid_z = addVect(nt, "nuKid_z", "float")
#[1754,]:
# [8418,]:#
for entry in xrange(entries):
#entry = 5429
if not (entry%1000):
print "%s / %s"%(entry,entries)
#if debug:
#if (entry%1000):
# continue
#res = {}
t.GetEntry(entry)
event[0] = entry
particles = t.MCTrack
rpc = t.ShipRpcPoint
scint = t.vetoPoint
strawPoints = t.strawtubesPoint
vetoScint = t.vetoPoint
#print "checking rpc: entry %s"%entry
#PrintRPChits(rpc,pdg)
## initialization
putToZero(nNu)
putToZero(nPart_fromNu)
putToZero(nChargedPart_fromNu)
putToZero(nNeutrPart_fromNu)
putToZero(nNeutrPart_withKs)
putToZero(nChargedPart_withKs)
putToZero(nPart_withKs)
putToZero(nNeutrPart_withKL)
putToZero(nChargedPart_withKL)
putToZero(nPart_withKL)
putToZero(nNeutrPart_withN)
putToZero(nChargedPart_withN)
putToZero(nPart_withN)
putToZero(hasKs)
putToZero(hasKL)
putToZero(hasN)
idPart_fromNu.clear()
idChargedPart_fromNu.clear()
idNeutrPart_fromNu.clear()
idStrPart_fromNu.clear()
idStrChargedPart_fromNu.clear()
idStrNeutrPart_fromNu.clear()
idNeutrPart_withKs.clear()
idChargedPart_withKs.clear()
idPart_withKs.clear()
idStrNeutrPart_withKs.clear()
idStrChargedPart_withKs.clear()
idStrPart_withKs.clear()
idNeutrPart_withKL.clear()
idChargedPart_withKL.clear()
idPart_withKL.clear()
idStrNeutrPart_withKL.clear()
idStrChargedPart_withKL.clear()
idStrPart_withKL.clear()
idNeutrPart_withN.clear()
idChargedPart_withN.clear()
idPart_withN.clear()
idStrNeutrPart_withN.clear()
idStrChargedPart_withN.clear()
idStrPart_withN.clear()
nuKid_z.clear()
strawVetoPoint_x.clear()
strawVetoPoint_y.clear()
strawVetoPoint_z.clear()
strawPoint_x.clear()
strawPoint_y.clear()
strawPoint_z.clear()
rpcPoint_x.clear()
rpcPoint_y.clear()
rpcPoint_z.clear()
scintWPoint_x.clear()
scintWPoint_y.clear()
scintWPoint_z.clear()
primaryDone=False
isPrimary[0] = int(False)
for (ip,part) in enumerate(particles):
## exit if we have reached the empty part of the array
if not (type(part)==type(ShipMCTrack())):
break
if not pdg.GetParticle(part.GetPdgCode()):
#%print "--> to be added before continuing:", part.GetPdgCode()
#sys.exit(0)
pdg.AddParticle("%s"%part.GetPdgCode(),"%s"%part.GetPdgCode(),0.1,kFALSE,0,0,"Isotope",part.GetPdgCode());
print "ERROR: you should insert proper values for: %s"%part.GetPdgCode()
pdgPart = pdg.GetParticle(part.GetPdgCode())
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):
## Starting the counter of how many particles were produced by the interaction of this specific nu
#fromThisNu.push_back(0)
nNu[0]+=1
if part.GetMotherId()==-1:
assert(primaryDone==False)
primaryDone=True
isPrimary[0] = int(True)
assert(len(idStrPart_fromNu)==0)
assert(len(idStrNeutrPart_fromNu)==0)
else:
continue
startZ_nu[0] = part.GetStartZ()
startY_nu[0] = part.GetStartY()
startX_nu[0] = part.GetStartX()
nuE[0] = part.GetEnergy()
nuP[0] = part.GetP()
# Re-initialization for each interacting neutrino
## ---> <--- ##
## Looking for particles produced by the neutrino interaction
interacted = False
partKidsId = []
maxPcharged[0] = 0
maxPneutr[0] = 0
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)
tmpP = part2.GetP()
#print "pID",part2Id, pdg.GetParticle(part2Id).GetName()
nPart_fromNu[0]+=1
idPart_fromNu.push_back(part2Id)
idStrPart_fromNu.push_back(getPartName(part2Id))
## Counting how many particles produced by the nu-interaction are charged or neutral
if int(fabs(pdg.GetParticle(part2Id).Charge()))==int(0):
nNeutrPart_fromNu[0]+=1
idNeutrPart_fromNu.push_back(part2Id)
idStrNeutrPart_fromNu.push_back(getPartName(part2Id))
if fabs(maxPcharged[0])<fabs(tmpP):
maxPcharged[0]=tmpP
else:
nChargedPart_fromNu[0]+=1
idChargedPart_fromNu.push_back(part2Id)
idStrChargedPart_fromNu.push_back(getPartName(part2Id))
if fabs(maxPneutr[0])<fabs(tmpP):
maxPneutr[0]=tmpP
## How many are produced by the interaction in the last passive element of the "opera-mu system"?
part2Z = part2.GetStartZ()
nuKid_z.push_back(part2Z)
somewhere = False
if part2Z>lastPassiveEl[0] and part2Z<lastPassiveEl[1]:
somewhere = True
interactionElement[0] = 0
intElName = 'OPERA'
#RpcPassiveFlag = True
lastPassive[0] = int(True)
## To know if it was in the full OPERA-system (excluded last passive)
if part2Z>OPERA[0] and part2Z<OPERA[1] and somewhere==False:
somewhere = True
interactionElement[0] = 3
intElName = 'OPERA'
## To account for things that are both in the correct and wrong OPERA place
#if part2Z>OPERA_wrong[0] and part2Z<OPERA_wrong[1]:
# interactionElement.push_back(999)
## To account for things that are only in the wrong OPERA place
if part2Z>OPERA_wrong[0] and part2Z<OPERA_wrong[1] and somewhere==False:
somewhere = True
interactionElement[0] = 999
intElName = 'OPERA'
## To know if it was between the two windows
if part2Z>entranceWindows[0][1] and part2Z<entranceWindows[1][0]:
somewhere = True
interactionElement[0] = 4
#scintElement = False
for (iew,ew) in enumerate(entranceWindows):
if part2Z>=ew[0] and part2Z<=ew[1]:
somewhere = True
interactionElement[0] = 1
#scintElement = True
if iew==0:
intElName = 'window1'
else:
intElName = 'window2'
#scintFlag = False
'''if scintElement:
for scintEl in scint:
if scintTankW[0]<scintEl.GetZ()<scintTankW[1]:
scintFlag=True
if scintFlag:
scintW[0] = int(True)
eff_flag = gRandom.Uniform(0.,1.)
if eff_flag<0.9:
scintW_eff[0] = int(True)
else:
scintW_eff[0] = int(False)
else:
scintW[0] = int(False)
'''
#scintElement = False
if part2Z>volume[0] and part2Z<volume[1] and somewhere==False: ## extra check needed for the windows, and for be sure that x and y are as expected
interactionElement[0] = 2
#scintElement = True
somewhere = True
intElName = 'volume1'
'''if scintElement:
for scintEl in scint:
if scintTankV[0]<scintEl.GetZ()<scintTankV[1]:
scintFlag=True
if scintFlag:
scintV[0] = int(True)
eff_flag = gRandom.Uniform(0.,1.)
if eff_flag<0.9:
scintV_eff[0] = int(True)
else:
scintV_eff[0] = int(False)
else:
scintV[0] = int(False)
'''
if somewhere==False:
interactionElement[0] = -1
## Knowing how many particles comes with the KL
if fabs(part2Id) == KLID:
hasKL[0] = int(True)
for (ip3,part3) in enumerate(particles):
if not (type(part)==type(ShipMCTrack())):
break
## same mother has the KL and is not the KL itself
if part3.GetMotherId()==ip and part3 is not part2:
nPart_withKL[0] += 1
part3Id = part3.GetPdgCode()
idPart_withKL.push_back(part3Id)
idStrPart_withKL.push_back(getPartName(part3Id))
if (int(pdg.GetParticle(part3Id).Charge())==int(0)):
nNeutrPart_withKL[0]+=1
idNeutrPart_withKL.push_back(part3Id)
idStrNeutrPart_withKL.push_back(getPartName(part3Id))
else:
nChargedPart_withKL[0]+=1
idChargedPart_withKL.push_back(part3Id)
idStrChargedPart_withKL.push_back(getPartName(part3Id))
## Knowing how many particles comes with the Ks
elif fabs(part2Id) == KsID:
hasKs[0] = int(True)
for (ip3,part3) in enumerate(particles):
if not (type(part)==type(ShipMCTrack())):
break
## same mother has the Ks and is not the Ks itself
if part3.GetMotherId()==ip and part3 is not part2:
nPart_withKs[0] += 1
part3Id = part3.GetPdgCode()
idPart_withKs.push_back(part3Id)
idStrPart_withKs.push_back(getPartName(part3Id))
if (int(pdg.GetParticle(part3Id).Charge())==int(0)):
nNeutrPart_withKs[0]+=1
idNeutrPart_withKs.push_back(part3Id)
idStrNeutrPart_withKs.push_back(getPartName(part3Id))
else:
nChargedPart_withKs[0]+=1
idChargedPart_withKs.push_back(part3Id)
idStrChargedPart_withKs.push_back(getPartName(part3Id))
## Knowing how many particles comes with the N
elif fabs(part2Id) == NID:
hasN[0] = int(True)
for (ip3,part3) in enumerate(particles):
if not (type(part)==type(ShipMCTrack())):
break
## same mother has the N and is not the N itself
if part3.GetMotherId()==ip and part3 is not part2:
nPart_withN[0] += 1
part3Id = part3.GetPdgCode()
idPart_withN.push_back(part3Id)
idStrPart_withN.push_back(getPartName(part3Id))
if (int(pdg.GetParticle(part3Id).Charge())==int(0)):
nNeutrPart_withN[0]+=1
idNeutrPart_withN.push_back(part3Id)
idStrNeutrPart_withN.push_back(getPartName(part3Id))
else:
nChargedPart_withN[0]+=1
idChargedPart_withN.push_back(part3Id)
idStrChargedPart_withN.push_back(getPartName(part3Id))
####
#assert(interacted)
## In case it was in the OPERA system (wrong && correct place)
''' if somewhere and interactionElement[-1] in [0,3,999]:
RPC[0], RPC_eff[0], nRPCstations[0], nRPCstations_eff[0] = wasFired(partKidsId, rpc, RPCstationsPos, checkOn=False)
RPCany[0], RPCany_eff[0], nRPCstationsAny[0], nRPCstationsAny_eff[0] = wasFired(None, rpc, RPCstationsPos, checkOn=False, pointVects=[rpcPoint_x, rpcPoint_y, rpcPoint_z])
if somewhere and interactionElement[-1] in [1,]:
scintW[0], scintW_eff[0], nScintWstations[0], nScintWstations_eff[0] = wasFired(partKidsId, vetoScint, scintWPos, checkOn=False)
scintWany[0], scintWany_eff[0], nScintWstationsAny[0], nScintWstationsAny_eff[0] = wasFired(None, vetoScint, scintWPos, checkOn=False, pointVects=[scintWPoint_x, scintWPoint_y, scintWPoint_z])'''
#if entry==744:
# PrintEventPart(particles,pdg)
#assert(False)
## Looking at hits in the tracking stations
TrackSyst[0],TrackSyst_eff[0],nTrackStations[0],nTrackStations_eff[0] = wasFired(None, strawPoints, trackStationsPos, checkOn=False, pointVects=[strawPoint_x, strawPoint_y, strawPoint_z] )
strawVeto[0],strawVeto_eff[0],nStrawVetoStations[0],nStrawVetoStations_eff[0] = wasFired(partKidsId, strawPoints, strawVetoPos, checkOn=False )
strawVetoAny[0],strawVetoAny_eff[0],nStrawVetoStationsAny[0],nStrawVetoStationsAny_eff[0] = wasFired(None, strawPoints, strawVetoPos, checkOn=False, pointVects=[strawVetoPoint_x, strawVetoPoint_y, strawVetoPoint_z])
scintW[0], scintW_eff[0], nScintWstations[0], nScintWstations_eff[0] = wasFired(partKidsId, vetoScint, scintWPos, checkOn=False)
scintWany[0], scintWany_eff[0], nScintWstationsAny[0], nScintWstationsAny_eff[0] = wasFired(None, vetoScint, scintWPos, checkOn=False, pointVects=[scintWPoint_x, scintWPoint_y, scintWPoint_z])
RPC[0], RPC_eff[0], nRPCstations[0], nRPCstations_eff[0] = wasFired(partKidsId, rpc, RPCstationsPos, checkOn=False)
RPCany[0], RPCany_eff[0], nRPCstationsAny[0], nRPCstationsAny_eff[0] = wasFired(None, rpc, RPCstationsPos, checkOn=False, pointVects=[rpcPoint_x, rpcPoint_y, rpcPoint_z])
assert(len(idStrPart_fromNu)==len(idPart_fromNu))
assert(len(idStrChargedPart_fromNu)==len(idChargedPart_fromNu))
assert(len(idStrNeutrPart_fromNu)==len(idNeutrPart_fromNu))
assert(len(idStrPart_fromNu)==nChargedPart_fromNu[0]+nNeutrPart_fromNu[0])
assert(len(idStrChargedPart_fromNu)==nChargedPart_fromNu[0])
assert(len(idStrNeutrPart_fromNu)==nNeutrPart_fromNu[0])
if not primaryDone:
print "It looks like there is not primary nu interacting"
PrintEventPart(particles,pdg)
assert(False)
nt.SetWeight(calcWeight(intElName, nuE[0], params, entries, ON))
#PrintEventPart(particles,pdg)
#assert(False)
nt.Fill()
nt.Write()
nf.Save()
nf.Close()
print "Wrote file %s"%outputFileName
'''f2 = TFile("prova.root")
t2 = f2.Get("t")
t2.Print()
t2.Draw("nuE", "nNu==1 && isPrimary==1")
#t2.Draw("startZ_nu", "nNu==1 && isPrimary==1 && startZ_nu>=-2501.8 && startZ_nu<=-2478.8")
'''
'''### Looking for the first RPC station that can be fired (i.e. after the interaction point)
leftRPC = None#tmp_firstRPCStation = None
rightRPC = None
for (indexRPC,RPCs) in enumerate(RPCstationsPos):
if part2Z<RPCs[1]:
#tmp_firstRPCStation = indexRPC
leftRPC = indexRPC-1
rightRPC = indexRPC
break
#if tmp_firstRPCStation is None:
if leftRPC is None:
print "ERROR: need to debug event %s"%entry
print "part2z = ",part2Z
print "RPCstationPos = ",RPCstationsPos
sys.exit(0)
#if tmp_firstRPCStation is not 0:
# print "Did you change and you are not looking anymore at the last passive element? If not it should be a bug in the code, if yes remove this part (used only as a safety check)"
#sys.exit(0)
print "###", leftRPC, rightRPC, RPCstationsPos[leftRPC], RPCstationsPos[rightRPC], part2Z
# I used a list in case we would like to include more than the two RPC layers around the passive element where we had the interaction
potentiallyFiredRPCindex = [leftRPC,rightRPC]
## checking if I have one of the two active stations around the interaction point with hits
RPCflag = False
RPCflag_eff = False
for RPCpt in rpc:
print RPCpt.GetZ(), pdg.GetParticle(RPCpt.PdgCode()).GetName(), RPCpt.GetTrackID()
assert(False)
for RPCpt in rpc:
for potRPCindex in potentiallyFiredRPCindex:
print RPCpt.GetZ(), RPCstationsPos[potRPCindex]
if RPCpt.GetZ()<RPCstationsPos[potRPCindex][1] and RPCpt.GetZ()>RPCstationsPos[potRPCindex][0]:
RPCflag = True
## For each station after the interaction point I simulate a 90% efficiency of the detector
## only if one station is still fired I put a positive flag
eff_flag = gRandom.Uniform(0.,1.)
if eff_flag<0.9:
RPCflag_eff = True
print RPCflag, RPCflag_eff
for (pid,part) in enumerate(particles):
print pid, pdg.GetParticle(part.GetName(), part.GetMotherId()
assert(False)
'''
Elena Graverini