shipPatRec_prev Namespace Reference

Functions
	initialize (fGeo)
	getReconstructibleTracks (iEvent, sTree, sGeo)
	SmearHits (iEvent, sTree, modules, SmearedHits, ReconstructibleMCTracks)
	Digitization (sTree, SmearedHits)
	PatRec (firsttwo, zlayer, zlayerv2, StrawRaw, StrawRawLink, ReconstructibleMCTracks)
	TrackFit (hitPosList, theTrack, charge, pinv)
	ptrack (zlayer, ptrackhits, nrwant, window)
	dorot (x, y, xc, yc, alpha)
	line2plane (fitt, fitc, uvview, zuv)
	fitline (indices, xhits, zhits, resolution)
	IP (s, t, b)
	fracMCsame (trackids)
	match_tracks (t1, t2, zmagnet, Bdl)
	dist2line (tan, cst, points)
	hit2wire (ahit, bot, top, no_amb=None)
	debugevent (stereo, y2, y3, ymin2, yother, z2, z3, zmin2, zother, foundtrackids)
	execute (SmearedHits, sTree, ReconstructibleMCTracks)
	finalize ()

Variables
int	reconstructiblehorizontalidsfound12 = 0
int	reconstructiblestereoidsfound12 = 0
int	reconstructiblehorizontalidsfound34 = 0
int	reconstructiblestereoidsfound34 = 0
int	reconstructibleidsfound12 = 0
int	reconstructibleidsfound34 = 0
int	morethan500 = 0
int	morethan100tracks = 0
int	falsenegative = 0
int	falsepositive = 0
int	reconstructibleevents = 0
int	cheated = 0
int	monitor = 0
int	printhelp = 0
int	reconstructiblerequired = 2
int	threeprong = 0
str	geoFile = ''
str	fgeo = ''
int	totalaftermatching = 0
int	totalafterpatrec = 0
list	ReconstructibleMCTracks = []
list	MatchedReconstructibleMCTracks = []
list	fittedtrackids = []
list	theTracks = []
	random = ROOT.TRandom()
	PDG = ROOT.TDatabasePDG.Instance()
	fitter = ROOT.genfit.DAF()
dict	h = {}
dict	rc = h['pinvvstruepinv'].SetMarkerStyle(8)
list	particles = ["e-","e+","mu-","mu+","pi-","pi+","other"]
int	i1 = 1
int	i2 = 16
dict	zlayer = {}
dict	zlayerv2 = {}
dict	z34layer = {}
dict	z34layerv2 = {}
int	TStation1StartZ = 0.
int	TStation4EndZ = 0.
int	VetoStationZ = 0.
int	VetoStationEndZ = 0.

Function Documentation

debugevent()

shipPatRec_prev.debugevent	(		stereo,
			y2,
			y3,
			ymin2,
			yother,
			z2,
			z3,
			zmin2,
			zother,
			foundtrackids
	)

Definition at line 1620 of file shipPatRec_prev.py.

def debugevent(stereo,y2,y3,ymin2,yother,z2,z3,zmin2,zother,foundtrackids):   
     c = ROOT.TCanvas("c","c",600, 400)
     if stereo==False:
        coord="y"
     else:
        coord="projected x" 
     mg=ROOT.TMultiGraph(coord+"-hit coord vs z",coord+"-hit coord vs z; evtnb="+str(iEvent))
     y3vector=ROOT.TVector(len(y3))
     z3vector=ROOT.TVector(len(z3))
     for i in range(len(z3)):
        y3vector[i]=y3[i]
        z3vector[i]=z3[i]
     y2vector=ROOT.TVector(len(y2))
     z2vector=ROOT.TVector(len(z2))
     for i in range(len(z2)):
        y2vector[i]=y2[i]
        z2vector[i]=z2[i]   
     ymin2vector=ROOT.TVector(len(ymin2))
     zmin2vector=ROOT.TVector(len(zmin2))
     for i in range(len(zmin2)):
        ymin2vector[i]=ymin2[i]
        zmin2vector[i]=zmin2[i] 
     yothervector=ROOT.TVector(len(yother))
     zothervector=ROOT.TVector(len(zother)) 
     for i in range(len(zother)):
        yothervector[i]=yother[i]
        zothervector[i]=zother[i]   
     if len(z3)>0:    
        g3=ROOT.TGraph(z3vector,y3vector)
        g3.SetName("g3")
        g3.SetTitle("TrackID=3")    
        g3.SetMarkerStyle(3)  
        g3.SetDrawOption("AP")   
        g3.SetFillStyle(0)       
        mg.Add(g3)
     if len(z2)>0:
        g2=ROOT.TGraph(z2vector,y2vector)   
        g2.SetName("g2")
        g2.SetTitle("TrackID=2")
        g2.SetMarkerStyle(2)  
        g2.SetDrawOption("AP")
        g2.SetFillStyle(0)   
        mg.Add(g2)
     if len(zmin2)>0:
        gmin2=ROOT.TGraph(zmin2vector,ymin2vector)
        gmin2.SetName("gmin2")
        gmin2.SetMarkerStyle(4)  
        gmin2.SetTitle("TrackID=-2")    
        gmin2.SetDrawOption("AP")
        gmin2.SetFillStyle(0)  
        mg.Add(gmin2)
     if len(zother) > 0: 
        gother=ROOT.TGraph(zothervector,yothervector)
    gother.SetName("gother")
        gother.SetMarkerStyle(5)  
        gother.SetDrawOption("AP")
    gother.SetTitle("TrackID=Other")
    gother.SetFillStyle(0)  
    mg.Add(gother)  
     mg.Draw("A P")
     c.BuildLegend(0.6,0.3,0.99,0.5,"Recble TrackIDs="+str(ReconstructibleMCTracks)+"Found TrackIDs:"+str(foundtrackids))
     c.Write()
     return
   

Digitization()

shipPatRec_prev.Digitization	(		sTree,
			SmearedHits
	)

Definition at line 660 of file shipPatRec_prev.py.

def Digitization(sTree,SmearedHits):
  #digitization
  #input: Smeared TrackingHits 
  #output: StrawRaw, StrawRawLink
      
  StrawRaw={} #raw hit dictionary: key=hit number, values=xtop,ytop,ztop,xbot,ybot,zbot,dist2Wire,strawname
  StrawRawLink={} #raw hit dictionary: key=hit number, value=the hit object
  j=0
  for i in range(len(SmearedHits)):  
    xtop=SmearedHits[i]['xtop']
    xbot=SmearedHits[i]['xbot']
    ytop=SmearedHits[i]['ytop']
    ybot=SmearedHits[i]['ybot']
    ztop=SmearedHits[i]['z'] 
    zbot=SmearedHits[i]['z']
    distance=SmearedHits[i]['dist']
    strawname=sTree.strawtubesPoint[i].GetDetectorID()
    a=[]   
    a.append(xtop)
    a.append(ytop)
    a.append(ztop)
    a.append(xbot)
    a.append(ybot)
    a.append(zbot)  
    a.append(float(distance))
    a.append(str(strawname))
    StrawRaw[j]=a
    StrawRawLink[j]=[sTree.strawtubesPoint[i]]
    j=j+1 
 
  if global_variables.debug:
    print("Nbr of digitized hits", j)  
  return StrawRaw,StrawRawLink
         

dist2line()

shipPatRec_prev.dist2line	(		tan,
			cst,
			points
	)

Definition at line 1602 of file shipPatRec_prev.py.

def dist2line(tan,cst,points):
  #points = (x,y)
  #tan, cst are the tangens & constant of the fitted track
  dist = tan*points[0][0] + cst - points[1]
  return dist
 

dorot()

shipPatRec_prev.dorot	(		x,
			y,
			xc,
			yc,
			alpha
	)

Definition at line 1472 of file shipPatRec_prev.py.

def dorot(x,y,xc,yc,alpha):
    #rotate x,y around xc,yc over angle alpha
    ca=cos(alpha)
    sa=sin(alpha)
    xout=ca*(x-xc)-sa*(y-yc)+xc
    yout=sa*(x-xc)+ca*(y-yc)+yc
    return xout,yout
 

execute()

shipPatRec_prev.execute	(		SmearedHits,
			sTree,
			ReconstructibleMCTracks
	)

Definition at line 1684 of file shipPatRec_prev.py.

def execute(SmearedHits,sTree,ReconstructibleMCTracks):
 global totalaftermatching,morethan500,falsepositive,falsenegative,totalafterpatrec
 global reconstructibleevents,morethan100tracks,theTracks
 
 fittedtrackids=[]  
 reconstructibles12=0
 reconstructibles34=0
 theTracks=[]
 
 if global_variables.debug:
   print("************* PatRect START **************")  
    
 nShits=sTree.strawtubesPoint.GetEntriesFast() 
 nMCTracks = sTree.MCTrack.GetEntriesFast() 
 
 #if nMCTracks < 100:
 if nShits < 500: 
   
     StrawRaw,StrawRawLink=Digitization(sTree,SmearedHits)
     reconstructibleevents+=1  
     nr12tracks,tracks12,hitids12,xtan12,xcst12,stereotan12,stereocst12,px12,py12,pz12,fraction12,trackid12=PatRec(True,zlayer,zlayerv2,StrawRaw,StrawRawLink,ReconstructibleMCTracks)
     tracksfound=[]
     if monitor==True:
       for item in ReconstructibleMCTracks:
     for value in trackid12.values():  
        if item == value and item not in tracksfound:
            reconstructibles12+=1
        tracksfound.append(item)
        
       if global_variables.debug:
     print(
         "tracksfound", tracksfound,
         "reconstructibles12", reconstructibles12,
         "ReconstructibleMCTracks", ReconstructibleMCTracks
         )  
       if len(tracksfound)==0 and len(ReconstructibleMCTracks)>0: 
         if global_variables.debug:
       print("No tracks found in event after stations 1&2. Rejecting event.")
         return
    
     nr34tracks,tracks34,hitids34,xtan34,xcst34,stereotan34,stereocst34,px34,py34,pz34,fraction34,trackid34=PatRec(False,z34layer,z34layerv2,StrawRaw,StrawRawLink,ReconstructibleMCTracks)
     
     tracksfound=[]      
     if monitor==True:
      for item in ReconstructibleMCTracks:
    for value in trackid34.values():  
        if item == value and item not in tracksfound:
            reconstructibles34+=1 
        tracksfound.append(item)
      if len(tracksfound)==0 and len(ReconstructibleMCTracks)>0: 
        if global_variables.debug:
      print("No tracks found in event after stations 3&4. Rejecting event.")
        return
      MatchedReconstructibleMCTracks=len(ReconstructibleMCTracks)*[0]  
 
 
     if global_variables.debug:
      if (nr12tracks>0) : 
       print(nr12tracks,"tracks12  ",tracks12)
       print("hitids12  ",hitids12)
       print("xtan  ",xtan12,"xcst",xcst12)
       print("stereotan  ",stereotan12,"stereocst  ",stereocst12)
       print("fraction12",fraction12,"trackid12",trackid12)
      else : 
       print("No tracks found in stations 1&2.")
     
      if (nr34tracks>0) : 
       print(nr34tracks,"tracks34  ",tracks34)
       print("hitids34  ",hitids34)   
       print("xtan34 ",xtan34,"xcst34 ",xcst34)
       print("stereotan34  ",stereotan34,"stereocst34  ",stereocst34)
       print("px34",px34)
      else : 
       print("No tracks found in stations 3&4.")
 
     if monitor==False: totalafterpatrec+=1    
         
     zmagnet=ShipGeo.Bfield.z
     tracksfound=[]
     matches=0
 
     for item in xtan12: 
      z1=0.
      tgy1=xtan12[item]
      tgx1=stereotan12[item]
      y1=xcst12[item]
      x1=stereocst12[item]
      t1=[x1,y1,z1,tgx1,tgy1]
      for ids in hitids12[item]:
         #cheated
         if ids != 999 :
            #loop until we get the particle of this track  
        particle12   = PDG.GetParticle(StrawRawLink[ids][0].PdgCode())
        try:
           charge12=particle12.Charge()/3.
           break
        except:
           continue   
         
      falsenegativethistrack=0
      falsepositivethistrack=0
      for item1 in xtan34:  
         if monitor==True:  totalafterpatrec+=1          
         for k,v in enumerate(ReconstructibleMCTracks):
        if v not in tracksfound: 
          tracksfound.append(v)      
         rawlinkindex=0            
         for ids in hitids34[item1]:
            if ids != 999 :  
           particle34   = PDG.GetParticle(StrawRawLink[ids][0].PdgCode())
           try:
             charge34=particle34.Charge()/3.
             break
           except:
             continue   
         z2=0.
         tgx2=stereotan34[item1]
         tgy2=xtan34[item1]
         y2=xcst34[item1]
         x2=stereocst34[item1]
         t2=[x2,y2,z2,tgx2,tgy2] 
     
         dx,dy,alpha,pinv=match_tracks(t1,t2,zmagnet,-0.75)       
         p2x=px34[item1]
         p2y=py34[item1]
         p2z=pz34[item1]
         p2true=1./math.sqrt(p2x**2+p2y**2+p2z**2)
         if global_variables.debug:
       print(
           "Matching 1&2 track", item,
           "with 3&4 track", item1,
           "dx", dx, "dy", dy, "alpha", alpha, "pinv", pinv, "1/p2true", p2true
           )
         rc=h['dx-matchedtracks'].Fill(dx)
         rc=h['dy-matchedtracks'].Fill(dy)
         if ((abs(dx)<20) & (abs(dy)<2)) :
            #match found between track nr item in stations 12 & item1 in stations 34 
        #get charge from deflection and check if it corresponds to the MC truth
        #field is horizontal (x) hence deflection in y
        tantheta=(tgy1-tgy2)/(1+tgy1*tgy2)
        if tantheta>0 : 
           charge="-1" 
           if charge34>0: 
             falsenegative+=1
             falsenegativethistrack=1
        else: 
           charge="1" 
           if charge34<0: 
              falsepositive+=1
              falsepositivethistrack=1
        
            #reject track (and event) if it doesn't correspond to MC truth      
        if monitor==True:
          if (falsepositivethistrack==0 & falsenegativethistrack==0):
            totalaftermatching+=1
          else: 
            if global_variables.debug:
          print("Charge from track deflection doesn't match MC truth. Rejecting it.")
            break
        else:
          if matches==0: 
            matches==1
            totalaftermatching+=1
 
        if global_variables.debug: 
           print("******* MATCH FOUND stations 12 track id",trackid12[item],"(fraction",fraction12[item]*100,"%) and stations 34 track id",trackid34[item1],"(fraction",fraction34[item1]*100,"%)")
           print("******* Reconstructible tracks ids",ReconstructibleMCTracks)
        
        rc=h['matchedtrackefficiency'].Fill(fraction12[item],fraction34[item1])
        for k,v in enumerate(ReconstructibleMCTracks):
          if global_variables.debug:
        print("MatchedReconstructibleMCTracks", MatchedReconstructibleMCTracks, "k", k, "v", v)
          if v not in MatchedReconstructibleMCTracks:
            if global_variables.debug:
          print("ReconstructibleMCTracks", ReconstructibleMCTracks, "trackid34[item1]", trackid34[item1])
            if v==trackid34[item1]: MatchedReconstructibleMCTracks[k]=v
        p2true=p2true*int(charge)
        rc=h['pinvvstruepinv'].Fill(p2true,pinv)
        if math.fabs(pinv) > 0.0 : rc=h['ptrue-p/ptrue'].Fill((pinv-p2true)/pinv)
      
        if cheated==False: 
           hitPosList=[]             
               for ids in range(0,len(hitids12[item])):
             #print "hitids12[",item,"][",ids,"]",hitids12[item][ids]
                 if hitids12[item][ids]!=999 : 
                    m=[]
                if ids in range(0,32):
                  for z in range(0,7):
                   m.append(StrawRaw[hitids12[item][ids]][z])
                hitPosList.append(m)
               for ids in range(0,len(hitids34[item1])):
                 if hitids34[item1][ids]!=999 : 
                m=[]
                if ids in range(0,32):
                  for z in range(0,7):
                    m.append(StrawRaw[hitids34[item1][ids]][z])
                hitPosList.append(m)
               nM = len(hitPosList)
               if nM<25:
                 if global_variables.debug:
           print("Only", nM, "hits on this track. Insufficient for fitting.")
                 return
           #all particles are assumed to be muons
               if int(charge)<0:
                 pdg=13
               else:
                 pdg=-13   
               rep = ROOT.genfit.RKTrackRep(pdg)   
               posM = ROOT.TVector3(0, 0, 0)
               #would be the correct way but due to uncertainties on small angles the sqrt is often negative
               if math.fabs(pinv) > 0.0 : momM = ROOT.TVector3(0,0,int(charge)/pinv)
               else: momM = ROOT.TVector3(0,0,999)   
               covM = ROOT.TMatrixDSym(6)
               resolution = ShipGeo.strawtubes.sigma_spatial  
               for  i in range(3):   covM[i][i] = resolution*resolution
               covM[0][0]=resolution*resolution*100.
               for  i in range(3,6): covM[i][i] = ROOT.TMath.Power(resolution / nM / ROOT.TMath.Sqrt(3), 2)
               # smeared start state  
               stateSmeared = ROOT.genfit.MeasuredStateOnPlane(rep)
               rep.setPosMomCov(stateSmeared, posM, momM, covM)
               # create track
               seedState = ROOT.TVectorD(6)
               seedCov   = ROOT.TMatrixDSym(6)
               rep.get6DStateCov(stateSmeared, seedState, seedCov)
               theTrack=ROOT.genfit.Track(rep, seedState, seedCov) 
           TrackFit(hitPosList,theTrack,charge,pinv)
               fittedtrackids.append(trackid34[item1])
         else :
        if global_variables.debug:
          print(
          "******* MATCH NOT FOUND stations 12 track id", trackid12[item],
          "(fraction", fraction12[item]*100, "%) and stations 34 track id", trackid34[item1],
          "(fraction", fraction34[item1]*100, "%)"
          )
            if trackid12[item]==trackid34[item1] : 
           if global_variables.debug:
         print("trackids the same, but not matched.")
       
 else: 
      #remove events with >500 hits
      morethan500+=1  
      return
      
 #else:
    #morethan100tracks+=1        
    #return
 return fittedtrackids 
 
 

finalize()

shipPatRec_prev.finalize

(

)

Definition at line 1933 of file shipPatRec_prev.py.

def finalize(): 
   if global_variables.debug:
     ut.writeHists(h,"recohists_patrec.root")

fitline()

shipPatRec_prev.fitline	(		indices,
			xhits,
			zhits,
			resolution
	)

Definition at line 1503 of file shipPatRec_prev.py.

def fitline(indices,xhits,zhits,resolution):
   #fit linear function (y=fitt*x+fitc) to list of hits, 
   #"x" is the z-coordinate, "y" is the cSoordinate in tracking plane perp to z
   #use equal weights for the time being, and outlier removal (yet)?
 
   n=0
   sumx=0.
   sumx2=0
   sumxy=0.
   sumy=0.
   sumweight=0.
   sumweightx=0.
   sumweighty=0. 
   Dw=0.     
   term=0.
   for ipl in range(1,len(indices)):
      indx= indices[ipl]
      if indx>-1:
         #n+=1
     #weigh points accordint to their distance to the wire
     weight=1/math.sqrt(xhits[ipl][3][indx]**2+resolution**2)
     x=zhits[ipl][0]
     y=xhits[ipl][0][indx]
     sumweight+=weight
     sumweightx+=weight*x
     sumweighty+=weight*y
   xmean=sumweightx/sumweight
   ymean=sumweighty/sumweight    
   for ipl in range(1,len(indices)):
      indx= indices[ipl]
      if indx>-1:
         n+=1
     weight=1/math.sqrt(xhits[ipl][3][indx]**2+resolution**2)
     x=zhits[ipl][0]
     y=xhits[ipl][0][indx]       
     Dw+=weight*(x-xmean)**2
     term+=weight*(x-xmean)*y
         sumx+=zhits[ipl][0]
         sumx2+=zhits[ipl][0]**2
         sumxy+=xhits[ipl][0][indx]*zhits[ipl][0]
         sumy+=xhits[ipl][0][indx]
   fitt=(1/Dw)*term
   fitc=ymean-fitt*xmean     
   unweightedfitt=(sumxy-sumx*sumy/n)/(sumx2-sumx**2/n)
   unweightedfitc=(sumy-fitt*sumx)/n
   return fitt,fitc
      

fracMCsame()

shipPatRec_prev.fracMCsame

(

trackids

)

Definition at line 1558 of file shipPatRec_prev.py.

def fracMCsame(trackids):
# determine largest fraction of hits which have benn generated by the same # MC true track
#hits: list of pointers to true MC TrackingHits
#return: fraction of hits from same track, and its track-ID
    track={}
    nh=len(trackids)
    for tid in trackids:
       if tid==999: 
         nh-=1
         continue
       if tid in track:
          track[tid]+=1
       else:
        track[tid]=1
    #now get track with largest number of hits
    tmax=max(track, key=track.get)
 
    frac=0.
    if nh>0: frac=float(track[tmax])/float(nh)
    return frac,tmax     
   

getReconstructibleTracks()

shipPatRec_prev.getReconstructibleTracks	(		iEvent,
			sTree,
			sGeo
	)

Definition at line 236 of file shipPatRec_prev.py.

def getReconstructibleTracks(iEvent,sTree,sGeo):
 
  #returns a list of reconstructible tracks for this event
  #call this routine once for each event before smearing
  MCTrackIDs=[]
  rc = sTree.GetEvent(iEvent) 
  nMCTracks = sTree.MCTrack.GetEntriesFast()   
 
  if global_variables.debug:
    print("event nbr", iEvent, "has", nMCTracks, "tracks")
  #1. MCTrackIDs: list of tracks decaying after the last tstation and originating before the first
  for i in reversed(range(nMCTracks)):
     atrack = sTree.MCTrack.At(i) 
     #for 3 prong decays check if its a nu
     if threeprong == 1:    
       if PDG.GetParticle(atrack.GetPdgCode()):          
         if PDG.GetParticle(atrack.GetPdgCode()).GetName()[:5]=="nu_mu":
         if (atrack.GetStartZ() < TStation1StartZ and  atrack.GetStartZ() > VetoStationEndZ) and i not in MCTrackIDs:
            MCTrackIDs.append(i)
         else:  
           if atrack.GetStartZ() > TStation4EndZ :         
             motherId=atrack.GetMotherId() 
         if motherId > -1 : 
           mothertrack=sTree.MCTrack.At(motherId)   
           mothertrackZ=mothertrack.GetStartZ() 
           #this mother track is a HNL decay          
           #track starts inside the decay volume? (after veto, before 1 st tstation)
               if mothertrackZ < TStation1StartZ and mothertrackZ > VetoStationEndZ:
             if motherId not in MCTrackIDs:
               MCTrackIDs.append(motherId)  
     else:    
       #track endpoint after tstations?
       if atrack.GetStartZ() > TStation4EndZ :       
         motherId=atrack.GetMotherId() 
     if motherId > -1 : 
         mothertrack=sTree.MCTrack.At(motherId) 
         mothertrackZ=mothertrack.GetStartZ() 
         #this mother track is a HNL decay        
         #track starts inside the decay volume? (after veto, before 1 st tstation)
             if mothertrackZ < TStation1StartZ and mothertrackZ > VetoStationEndZ:
           if motherId not in MCTrackIDs:
               MCTrackIDs.append(motherId)
  if global_variables.debug:
    print("Tracks with origin in decay volume", MCTrackIDs)  
  if len(MCTrackIDs)==0: return MCTrackIDs
    
  #2. hitsinTimeDet: list of tracks with hits in TimeDet       
  nVetoHits = sTree.vetoPoint.GetEntriesFast() 
  hitsinTimeDet=[]
  for i in range(nVetoHits):
     avetohit = sTree.vetoPoint.At(i)
     #hit in TimeDet?
     if sGeo.FindNode(avetohit.GetX(),avetohit.GetY(),avetohit.GetZ()).GetName() == 'TimeDet_1':
        if avetohit.GetTrackID() not in hitsinTimeDet:
       hitsinTimeDet.append(avetohit.GetTrackID())
     
  #3. Remove tracks from MCTrackIDs that are not in hitsinTimeDet   
  itemstoremove=[]
  for item in MCTrackIDs:
      if threeprong==1:
        #don't remove the nu
        if PDG.GetParticle(sTree.MCTrack.At(item).GetPdgCode()).GetName()[:5]!="nu_mu" and item not in hitsinTimeDet:  
          itemstoremove.append(item)       
      else :
        if item not in hitsinTimeDet:
          itemstoremove.append(item)
  for item in itemstoremove:
      MCTrackIDs.remove(item)             
 
  if global_variables.debug:
    print("Tracks with hits in timedet", MCTrackIDs) 
  if len(MCTrackIDs)==0: return MCTrackIDs
  #4. Find straws that have multiple hits
  nHits = sTree.strawtubesPoint.GetEntriesFast()  
  hitstraws={}
  duplicatestrawhit=[]
  if global_variables.debug:
    print("Nbr of Rawhits=", nHits)
 
  for i in range(nHits):
    ahit = sTree.strawtubesPoint[i]
    if (str(ahit.GetDetectorID())[:1]=="5") : 
       if global_variables.debug:
     print("Hit in straw Veto detector. Rejecting.")
       continue
    strawname=str(ahit.GetDetectorID())
    
    if strawname in hitstraws:
       #straw was already hit
       if ahit.GetX()>hitstraws[strawname][1]:
          #this hit has higher x, discard it
          duplicatestrawhit.append(i)
       else:
          #del hitstraws[strawname] 
      duplicatestrawhit.append(hitstraws[strawname][0])
          hitstraws[strawname]=[i,ahit.GetX()]          
    else: 
       hitstraws[strawname]=[i,ahit.GetX()]
  
  #5. Split hits up by station and outside stations
  hits1={}
  hits2={}
  hits3={}
  hits4={}
  trackoutsidestations=[]
  for i in range(nHits):
    if i in  duplicatestrawhit: 
       if global_variables.debug:
     print("Duplicate hit", i, "not reconstructible, rejecting.")
       continue  
    ahit = sTree.strawtubesPoint[i] 
    #is hit inside acceptance? if not mark the track as bad   
    if (((ahit.GetX()/245.)**2 + (ahit.GetY()/495.)**2) >= 1.): 
       if ahit.GetTrackID() not in trackoutsidestations:
          trackoutsidestations.append(ahit.GetTrackID())
    if ahit.GetTrackID() not in MCTrackIDs:
       #hit on not reconstructible track
       if global_variables.debug:
     print("Hit not on reconstructible track. Rejecting.")
       continue   
    #group hits per tracking station, key = trackid
    if str(ahit.GetDetectorID())[:1]=="1" :
       if ahit.GetTrackID() in hits1:
            hits1[ahit.GetTrackID()]=[hits1[ahit.GetTrackID()][0],i]
       else:  
            hits1[ahit.GetTrackID()]=[i]    
    if str(ahit.GetDetectorID())[:1]=="2" :
       if ahit.GetTrackID() in hits2:
            hits2[ahit.GetTrackID()]=[hits2[ahit.GetTrackID()][0],i]
       else:  
            hits2[ahit.GetTrackID()]=[i]   
    if str(ahit.GetDetectorID())[:1]=="3" :
       if ahit.GetTrackID() in hits3:
            hits3[ahit.GetTrackID()]=[hits3[ahit.GetTrackID()][0],i]
       else:  
            hits3[ahit.GetTrackID()]=[i]           
    if str(ahit.GetDetectorID())[:1]=="4" :
       if ahit.GetTrackID() in hits4:
            hits4[ahit.GetTrackID()]=[hits4[ahit.GetTrackID()][0],i]
       else:  
            hits4[ahit.GetTrackID()]=[i] 
        
  #6. Make list of tracks with hits in in station 1,2,3 & 4         
  tracks_with_hits_in_all_stations=[]  
  for key in hits1.keys():
      if (key in hits2 and key in hits3 ) and key in hits4:
         if key not in tracks_with_hits_in_all_stations and key not in trackoutsidestations:
            tracks_with_hits_in_all_stations.append(key) 
  for key in hits2.keys():
      if (key in hits1 and key in hits3 ) and key in hits4:
         if key not in tracks_with_hits_in_all_stations and key not in trackoutsidestations:
            tracks_with_hits_in_all_stations.append(key) 
  for key in hits3.keys():
      if ( key in hits2 and key in hits1 ) and key in hits4:
         if key not in tracks_with_hits_in_all_stations and key not in trackoutsidestations:
            tracks_with_hits_in_all_stations.append(key) 
  for key in hits4.keys():
      if (key in hits2 and key in hits3) and key in hits1:
         if key not in tracks_with_hits_in_all_stations and key not in trackoutsidestations:
            tracks_with_hits_in_all_stations.append(key) 
 
  #7. Remove tracks from MCTrackIDs with hits outside acceptance or doesn't have hits in all stations             
  itemstoremove=[]
  for item in MCTrackIDs:   
      if threeprong==1:
        #don't remove the nu
        if PDG.GetParticle(sTree.MCTrack.At(item).GetPdgCode()).GetName()[:5]!="nu_mu" and item not in tracks_with_hits_in_all_stations:  
          itemstoremove.append(item)      
      else:
        if item not in tracks_with_hits_in_all_stations:
          itemstoremove.append(item)
  for item in itemstoremove:
      MCTrackIDs.remove(item)   
 
  if global_variables.debug: 
     print("tracks_with_hits_in_all_stations",tracks_with_hits_in_all_stations)
     print("Tracks with hits in all stations & inside acceptance ellipse",MCTrackIDs)   
  if len(MCTrackIDs)==0: return MCTrackIDs   
  nbrechits=0       
  for i in range(nHits):
    if i in  duplicatestrawhit: 
       continue 
    nbrechits+=1    
    ahit = sTree.strawtubesPoint[i]     
    if ahit.GetTrackID()>-1 and ahit.GetTrackID() in MCTrackIDs:           
      atrack = sTree.MCTrack.At(ahit.GetTrackID())
      for j in range(ahit.GetTrackID()+1,nMCTracks) :
        childtrack = sTree.MCTrack.At(j)
        if childtrack.GetMotherId() == ahit.GetTrackID():      
        trackmomentum=atrack.GetP()
        trackweight=atrack.GetWeight()
        rc=h['reconstructiblemomentum'].Fill(trackmomentum,trackweight)
        motherId=atrack.GetMotherId() 
        if motherId==1 :
        HNLmomentum=sTree.MCTrack.At(1).GetP()
        rc=h['HNLmomentumvsweight'].Fill(trackweight,HNLmomentum) 
            if j==nMCTracks :
             trackmomentum=atrack.GetP()
             trackweight=atrack.GetWeight()
             rc=h['reconstructiblemomentum'].Fill(trackmomentum,trackweight)
             if atrack.GetMotherId()==1 :
               HNLmomentum=sTree.MCTrack.At(1).GetP()
               rc=h['HNLmomentumvsweight'].Fill(trackweight,HNLmomentum)
  itemstoremove=[]
  for item in MCTrackIDs:          
    atrack = sTree.MCTrack.At(item)
    motherId=atrack.GetMotherId() 
    if motherId != 1:
        itemstoremove.append(item)
  for item in itemstoremove:
      MCTrackIDs.remove(item)   
      if global_variables.debug:
    print("After removing the non HNL track, MCTrackIDs", MCTrackIDs)
  if global_variables.debug:
    print("Tracks with HNL mother", MCTrackIDs) 
  
  #8. check if the tracks are HNL children 
  mufound=0
  pifound=0 
  nu_mufound=0
  itemstoremove=[]
  if MCTrackIDs:
    for item in MCTrackIDs: 
      try: 
        if PDG.GetParticle(sTree.MCTrack.At(item).GetPdgCode()).GetName()[:2]=="mu"   : mufound+=1  
        if PDG.GetParticle(sTree.MCTrack.At(item).GetPdgCode()).GetName()[:2]=="pi"   : pifound+=1 
        if PDG.GetParticle(sTree.MCTrack.At(item).GetPdgCode()).GetName()[:5]=="nu_mu": 
       nu_mufound+=1  
       itemstoremove.append(item)
      except:
        if global_variables.debug:
      print("Unknown particle with pdg code:", sTree.MCTrack.At(item).GetPdgCode()) 
    if reconstructiblerequired == 1 :
      if mufound!=1  and pifound!=1: 
          if global_variables.debug:
        print("No reconstructible pion or muon.") 
      MCTrackIDs=[]   
    if reconstructiblerequired == 2 : 
      if threeprong == 1 :       
          if mufound!=2 or nu_mufound!=1 : 
            if global_variables.debug:
          print("No reconstructible mu-mu-nu.")  
        MCTrackIDs=[]
      else:
        #remove the neutrino from MCTrackIDs for the rest
        for item in itemstoremove:
               MCTrackIDs.remove(item)  
      else:         
          if mufound!=1 or pifound!=1 : 
            if global_variables.debug:
          print("No reconstructible pion and muon.")  
        MCTrackIDs=[]     
  if len(MCTrackIDs)>0:
     rc=h['nbrhits'].Fill(nHits)
     rc=h['nbrtracks'].Fill(nMCTracks)
  if global_variables.debug:
    print("Tracks with required HNL decay particles", MCTrackIDs)        
  return MCTrackIDs
 

hit2wire()

shipPatRec_prev.hit2wire	(		ahit,
			bot,
			top,
			no_amb = None
	)

Definition at line 1608 of file shipPatRec_prev.py.

def hit2wire(ahit,bot,top,no_amb=None):
     detID = ahit.GetDetectorID()
     ex = ahit.GetX()
     ey = ahit.GetY()
     ez = ahit.GetZ()
   #distance to wire, and smear it.
     dw  = ahit.dist2Wire()
     smear = dw
     if not no_amb: smear = ROOT.fabs(random.Gaus(dw,ShipGeo.strawtubes.sigma_spatial))
     smearedHit = {'mcHit':ahit,'xtop':top.x(),'ytop':top.y(),'z':top.z(),'xbot':bot.x(),'ybot':bot.y(),'z':bot.z(),'dist':smear}
     return smearedHit
 

initialize()

shipPatRec_prev.initialize

(

fGeo

)

Definition at line 127 of file shipPatRec_prev.py.

def initialize(fGeo):
   #creates a dictionary with z coordinates of layers
   #and variables with station start/end coordinates
   #to be called once at the beginning of the eventloop
   global i1,i2,zlayer,zlayerv2,z34layer,z34layerv2,TStation1StartZ,TStation4EndZ,VetoStationZ,VetoStationEndZ
   global fgeo
   fgeo=fGeo 
   #z-positions of Y-view tracking
   #4 stations, 4 views (Y,u,v,Y); each view has 2 planes and each plane has 2 layers
 
   for i in range(i1,i2+1):
     TStationz = ShipGeo.TrackStation1.z
     if (i>8) : 
        TStationz = ShipGeo.TrackStation2.z  
     # Y: vnb=0 or 3   
     vnb=0.
     if (i>4): vnb=3.
     if (i>8): vnb=0.
     if (i>12): vnb=3.
     lnb = 0.  
     if (i % 2 == 0) : lnb=1.
     pnb=0.
     if (i==3 or i==4 or i==7 or i==8 or i==11 or i==12 or i==15 or i==16) : pnb=1.
 
     #z positions of Y view of stations   
     Zpos = TStationz+(vnb-3./2.)*ShipGeo.strawtubes.DeltazView+(float(pnb)-1./2.)*ShipGeo.strawtubes.DeltazPlane+(float(lnb)-1./2.)*ShipGeo.strawtubes.DeltazLayer 
     zlayer[i]=[Zpos]
 
   #z-positions for stereo views
 
   for i in range(i1,i2+1):
     TStationz = ShipGeo.TrackStation1.z
     if (i>8) : 
        TStationz = ShipGeo.TrackStation2.z   
     #stereo views: vnb=1 or 2  
     vnb=1.
     if (i>4): vnb=2.
     if (i>8): vnb=1.
     if (i>12): vnb=2. 
     lnb = 0.  
     if (i % 2 == 0) : lnb=1.
     pnb=0.
     if (i==3 or i==4 or i==7 or i==8 or i==11 or i==12 or i==15 or i==16) : pnb=1.
 
     #z positions of u,v view of stations   
     Zpos_u = TStationz+(vnb-3./2.)*ShipGeo.strawtubes.DeltazView+(float(pnb)-1./2.)*ShipGeo.strawtubes.DeltazPlane+(float(lnb)-1./2.)*ShipGeo.strawtubes.DeltazLayer 
     zlayerv2[i]=[Zpos_u]
 
 
   for i in range(i1,i2+1):
     TStationz = ShipGeo.TrackStation3.z
     if (i>8) : 
        TStationz = ShipGeo.TrackStation4.z  
     # Y: vnb=0 or 3   
     vnb=0.
     if (i>4): vnb=3.
     if (i>8): vnb=0.
     if (i>12): vnb=3. 
     lnb = 0.  
     if (i % 2 == 0) : lnb=1.
     pnb=0.
     if (i==3 or i==4 or i==7 or i==8 or i==11 or i==12 or i==15 or i==16) : pnb=1.
  
     #z positions of x1 view of stations   
     Zpos = TStationz+(vnb-3./2.)*ShipGeo.strawtubes.DeltazView+(float(pnb)-1./2.)*ShipGeo.strawtubes.DeltazPlane+(float(lnb)-1./2.)*ShipGeo.strawtubes.DeltazLayer 
     z34layer[i]=[Zpos]
 
 
   for i in range(i1,i2+1):
     #zlayerv2[i]=[i*100.+50.]
     TStationz = ShipGeo.TrackStation3.z
     if (i>8) : 
        TStationz = ShipGeo.TrackStation4.z   
     #stereo views: vnb=1 or 2  
     vnb=1.
     if (i>4): vnb=2.
     if (i>8): vnb=1.
     if (i>12): vnb=2. 
     lnb = 0.  
     if (i % 2 == 0) : lnb=1.
     pnb=0.
     if (i==3 or i==4 or i==7 or i==8 or i==11 or i==12 or i==15 or i==16) : pnb=1.
 
     #z positions of u,v view of stations   
     Zpos_u = TStationz+(vnb-3./2.)*ShipGeo.strawtubes.DeltazView+(float(pnb)-1./2.)*ShipGeo.strawtubes.DeltazPlane+(float(lnb)-1./2.)*ShipGeo.strawtubes.DeltazLayer 
     z34layerv2[i]=[Zpos_u]
 
   VetoStationZ = ShipGeo.vetoStation.z
   if global_variables.debug:
     print("VetoStation midpoint z=", VetoStationZ)
   VetoStationEndZ=VetoStationZ+(ShipGeo.strawtubes.DeltazView+ShipGeo.strawtubes.OuterStrawDiameter)/2
   for i in range(1,5):   
     if i==1: TStationz = ShipGeo.TrackStation1.z
     if i==2: TStationz = ShipGeo.TrackStation2.z  
     if i==3: TStationz = ShipGeo.TrackStation3.z  
     if i==4: TStationz = ShipGeo.TrackStation4.z 
     if global_variables.debug:
       print("TrackStation",i," midpoint z=",TStationz) 
       for vnb in range(0,4):
         for pnb in range (0,2):
           for lnb in range (0,2):
              Zpos = TStationz+(vnb-3./2.)*ShipGeo.strawtubes.DeltazView+(float(pnb)-1./2.)*ShipGeo.strawtubes.DeltazPlane+(float(lnb)-1./2.)*ShipGeo.strawtubes.DeltazLayer 
              print("TStation=",i,"view=",vnb,"plane=",pnb,"layer=",lnb,"z=",Zpos)
 
   TStation1StartZ=zlayer[1][0]-ShipGeo.strawtubes.OuterStrawDiameter/2
   TStation4EndZ=z34layer[16][0]+ShipGeo.strawtubes.OuterStrawDiameter/2
 
   return 
             

IP()

shipPatRec_prev.IP	(		s,
			t,
			b
	)

Definition at line 1550 of file shipPatRec_prev.py.

def IP(s, t, b):
     #some vector manipulations to get impact point of t-b to s
     dir=ROOT.TVector3(t-b)
     udir = ROOT.TVector3(dir.Unit())
     sep = ROOT.TVector3(t-s)
     ip = ROOT.TVector3(udir.Cross(sep.Cross(udir) ))
     return ip.Mag()
 

line2plane()

shipPatRec_prev.line2plane	(		fitt,
			fitc,
			uvview,
			zuv
	)

Definition at line 1480 of file shipPatRec_prev.py.

def line2plane(fitt,fitc,uvview,zuv):
  #intersect horizontal plane, defined by tangent, constant in yz plane, with
  #stereo hits given in uvview,zuv with top/bottom hits of "straw"
  xb=uvview[0]
  yb=uvview[1]
  xt=uvview[2]
  yt=uvview[3]
  x=[]
  items=[]
  term=fitc+zuv*fitt
  #loop over all hits in this view
  for i in range(len(yb)):
    #f2=(term-yb[i])/(yt[i]-yb[i])
    #xint=xb[i]+f2*(xt[i]-xb[i])
    f2=(term-yb[i])/(yb[i]-yt[i])
    xint=xb[i]+f2*(xb[i]-xt[i])
    c=uvview[4][i]
    #do they cross inside sensitive volume defined by top/bottom of straw?
    if xint<max(xt[i],xb[i]) and xint>min(xb[i],xt[i]): 
      x.append(xint)
      items.append(c)
  return x,items
 

match_tracks()

shipPatRec_prev.match_tracks	(		t1,
			t2,
			zmagnet,
			Bdl
	)

Definition at line 1579 of file shipPatRec_prev.py.

def match_tracks(t1,t2,zmagnet,Bdl):
# match tracks before to after magnet
# t1,t2: x,y,z,tgx,tgy of two input tracks
# zmagnet: middle of the magnet
# BdL: Bdl in Tm
# returns:
#     dx,dy: distance between tracks at zmagnet
#     1./pmom:  innverse momentum*q estimated
 
#linear extrapoaltion to centre of the magnet
   dz1=zmagnet-t1[2]
   x1m=t1[0]+dz1*t1[3]
   y1m=t1[1]+dz1*t1[4]
   dz2=zmagnet-t2[2]
   x2m=t2[0]+dz2*t2[3]
   y2m=t2[1]+dz2*t2[4]
   dx=x1m-x2m
   dy=y1m-y2m
 
   alpha=math.atan(t1[4])-math.atan(t2[4])
   pinv=math.sin(alpha)/(Bdl*0.3)
   return dx,dy,alpha,pinv
 

PatRec()

shipPatRec_prev.PatRec	(		firsttwo,
			zlayer,
			zlayerv2,
			StrawRaw,
			StrawRawLink,
			ReconstructibleMCTracks
	)

Definition at line 694 of file shipPatRec_prev.py.

def PatRec(firsttwo,zlayer,zlayerv2,StrawRaw,StrawRawLink,ReconstructibleMCTracks): 
  global reconstructiblehorizontalidsfound12,reconstructiblestereoidsfound12,reconstructiblehorizontalidsfound34,reconstructiblestereoidsfound34
  global reconstructibleidsfound12,reconstructibleidsfound34,rawhits,totalaftermatching
 
  hits={}
  rawhits={}
  stereohits={}
  hitids={}
  ytan={}
  ycst={}
  horpx={}
  horpy={}
  horpz={}
  stereomom={}
  stereotan={}
  stereocst={}
  fraction={}
  trackid={}
  duplicates=[]
  j=0
  resolution=ShipGeo.strawtubes.sigma_spatial
 
  for item in StrawRaw:  
     #y hits for horizontal straws
     rawhits[j]=copy.deepcopy(((StrawRaw[item][1]+StrawRaw[item][4])/2,(StrawRaw[item][2]+StrawRaw[item][5])/2,StrawRaw[item][6]))
     if firsttwo==True: 
       if global_variables.debug:
     print(
         "rawhits[", j, "]=", rawhits[j],
         "trackid", StrawRawLink[item][0].GetTrackID(),
         "strawname", StrawRawLink[item][0].GetDetectorID(),
         "true x", StrawRawLink[item][0].GetX(),
         "true y", StrawRawLink[item][0].GetY(),
         "true z", StrawRawLink[item][0].GetZ()
         )
     j=j+1    
      
  sortedrawhits=OrderedDict(sorted(rawhits.items(),key=lambda t:t[1][1])) 
  if global_variables.debug: 
     print(" ")
     print("horizontal view (y) hits ordered by plane: plane nr, zlayer, hits")
 
  y2=[]
  y3=[]
  yother=[]
  ymin2=[]
  z2=[]
  z3=[]
  zother=[]
  zmin2=[]
  for i in range(i1,i2+1):
    a=[] #the hits
    c=[] #the keys pointing to the raw hits
    d=[] #the distance to the wire
    for item in sortedrawhits:
      if (float(sortedrawhits[item][1])==float(zlayer[i][0])) : 
        #difference with previous version: make each hit a dictionary so we can later get back the mc trackid
    if sortedrawhits[item][0] not in a: 
       a.append(float(sortedrawhits[item][0]))
    else: 
       #This should never happen - duplicate hits are already removed at digitization
       if global_variables.debug:
         print(" wire hit again", sortedrawhits[item], "strawname=", StrawRawLink[item][0].GetDetectorID())
       a.append(float(sortedrawhits[item][0]))
       duplicates.append(float(zlayer[i][0]))  
    a.sort()   
    #find the key of the sorted hit   
    for value in a:
       for item in sortedrawhits:
          if ((float(sortedrawhits[item][0])==value) & (float(sortedrawhits[item][1])==float(zlayer[i][0]))) :
            if item not in c : 
          c.append(item)
          d.append(float(sortedrawhits[item][2]))
        else : 
           continue
        break   
    #indicate which hit has been used    
    b=len(a)*[0]
    hits[i]=[a,b,c,d]
    if global_variables.debug:
      print(i, zlayer[i], hits[i]) 
    
    # split hits up by trackid for debugging 
 
    for item in range(len(hits[i][0])): 
    if StrawRawLink[hits[i][2][item]][0].GetTrackID() == 3:
        y3.append(float(hits[i][0][item]))
        z3.append(float(zlayer[i][0]))
    else:    
        if StrawRawLink[hits[i][2][item]][0].GetTrackID() == 2:
           y2.append(float(hits[i][0][item]))
           z2.append(float(zlayer[i][0]))
        else:           
           if StrawRawLink[hits[i][2][item]][0].GetTrackID() == -2:
              ymin2.append(float(hits[i][0][item]))
          zmin2.append(float(zlayer[i][0]))
           else:
              yother.append(float(hits[i][0][item]))
          zother.append(float(zlayer[i][0]))                   
 
  if cheated==True :
      if threeprong==1 : nrtrcandv1,trcandv1=ptrack(zlayer,hits,6,0.6)  
      else : nrtrcandv1,trcandv1=ptrack(zlayer,hits,7,0.5) 
  else :
      if threeprong==1 : nrtrcandv1,trcandv1=ptrack(zlayer,hits,6,0.6)  
      else: nrtrcandv1,trcandv1=ptrack(zlayer,hits,7,0.7) 
  
  foundhorizontaltrackids=[]
  horizontaltrackids=[]
  removetrackids=[]
  for t in trcandv1:
    trackids=[]
    if global_variables.debug: 
       print(' ')
       print('track found in Y view:',t,' indices of hits in planes',trcandv1[t])
    for ipl in range(len(trcandv1[t])):      
      indx= trcandv1[t][ipl]
      if indx>-1:     
    if global_variables.debug:
      print(
          '   plane nr, z-position, y of hit, trackid:',
          ipl, zlayer[ipl], hits[ipl][0][indx],
          StrawRawLink[hits[ipl][2][indx]][0].GetTrackID()
      )
    trackids.append(StrawRawLink[hits[ipl][2][indx]][0].GetTrackID())               
    if global_variables.debug:
      print(
      "   Largest fraction of hits in Y view:", fracMCsame(trackids)[0],
      "on MC track with id",fracMCsame(trackids)[1]
      )
    if firsttwo==True: 
       h['fracsame12-y'].Fill(fracMCsame(trackids)[0])
    else: 
       h['fracsame34-y'].Fill(fracMCsame(trackids)[0])
    if monitor==1:   
       if fracMCsame(trackids)[1] in ReconstructibleMCTracks:
          horizontaltrackids.append(fracMCsame(trackids)[1])
          if fracMCsame(trackids)[1] not in foundhorizontaltrackids:
             foundhorizontaltrackids.append(fracMCsame(trackids)[1])
       else:
          #this track is not reconstructible, remove it
          if global_variables.debug:
        print("Y view track with trackid", fracMCsame(trackids)[1], "is not reconstructible. Removing it.")
          removetrackids.append(t) 
 
       if (len(foundhorizontaltrackids) == 0 or len(horizontaltrackids)==0 ): 
          if global_variables.debug:
        print("No reconstructible Y view tracks found.")
          if monitor==True : return 0,[],[],[],[],[],[],{},{},{},{},{}
  
  for i in range(0,len(removetrackids)):
      del trcandv1[removetrackids[i]]
      nrtrcandv1=nrtrcandv1-1
  
  #reorder trcandv1 after removing tracks    
  if len(removetrackids)>0:
     j=0
     for key in trcandv1:
        j=j+1
        if key!=j:
       trcandv1[j]=trcandv1[key]
       del trcandv1[key]
  
  if monitor==1:          
     if len(ReconstructibleMCTracks)>0:    
        if firsttwo==True: 
          if len(foundhorizontaltrackids)>=len(ReconstructibleMCTracks) : reconstructiblehorizontalidsfound12+=1
          else : 
             if global_variables.debug: 
           print("!!!!!!!!!!!!!!!! Difference between Y-view tracks found and reconstructible tracks (station 1&2). Quitting patrec.")
               debugevent(iEvent,False,y2,y3,ymin2,yother,z2,z3,zmin2,zother,foundhorizontaltrackids)     
         return 0,[],[],[],[],[],[],{},{},{},{},{}
        else: 
          if len(foundhorizontaltrackids)>=len(ReconstructibleMCTracks) : reconstructiblehorizontalidsfound34+=1
          else : 
             if global_variables.debug: 
           print("!!!!!!!!!!!!!!!! Difference between Y-view tracks found and reconstructible tracks (station 3&4). Quitting patrec.")
           debugevent(iEvent,False,y2,y3,ymin2,yother,z2,z3,zmin2,zother,foundhorizontaltrackids)
         return 0,[],[],[],[],[],[],{},{},{},{},{}
 
     if len(foundhorizontaltrackids) != reconstructiblerequired :
       if global_variables.debug:
     print(len(foundhorizontaltrackids), "Y view tracks found, but ", reconstructiblerequired, "required.")
       return 0,[],[],[],[],[],[],{},{},{},{},{} 
  else:
     if firsttwo==True: reconstructiblehorizontalidsfound12+=1
     else: reconstructiblehorizontalidsfound34+=1
 
 
  if nrtrcandv1==0 : 
    if global_variables.debug:
      print("0 tracks found in Y view. Reconstructible:", len(ReconstructibleMCTracks))
    if monitor==True: return 0,[],[],[],[],[],[],{},{},{},{},{}
  else : 
    if global_variables.debug:
      print(nrtrcandv1, "tracks found in Y view. Reconstructible:", len(ReconstructibleMCTracks))
   
  if global_variables.debug:
    print("***************** Start of Stereo PatRec **************************")  
 
  v2hits={}
  v2hitsMC={}
 
  uvview={}
  j=0
  rawxhits={}
  for item in StrawRaw:  
     rawxhits[j]=copy.deepcopy((StrawRaw[item][0],StrawRaw[item][1],StrawRaw[item][3],StrawRaw[item][4],StrawRaw[item][2],StrawRaw[item][6]))
     if firsttwo==True: 
       if global_variables.debug:
     print(
         "rawxhits[", j, "]=", rawxhits[j],
         "trackid", StrawRawLink[item][0].GetTrackID(),
         "true x", StrawRawLink[item][0].GetX(),
         "true y",StrawRawLink[item][0].GetY()
         )
     j=j+1  
 
  sortedrawxhits=OrderedDict(sorted(rawxhits.items(),key=lambda t:t[1][4])) 
 
  if global_variables.debug:
    print("stereo view hits ordered by plane:")
  for i in range(i1,i2+1):
    xb=[]
    yb=[]
    xt=[]  
    yt=[]  
    c=[]
    d=[]
    for item in sortedrawxhits:
      if (float(sortedrawxhits[item][4])==float(zlayerv2[i][0])) :
        xt.append(float(sortedrawxhits[item][0]))
        yt.append(float(sortedrawxhits[item][1]))
        xb.append(float(sortedrawxhits[item][2]))
        yb.append(float(sortedrawxhits[item][3]))  
    #c is the rawxhits hit number
    c.append(item)   
    #d is the distance to the wire  
    d.append(float(sortedrawxhits[item][5]))
    uvview[i]=[xb,yb,xt,yt,c,d]
 
    if global_variables.debug:
      print(
      '   uv hits, z,xb,yb,xt,yt,dist    ',
      i, zlayerv2[i], uvview[i][0], uvview[i][1], uvview[i][2],
      uvview[i][3], uvview[i][4], uvview[i][5]
      )
 
  # now do pattern recognition in view perpendicular to first search view
  # loop over tracks found in Y view, and intersect this "plane"
  # with hits in other views, and project in "x", then ptrack in "x", etc..
 
  #loop over tracks found in Y view 
  ntracks=0 
  trackkey=0 
  tracks={}
  if global_variables.debug:
    if (firsttwo==True) :  
      print('Loop over tracks found in Y view, stations 1&2')
    else :  
      print('Loop over tracks found in Y view, stations 3&4')
 
  foundstereotrackids=[]
  foundtrackids=[]
 
  for t in trcandv1:
    alltrackids=[]
    y11trackids=[]
    y12trackids=[]
    y21trackids=[]
    y22trackids=[]
    y11hitids=[]
    y12hitids=[]
    y21hitids=[]
    y22hitids=[]
    trackkey+=1      
    
    #linear "fit" to track found in this view 
 
    fitt,fitc=fitline(trcandv1[t],hits,zlayer,resolution)
    if global_variables.debug:
      print('   Track nr', t, 'in Y view: tangent, constant=', fitt, fitc)
 
    tan=0.
    cst=0.
    px=0.
    py=0.
    pz=0.
    m=0
    if firsttwo==True: looplist=reversed(list(range(len(trcandv1[t]))))
    else: looplist=list(range(len(trcandv1[t])))
    for ipl in looplist:      
      indx= trcandv1[t][ipl]
      if indx>-1: 
    if global_variables.debug:
      print('      Plane nr, z-position, y of hit:', ipl, zlayer[ipl], hits[ipl][0][indx])
        hitpoint=[zlayer[ipl],hits[ipl][0][indx]]
    rc=h['disthittoYviewtrack'].Fill(dist2line(fitt,fitc,hitpoint))
    #if px==0. : px=StrawRawLink[hits[ipl][2][indx]][0].GetPx()
    #if py==0. : py=StrawRawLink[hits[ipl][2][indx]][0].GetPy()
    #if pz==0. : pz=StrawRawLink[hits[ipl][2][indx]][0].GetPz()
    px=StrawRawLink[hits[ipl][2][indx]][0].GetPx()
    py=StrawRawLink[hits[ipl][2][indx]][0].GetPy()
    pz=StrawRawLink[hits[ipl][2][indx]][0].GetPz()
    ptmp=math.sqrt(px**2+py**2+pz**2)
    if global_variables.debug:
        print("      p",ptmp,"px",px,"py",py,"pz",pz)
    if tan==0. : tan=py/pz
    if cst==0. : cst=StrawRawLink[hits[ipl][2][indx]][0].GetY()-tan*zlayer[ipl][0]
    rc=h['disthittoYviewMCtrack'].Fill(dist2line(tan,cst,hitpoint))
 
    if global_variables.debug:
      print('   Track nr', t, 'in Y view: MC tangent, MC constant=', tan, cst)  
    
    for ipl in range(len(trcandv1[t])):      
      indx= trcandv1[t][ipl]
      if indx>-1:   
        diffy=hits[ipl][0][indx]-StrawRawLink[hits[ipl][2][indx]][0].GetY()
    h['digi-truevstruey'].Fill(diffy)
        if ipl<5:   
        y11trackids.append(StrawRawLink[hits[ipl][2][indx]][0].GetTrackID())
        y11hitids.append([hits[ipl][2][indx]][0])
    if (ipl>4 and ipl<9):   
        y12trackids.append(StrawRawLink[hits[ipl][2][indx]][0].GetTrackID())
        y12hitids.append([hits[ipl][2][indx]][0])      
    if (ipl>9 and ipl<13):   
        y21trackids.append(StrawRawLink[hits[ipl][2][indx]][0].GetTrackID())
        y21hitids.append([hits[ipl][2][indx]][0])   
    if ipl>12:   
        y22trackids.append(StrawRawLink[hits[ipl][2][indx]][0].GetTrackID())
        y22hitids.append([hits[ipl][2][indx]][0])
      else: 
        if ipl<5:   
        y11trackids.append(999)
        y11hitids.append(999)
    if (ipl>4 and ipl<9):   
        y12trackids.append(999)
        y12hitids.append(999)
    if (ipl>9 and ipl<13):   
        y21trackids.append(999)
        y21hitids.append(999)
    if ipl>12: 
        y22trackids.append(999)
        y22hitids.append(999)
  
    #intersect this "plane", with all hits in all other views...
 
    v2y2=[]
    v2y3=[]
    v2yother=[]
    v2ymin2=[]
    v2z2=[]
    v2z3=[]
    v2zother=[]
    v2zmin2=[]
    for ipl in range(1,len(zlayerv2)+1):
      items=[]
      xclean=[]
      distances=[]
      if cheated==False:          
     xclean,items=line2plane(fitt,fitc,uvview[ipl],zlayerv2[ipl][0])
 
      else:
         xclean=copy.deepcopy(uvview[ipl][0])      
         items=copy.deepcopy(uvview[ipl][4]) 
     
      distances=copy.deepcopy(uvview[ipl][5])           
      xcleanunsorted=[]
      xcleanunsorted=list(xclean)
 
      xclean.sort()
      b=len(xclean)*[0]
      d=[]
      e=[]
      for item in xclean:  
          d.append(items[xcleanunsorted.index(item)]) 
      e.append(distances[xcleanunsorted.index(item)]) 
      #fill hits info for ptrack, "b" records if hit has been used in ptrack
      v2hits[ipl]=[xclean,b,d,e]
      if global_variables.debug:
    print('      Plane,z, projected hits:', ipl , zlayerv2[ipl], v2hits[ipl])
      
      for item in range(len(v2hits[ipl][0])): 
    if StrawRawLink[v2hits[ipl][2][item]][0].GetTrackID() == 3:
        v2y3.append(float(v2hits[ipl][0][item]))
        v2z3.append(float(zlayerv2[ipl][0]))
    else:    
        if StrawRawLink[v2hits[ipl][2][item]][0].GetTrackID() == 2:
           v2y2.append(float(v2hits[ipl][0][item]))
           v2z2.append(float(zlayerv2[ipl][0]))
        else:           
           if StrawRawLink[v2hits[ipl][2][item]][0].GetTrackID() == -2:
              v2ymin2.append(float(v2hits[ipl][0][item]))
          v2zmin2.append(float(zlayerv2[ipl][0]))
           else:
              v2yother.append(float(v2hits[ipl][0][item]))
          v2zother.append(float(zlayerv2[ipl][0]))
            
    #pattern recognition for this track, blow up y-window, maybe make dependend of stereo angle,
    #constant now, should be 5*error/sin(alpha) :-), i.e. variable/plane if different "alphas"
    #Hence, for "real" y-view, should be small.
    if cheated==True:
       nrtrcandv2,trcandv2=ptrack(zlayerv2,v2hits,7,0.7)
    else:
       nrtrcandv2,trcandv2=ptrack(zlayerv2,v2hits,6,15.)
    if len(trcandv2)>1: 
      if global_variables.debug:
    print("   len(trcandv2) in stereo", len(trcandv2))
    
    nstereotracks=0
    if nrtrcandv2==0 : 
      if global_variables.debug:
    print("   0 tracks found in stereo view, for Y-view track nr", t)
    
 
    #if firsttwo==True: totalstereo12=reconstructiblestereoidsfound12
    #else: totalstereo34==reconstructiblestereoidsfound34
    
    for t1 in trcandv2:
      if global_variables.debug:
    print('   Track belonging to Y-view view track', t, 'found in stereo view:', t1, trcandv2[t1])      
      
      stereofitt,stereofitc=fitline(trcandv2[t1],v2hits,zlayerv2,resolution)
 
      pxMC=0.
      pzMC=0.
      stereotanMCv=0.
      stereocstMCv=0.
      if firsttwo==True: looplist=reversed(list(range(len(trcandv2[t1]))))
      else: looplist=list(range(len(trcandv2[t1])))
      for ipl in looplist:      
        indx= trcandv2[t1][ipl]
        if indx>-1:       
          hitpointx=[zlayerv2[ipl],v2hits[ipl][0][indx]]
          rc=h['disthittostereotrack'].Fill(dist2line(stereofitt,stereofitc,hitpointx)) 
      if pxMC==0. :  pxMC=StrawRawLink[v2hits[ipl][2][indx]][0].GetPx()
      if pzMC ==0. : pzMC=StrawRawLink[v2hits[ipl][2][indx]][0].GetPz()
      if stereotanMCv==0. : stereotanMCv=pxMC/pzMC
      if stereocstMCv==0. : stereocstMCv=StrawRawLink[v2hits[ipl][2][indx]][0].GetX()-stereotanMCv*zlayerv2[ipl][0]
      rc=h['disthittostereoMCtrack'].Fill(dist2line(stereotanMCv,stereocstMCv,hitpointx))
                 
 
      stereotrackids=[]
 
      stereo11trackids=[]
      stereo12trackids=[]
      stereo21trackids=[]
      stereo22trackids=[]
      stereo11hitids=[]
      stereo12hitids=[]
      stereo21hitids=[]
      stereo22hitids=[]
      nstereotracks+=1
     
 
      for ipl in range(len(trcandv2[t1])):      
        indx= trcandv2[t1][ipl]
        if indx>-1: 
      stereotrackids.append(StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID())     
      if global_variables.debug:
        print(
        "      plane nr, zpos, x of hit, hitid, trackid",
        ipl, zlayerv2[ipl], v2hits[ipl][0][indx], v2hits[ipl][2][indx],
        StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID()
        )   
      xdiff=v2hits[ipl][0][indx]-StrawRawLink[v2hits[ipl][2][indx]][0].GetX()
      h['digi-truevstruex'].Fill(xdiff)
    if ipl<5:   
             if indx>-1: 
            stereo11trackids.append(StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID())
        stereo11hitids.append(v2hits[ipl][2][indx]) 
         else: 
            stereo11trackids.append(999)
            stereo11hitids.append(999)      
    if (ipl>4 and ipl<9):   
         if indx>-1: 
            stereo12trackids.append(StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID())
        stereo12hitids.append(v2hits[ipl][2][indx])         
         else: 
            stereo12trackids.append(999)
            stereo12hitids.append(999)      
    if (ipl>9 and ipl<13):   
         if indx>-1: 
            stereo21trackids.append(StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID())
        stereo21hitids.append(v2hits[ipl][2][indx])         
         else: 
            stereo21trackids.append(999)
            stereo21hitids.append(999)      
    if ipl>12:   
         if indx>-1: 
            stereo22trackids.append(StrawRawLink[v2hits[ipl][2][indx]][0].GetTrackID())  
        stereo22hitids.append(v2hits[ipl][2][indx])                    
             else: 
            stereo22trackids.append(999)
            stereo22hitids.append(999)      
      
      if firsttwo==True: h['fracsame12-stereo'].Fill(fracMCsame(stereotrackids)[0])
      else: h['fracsame34-stereo'].Fill(fracMCsame(stereotrackids)[0])
      
      if global_variables.debug: 
         print("      Largest fraction of hits in stereo view:",fracMCsame(stereotrackids)[0],"on MC track with id",fracMCsame(stereotrackids)[1])    
         #check if this trackid is the same as the track from the Y-view it belongs to
         if monitor==True: print("      fracMCsame(stereotrackids)[1]", fracMCsame(stereotrackids)[1],"foundhorizontaltrackids[",t-1,"]",foundhorizontaltrackids[t-1])
      
      if monitor==True: 
         if fracMCsame(stereotrackids)[1] != horizontaltrackids[t-1] :
            if global_variables.debug:
          print(
          "      Stereo track with trackid",
          fracMCsame(stereotrackids)[1],
          "does not belong to the Y-view track with id=",
          horizontaltrackids[t-1]
          )
            continue
      
         #check if this trackid belongs to a reconstructible track
         if fracMCsame(stereotrackids)[1] in ReconstructibleMCTracks:         
            if fracMCsame(stereotrackids)[1] not in foundstereotrackids: 
           foundstereotrackids.append(fracMCsame(stereotrackids)[1])
         else:
            if global_variables.debug:
          print(
          "      Stereo track with trackid",
          fracMCsame(stereotrackids)[1], 
          "is not reconstructible. Removing it."
          )
        continue
      else:
          if fracMCsame(stereotrackids)[1] not in foundstereotrackids: 
         foundstereotrackids.append(fracMCsame(stereotrackids)[1])   
     
      ntracks=ntracks+1     
      alltrackids=y11trackids+stereo11trackids+stereo12trackids+y12trackids+y21trackids+stereo21trackids+stereo22trackids+y22trackids 
      if firsttwo==True: h['fracsame12'].Fill(fracMCsame(alltrackids)[0])
      else: h['fracsame34'].Fill(fracMCsame(alltrackids)[0])      
      if global_variables.debug:
    print(
        "      Largest fraction of hits in horizontal and stereo view:", fracMCsame(alltrackids)[0],
        "on MC track with id", fracMCsame(alltrackids)[1]
        )
      
      #calculate efficiency after joining horizontal and stereo tracks
      if monitor==True: 
         if fracMCsame(alltrackids)[1] in ReconstructibleMCTracks:
       if fracMCsame(alltrackids)[1] not in foundtrackids:
         foundtrackids.append(fracMCsame(alltrackids)[1])
                     
         tracks[trackkey*1000+nstereotracks]=alltrackids
         hitids[trackkey*1000+nstereotracks]=y11hitids+stereo11hitids+stereo12hitids+y12hitids+y21hitids+stereo21hitids+stereo22hitids+y22hitids    
         if cheated==False :
               ytan[trackkey*1000+nstereotracks]=fitt
               ycst[trackkey*1000+nstereotracks]=fitc
         else:
               ytan[trackkey*1000+nstereotracks]=tan
               ycst[trackkey*1000+nstereotracks]=cst      
           
    
         fraction[trackkey*1000+nstereotracks]=fracMCsame(alltrackids)[0]
         trackid[trackkey*1000+nstereotracks]=fracMCsame(alltrackids)[1]
 
  
         horpx[trackkey*1000+nstereotracks]=px
         horpy[trackkey*1000+nstereotracks]=py
         horpz[trackkey*1000+nstereotracks]=pz
                
         if cheated==False:
               stereotan[trackkey*1000+nstereotracks]=stereofitt
               stereocst[trackkey*1000+nstereotracks]=stereofitc 
         else:
               stereotan[trackkey*1000+nstereotracks]=stereotanMCv
               stereocst[trackkey*1000+nstereotracks]=stereocstMCv 
         else:
            if global_variables.debug:
          print("Track with trackid", fracMCsame(alltrackids)[1], "is not reconstructible. Removing it.")
        continue
      else:
       if fracMCsame(alltrackids)[1] not in foundtrackids:
         foundtrackids.append(fracMCsame(alltrackids)[1])
                     
         tracks[trackkey*1000+nstereotracks]=alltrackids
         hitids[trackkey*1000+nstereotracks]=y11hitids+stereo11hitids+stereo12hitids+y12hitids+y21hitids+stereo21hitids+stereo22hitids+y22hitids    
         if cheated==False :
               ytan[trackkey*1000+nstereotracks]=fitt
               ycst[trackkey*1000+nstereotracks]=fitc
         else:
               ytan[trackkey*1000+nstereotracks]=tan
               ycst[trackkey*1000+nstereotracks]=cst      
           
    
         fraction[trackkey*1000+nstereotracks]=fracMCsame(alltrackids)[0]
         trackid[trackkey*1000+nstereotracks]=fracMCsame(alltrackids)[1]
 
  
         horpx[trackkey*1000+nstereotracks]=px
         horpy[trackkey*1000+nstereotracks]=py
         horpz[trackkey*1000+nstereotracks]=pz
                
         if cheated==False:
               stereotan[trackkey*1000+nstereotracks]=stereofitt
               stereocst[trackkey*1000+nstereotracks]=stereofitc 
         else:
               stereotan[trackkey*1000+nstereotracks]=stereotanMCv
               stereocst[trackkey*1000+nstereotracks]=stereocstMCv
            
  if monitor==True: 
    if len(foundstereotrackids)>=len(ReconstructibleMCTracks) :
      if firsttwo==True: reconstructiblestereoidsfound12+=1
      else: reconstructiblestereoidsfound34+=1    
    else:     
      if global_variables.debug:
         debugevent(iEvent,True,v2y2,v2y3,v2ymin2,v2yother,v2z2,v2z3,v2zmin2,v2zother,foundstereotrackids)
         print("Nbr of reconstructible tracks after stereo",len(foundstereotrackids)," but ",len(ReconstructibleMCTracks)," reconstructible tracks in this event. Quitting.")
      return 0,[],[],[],[],[],[],{},{},{},{},{}
     
    if len(foundtrackids)>=len(ReconstructibleMCTracks):
      if firsttwo==True: reconstructibleidsfound12+=1
      else: reconstructibleidsfound34+=1  
    else: 
      if global_variables.debug: 
         debugevent(iEvent,True,v2y2,v2y3,v2ymin2,v2yother,v2z2,v2z3,v2zmin2,v2zother,foundstereotrackids)
         print("Nbr of reconstructed tracks ",len(foundtrackids)," but",len(ReconstructibleMCTracks)," reconstructible tracks in this event. Quitting.")
      return 0,[],[],[],[],[],[],{},{},{},{},{}
  else:
    if firsttwo==True: reconstructiblestereoidsfound12+=1
    else: reconstructiblestereoidsfound34+=1  
    if firsttwo==True: reconstructibleidsfound12+=1
    else: reconstructibleidsfound34+=1      
 
  #now Kalman fit, collect "all hits" around fitted track, outlier removal etc..  
  return ntracks,tracks,hitids,ytan,ycst,stereotan,stereocst,horpx,horpy,horpz,fraction,trackid
 

ptrack()

shipPatRec_prev.ptrack	(		zlayer,
			ptrackhits,
			nrwant,
			window
	)

Definition at line 1389 of file shipPatRec_prev.py.

def ptrack(zlayer,ptrackhits,nrwant,window):
 
# basic pattern recognition in one view
# zlayer= dictionary with z-position of planes: zlayer[iplane][0]==z-position (list of length 1 :-).
# ptrackhits is dictionary with hits in a projection:  ptrackhits[plane number][0]=list-of-hits
# ndrop= nr of hits that can be missing in this projection
# nrwant= min nr of hits on a track in this projection
# window: is the size of the search window
 
# get first and last plane number (needs to be consecutive, and also contains empty planes
# should be included in the list!
 planes=list(zlayer.keys())
 planes.sort()
 i_1=planes[0]
 i_2=planes[len(planes)-1]
 ndrop=len(planes)-nrwant
 #print 'ptrack input: planes=',i_1,i_2,ndrop,window,"   ptrackhits ",ptrackhits
 
 nrtracks=0
 tracks={}
 
 #loop over maxnr hits, to max-ndrop
 for nhitw in range(i_2-i_1+1,i_2-i_1-ndrop,-1):
  # nhitw: wanted number of hits when looking for a track
  for idrop in range(ndrop):
    #only start if wanted nr hits <= the nr of planes 
    nrhitmax=i_2-i_1-idrop+1
    if nhitw<=nrhitmax:
     for k in range(idrop+1):
        #calculate the id of the first and last plane for this try.
        ifirst=i_1+k
        ilast=i_2-(idrop-k)
        #now loop over hits in first/last planes, and construct line
        dz=zlayer[ilast][0]-zlayer[ifirst][0]
        #hits in first plane
        for ifr in range(len(ptrackhits[ifirst][0])):
      #for ifr in range(len(ptrackhits[ifirst][0])):
          #has this hit been used already for a track: skip
      if  ptrackhits[ifirst][1][ifr]==0:
           xfirst= ptrackhits[ifirst][0][ifr]
       #hits in last plane
           for il in range(len(ptrackhits[ilast][0])):
            #has this hit been used already for a track: skip
            if  ptrackhits[ilast][1][il]==0:         
             xlast= ptrackhits[ilast][0][il]
             nrhitsfound=2
             tancand=(xlast-xfirst)/dz
             #fill temporary hit list for track-candidate with -1
             trcand=(i_2-i_1+2)*[-1]
             #loop over in between planes
             for inbetween in range(ifirst+1,ilast):
              #can wanted number of hits be satisfied?
              if nrhitsfound+ilast-inbetween>=nhitw:
               #calculate candidate hit
               xwinlow=xfirst+tancand*(zlayer[inbetween][0]-zlayer[ifirst][0])-window
               for im in range(len(ptrackhits[inbetween][0])):
                #has this hit been used?
                if  ptrackhits[inbetween][1][im]==0:
                 xin= ptrackhits[inbetween][0][im]
                 #hit belwo window, try next one
                 if xin>xwinlow:
                  #if hit larger than upper edge of window: goto next plane
                  if xin>xwinlow+2*window: 
             break
                  #hit found, do administation
                  if trcand[inbetween]==-1: 
            nrhitsfound+=1
                  trcand[inbetween]=im
             #looped over in between planes, collected enough hits?
             if nrhitsfound==nhitw:
               #mark used hits
               trcand[ifirst]=ifr
               trcand[ilast]=il
               for i in range(ifirst,ilast+1):
                  if trcand[i]>=0: 
             ptrackhits[i][1][trcand[i]]=1
               #store the track
               nrtracks+=1
           #print "ptrack nrtracks",nrtracks
               tracks[nrtracks]=trcand
 #print "ptrack: tracks",tracks        
 return nrtracks,tracks  
 

SmearHits()

shipPatRec_prev.SmearHits	(		iEvent,
			sTree,
			modules,
			SmearedHits,
			ReconstructibleMCTracks
	)

Definition at line 495 of file shipPatRec_prev.py.

def SmearHits(iEvent,sTree,modules,SmearedHits,ReconstructibleMCTracks):
  #smears hits (when not cheated) 
  #apply cuts for >500 hits, duplicate straw hits and acceptance
  #call this routine once for each event, before the digitization
  
  top=ROOT.TVector3()
  bot=ROOT.TVector3()
  random = ROOT.TRandom()
  ROOT.gRandom.SetSeed(13)
  
  rc = sTree.GetEvent(iEvent) 
  nHits = sTree.strawtubesPoint.GetEntriesFast()
  withNoStrawSmearing=None
  hitstraws={}
  duplicatestrawhit=[]
 
  for i in range(nHits):
    ahit = sTree.strawtubesPoint[i]
    if (str(ahit.GetDetectorID())[:1]=="5") : continue
    strawname=str(ahit.GetDetectorID())    
    if strawname in hitstraws:
       #straw was already hit
       if ahit.GetX()>hitstraws[strawname][1]:
          #this hit has higher x, discard it
          duplicatestrawhit.append(i)
       else:
          #del hitstraws[strawname] 
      duplicatestrawhit.append(hitstraws[strawname][0])
          hitstraws[strawname]=[i,ahit.GetX()]          
    else: 
       hitstraws[strawname]=[i,ahit.GetX()]
          
  #the following code prints some histograms related to the MC hits  
  strawname=''
  if global_variables.debug:
    print("nbr of hits=", nHits, "in event", iEvent)
  station1hits={}
  station12xhits={}
  station12yhits={}
  station2hits={}
  station3hits={}
  station4hits={}  
  
  for i in range(nHits):
    ahit = sTree.strawtubesPoint[i]
    strawname=str(ahit.GetDetectorID())
    #is it a duplicate hit? if so, ignore it       
    if i in duplicatestrawhit: 
      continue
    #only look at hits in the strawtracking stations
    if (str(ahit.GetDetectorID())[:1]=="5") : continue
    #is hit inside acceptance?  
    if (((ahit.GetX()/245.)**2 + (ahit.GetY()/495.)**2) >= 1.): continue
    
    modules["Strawtubes"].StrawEndPoints(ahit.GetDetectorID(),bot,top)    
    if cheated==False : 
       #this is where the smearing takes place. the hit coordinates can also be read from somewhere else.
       sm   = hit2wire(ahit,bot,top,withNoStrawSmearing)
       m = array('d',[i,sm['xtop'],sm['ytop'],sm['z'],sm['xbot'],sm['ybot'],sm['z'],sm['dist'],ahit.GetDetectorID()])
    else :   
       #for MC truth 
       m = array('d',[i,ahit.GetX(),ahit.GetY(),top[2],ahit.GetX(),ahit.GetY(),top[2],ahit.dist2Wire(),ahit.GetDetectorID()])
    
    measurement = ROOT.TVectorD(9,m)
 
    smHits = SmearedHits.GetEntries()
    if SmearedHits.GetSize() == smHits: SmearedHits.Expand(smHits+1000)
    SmearedHits[smHits] = measurement
 
   
    angle=0.    
    if (str(ahit.GetDetectorID())[1:2]=="1"): angle=ShipGeo.strawtubes.ViewAngle
    if (str(ahit.GetDetectorID())[1:2]=="2"): angle=ShipGeo.strawtubes.ViewAngle*-1.
    if (str(ahit.GetDetectorID())[:1]=="1") :  
       if ahit.GetTrackID() in station1hits:
          station1hits[ahit.GetTrackID()]+=1
      rc=h['hits1xy'].Fill(ahit.GetX(),ahit.GetY())       
       else:
          station1hits[ahit.GetTrackID()]=1
    if (str(ahit.GetDetectorID())[:1]=="2") :  
       if ahit.GetTrackID() in station2hits:
          station2hits[ahit.GetTrackID()]+=1
       else:
          station2hits[ahit.GetTrackID()]=1
    if (str(ahit.GetDetectorID())[:1]=="3") :  
       if ahit.GetTrackID() in station3hits:
          station3hits[ahit.GetTrackID()]+=1
       else:
          station3hits[ahit.GetTrackID()]=1
    if (str(ahit.GetDetectorID())[:1]=="4") :   
       if ahit.GetTrackID() in station4hits:
          station4hits[ahit.GetTrackID()]+=1
       else:
          station4hits[ahit.GetTrackID()]=1
      
    if ((str(ahit.GetDetectorID())[:2]=="11" or str(ahit.GetDetectorID())[:2]=="12") or (str(ahit.GetDetectorID())[:2]=="21" or str(ahit.GetDetectorID())[:2]=="22")):
       if ahit.GetTrackID() in station12xhits:
      station12xhits[ahit.GetTrackID()]+=1        
       else:
          station12xhits[ahit.GetTrackID()]=1
      
    if ((str(ahit.GetDetectorID())[:2]=="10" or str(ahit.GetDetectorID())[:2]=="13") or (str(ahit.GetDetectorID())[:2]=="20" or str(ahit.GetDetectorID())[:2]=="23")):
       if ahit.GetTrackID() in station12yhits:
      station12yhits[ahit.GetTrackID()]+=1        
       else:
          station12yhits[ahit.GetTrackID()]=1  
        
  total1hits=0 
  total12xhits=0 
  total12yhits=0  
  hits1pertrack=0
  hits12xpertrack=0
  hits12ypertrack=0  
  total2hits=0  
  hits2pertrack=0
  total3hits=0  
  hits3pertrack=0
  total4hits=0  
  hits4pertrack=0
  for items in station1hits:
     # only count the hits on reconstructible tracks
     if monitor==True:
       if items in ReconstructibleMCTracks: total1hits=total1hits+station1hits[items]     
     else : total1hits=total1hits+station1hits[items]
  if len(station1hits) > 0 : 
     hits1pertrack=total1hits/len(station1hits)
  for items in station12xhits:
     if monitor==True:
        if items in ReconstructibleMCTracks: total12xhits=total12xhits+station12xhits[items]     
     else: total12xhits=total12xhits+station12xhits[items]
  if len(station12xhits) > 0 : hits12xpertrack=total12xhits/len(station12xhits)
  for items in station12yhits:
     if monitor==True:
        if items in ReconstructibleMCTracks: total12yhits=total12yhits+station12yhits[items]        
     else: total12yhits=total12yhits+station12yhits[items]
  if len(station12yhits) > 0 : hits12ypertrack=total12yhits/len(station12yhits)
  for items in station2hits:
     if monitor==True:
        if items in ReconstructibleMCTracks: total2hits=total2hits+station2hits[items]
     else: total2hits=total2hits+station2hits[items]
  if len(station2hits) > 0 : 
     hits2pertrack=total2hits/len(station2hits)
  for items in station3hits:
     if monitor==True:
        if items in ReconstructibleMCTracks: total3hits=total3hits+station3hits[items]     
     else: total3hits=total3hits+station3hits[items]
  if len(station3hits) > 0 : 
     hits3pertrack=total3hits/len(station3hits)
  for items in station4hits:
     if monitor==True:
        if items in ReconstructibleMCTracks: total4hits=total4hits+station4hits[items]  
     else:  total4hits=total4hits+station4hits[items]
  if len(station4hits) > 0 : 
     hits4pertrack=total4hits/len(station4hits)  
  
  rc=h['hits1-4'].Fill(hits1pertrack+hits2pertrack+hits3pertrack+hits4pertrack)  
  rc=h['hits1'].Fill(hits1pertrack)  
  rc=h['hits12x'].Fill(hits12xpertrack)  
  rc=h['hits12y'].Fill(hits12ypertrack)  
  rc=h['hits2'].Fill(hits2pertrack)    
  rc=h['hits3'].Fill(hits3pertrack)    
  rc=h['hits4'].Fill(hits4pertrack)        
  return SmearedHits
 
  

TrackFit()

shipPatRec_prev.TrackFit	(		hitPosList,
			theTrack,
			charge,
			pinv
	)

Definition at line 1325 of file shipPatRec_prev.py.

def TrackFit(hitPosList,theTrack,charge,pinv): 
   global theTracks
   if global_variables.debug:
     fitter.setDebugLvl(1)
   resolution = ShipGeo.strawtubes.sigma_spatial    
   hitCov = ROOT.TMatrixDSym(7)
   hitCov[6][6] = resolution*resolution
   for item in hitPosList:
     itemarray=array('d',[item[0],item[1],item[2],item[3],item[4],item[5],item[6]])
     ms=ROOT.TVectorD(7,itemarray) 
     tp = ROOT.genfit.TrackPoint(theTrack) # note how the point is told which track it belongs to 
     measurement = ROOT.genfit.WireMeasurement(ms,hitCov,1,6,tp) # the measurement is told which trackpoint it belongs to
     measurement.setMaxDistance(0.5*u.cm)     
     tp.addRawMeasurement(measurement) # package measurement in the TrackPoint                                          
     theTrack.insertPoint(tp)  # add point to Track
   theTracks.append(theTrack)
   if not global_variables.debug:
     return # leave track fitting to shipDigiReco
#check
   if not theTrack.checkConsistency():
     print('Problem with track before fit, not consistent', theTrack)
     return
# do the fit
   try:  fitter.processTrack(theTrack)
   except: 
     print("genfit failed to fit track")
     return
#check
   if not theTrack.checkConsistency():
     print('Problem with track after fit, not consistent', theTrack)
     return  
    
      
   fitStatus   = theTrack.getFitStatus()
   theTrack.prune("CFL")  #  http://sourceforge.net/p/genfit/code/HEAD/tree/trunk/core/include/Track.h#l280 
 
   nmeas = fitStatus.getNdf()
   pval = fitStatus.getPVal()
   
   #pval close to 0 indicates a bad fit
   chi2        = fitStatus.getChi2()/nmeas
 
   rc=h['chi2fittedtracks'].Fill(chi2)
   rc=h['pvalfittedtracks'].Fill(pval) 
 
   
   fittedState = theTrack.getFittedState()
   fittedMom = fittedState.getMomMag()  
   fittedMom = fittedMom*int(charge) 
   
   if math.fabs(pinv) > 0.0 : rc=h['pvspfitted'].Fill(1./pinv,fittedMom) 
   fittedtrackDir = fittedState.getDir()
   fittedx=math.degrees(math.acos(fittedtrackDir[0]))
   fittedy=math.degrees(math.acos(fittedtrackDir[1]))
   fittedz=math.degrees(math.acos(fittedtrackDir[2]))      
   fittedmass = fittedState.getMass()
   rc=h['momentumfittedtracks'].Fill(fittedMom)
   rc=h['xdirectionfittedtracks'].Fill(fittedx)
   rc=h['ydirectionfittedtracks'].Fill(fittedy)
   rc=h['zdirectionfittedtracks'].Fill(fittedz)
   rc=h['massfittedtracks'].Fill(fittedmass)   
      
   return
 

Variable Documentation

cheated

int shipPatRec_prev.cheated = 0

Definition at line 36 of file shipPatRec_prev.py.

falsenegative

int shipPatRec_prev.falsenegative = 0

Definition at line 32 of file shipPatRec_prev.py.

falsepositive

int shipPatRec_prev.falsepositive = 0

Definition at line 33 of file shipPatRec_prev.py.

fgeo

str shipPatRec_prev.fgeo = ''

Definition at line 42 of file shipPatRec_prev.py.

fittedtrackids

list shipPatRec_prev.fittedtrackids = []

Definition at line 48 of file shipPatRec_prev.py.

fitter

shipPatRec_prev.fitter = ROOT.genfit.DAF()

Definition at line 55 of file shipPatRec_prev.py.

geoFile

str shipPatRec_prev.geoFile = ''

Definition at line 41 of file shipPatRec_prev.py.

h

dict shipPatRec_prev.h = {}

Definition at line 57 of file shipPatRec_prev.py.

i1

int shipPatRec_prev.i1 = 1

Definition at line 115 of file shipPatRec_prev.py.

i2

int shipPatRec_prev.i2 = 16

Definition at line 116 of file shipPatRec_prev.py.

MatchedReconstructibleMCTracks

list shipPatRec_prev.MatchedReconstructibleMCTracks = []

Definition at line 47 of file shipPatRec_prev.py.

monitor

int shipPatRec_prev.monitor = 0

Definition at line 37 of file shipPatRec_prev.py.

morethan100tracks

int shipPatRec_prev.morethan100tracks = 0

Definition at line 31 of file shipPatRec_prev.py.

morethan500

int shipPatRec_prev.morethan500 = 0

Definition at line 30 of file shipPatRec_prev.py.

particles

list shipPatRec_prev.particles = ["e-","e+","mu-","mu+","pi-","pi+","other"]

Definition at line 102 of file shipPatRec_prev.py.

PDG

shipPatRec_prev.PDG = ROOT.TDatabasePDG.Instance()

Definition at line 54 of file shipPatRec_prev.py.

printhelp

int shipPatRec_prev.printhelp = 0

Definition at line 38 of file shipPatRec_prev.py.

random

shipPatRec_prev.random = ROOT.TRandom()

Definition at line 51 of file shipPatRec_prev.py.

rc

dict shipPatRec_prev.rc = h['pinvvstruepinv'].SetMarkerStyle(8)

Definition at line 99 of file shipPatRec_prev.py.

reconstructibleevents

int shipPatRec_prev.reconstructibleevents = 0

Definition at line 34 of file shipPatRec_prev.py.

reconstructiblehorizontalidsfound12

int shipPatRec_prev.reconstructiblehorizontalidsfound12 = 0

Definition at line 24 of file shipPatRec_prev.py.

reconstructiblehorizontalidsfound34

int shipPatRec_prev.reconstructiblehorizontalidsfound34 = 0

Definition at line 26 of file shipPatRec_prev.py.

reconstructibleidsfound12

int shipPatRec_prev.reconstructibleidsfound12 = 0

Definition at line 28 of file shipPatRec_prev.py.

reconstructibleidsfound34

int shipPatRec_prev.reconstructibleidsfound34 = 0

Definition at line 29 of file shipPatRec_prev.py.

ReconstructibleMCTracks

list shipPatRec_prev.ReconstructibleMCTracks = []

Definition at line 46 of file shipPatRec_prev.py.

reconstructiblerequired

int shipPatRec_prev.reconstructiblerequired = 2

Definition at line 39 of file shipPatRec_prev.py.

reconstructiblestereoidsfound12

int shipPatRec_prev.reconstructiblestereoidsfound12 = 0

Definition at line 25 of file shipPatRec_prev.py.

reconstructiblestereoidsfound34

int shipPatRec_prev.reconstructiblestereoidsfound34 = 0

Definition at line 27 of file shipPatRec_prev.py.

theTracks

list shipPatRec_prev.theTracks = []

Definition at line 49 of file shipPatRec_prev.py.

threeprong

int shipPatRec_prev.threeprong = 0

Definition at line 40 of file shipPatRec_prev.py.

totalaftermatching

int shipPatRec_prev.totalaftermatching = 0

Definition at line 44 of file shipPatRec_prev.py.

totalafterpatrec

int shipPatRec_prev.totalafterpatrec = 0

Definition at line 45 of file shipPatRec_prev.py.

TStation1StartZ

int shipPatRec_prev.TStation1StartZ = 0.

Definition at line 121 of file shipPatRec_prev.py.

TStation4EndZ

int shipPatRec_prev.TStation4EndZ = 0.

Definition at line 122 of file shipPatRec_prev.py.

VetoStationEndZ

int shipPatRec_prev.VetoStationEndZ = 0.

Definition at line 124 of file shipPatRec_prev.py.

VetoStationZ

int shipPatRec_prev.VetoStationZ = 0.

Definition at line 123 of file shipPatRec_prev.py.

z34layer

dict shipPatRec_prev.z34layer = {}

Definition at line 119 of file shipPatRec_prev.py.

z34layerv2

dict shipPatRec_prev.z34layerv2 = {}

Definition at line 120 of file shipPatRec_prev.py.

zlayer

dict shipPatRec_prev.zlayer = {}

Definition at line 117 of file shipPatRec_prev.py.

zlayerv2

dict shipPatRec_prev.zlayerv2 = {}

Definition at line 118 of file shipPatRec_prev.py.