Hello,
I am trying to optimise the speed to a PyRoot script that combines the invariant masses of muon tracks in a TTree and then combines the resulting invariant mass with pairwise combinations of oppositely charged pion tracks. I have written a function that performs the task at hand, however, it is crucial for me to try to decrease the runtime.
Here is my current attempt:
from ROOT import *
from math import *
import heapq
import numpy as np
import itertools
from array import array
def MM(tree, hist)::
    for entry in tree:
        #access the momenta and other vars
        buf_run = np.array(entry.run)
        buf_evt = np.array(entry.evt)
        buf_bunchID = np.array(entry.bunchid)
        bufpx = np.array(entry.px) 
        bufpy = np.array(entry.py)
        bufpz = np.array(entry.pz)
        bufpt = np.array(entry.pt)
        bufPIDp = np.array(entry.pidp)
        bufPIDk = np.array(entry.pidk)
        bufPIDpi = np.array(entry.pidpi)
        isMuon = np.array(entry.muon)
        longTrack = np.array(entry.nlong)
        Charge = np.array(entry.charge)
        
        
                      
        if sum(isMuon)>=1 and sum(Charge[np.where(isMuon==1)])==0: #at least one muon track and sum of charges of muon tracks = 0
            
            #create an empty 4-vector to track muon mass
            nul_mu = TLorentzVector()
            nul_mu.SetXYZM(0.,0.,0.,0.)
             
            #for every track in the event 
            for i in range(len(bufpx)):
                
                #add muon 4-vecs
                if isMuon[i]==1:
                  mubuf_particle = TLorentzVector()
                  mubuf_particle.SetXYZM(bufpx[i], bufpy[i], bufpz[i], m_mu)
                  nul_mu+=mubuf_particle
            
            #reconstruct mass around the J-Psi nominal mass value and at least 2 other non-muon (pion) tracks
            if nul_mu.M() > (m_JPsi - 65.)  and nul_mu.M() < (m_JPsi + 65.) and len(np.where(isMuon!=1)[0])>=2: 
                
                pi_idxs=np.where(isMuon!=1) #indices of pio tracks
                temp=np.array(list(itertools.combinations(pi_idxs[0], 2))) #combinatorics of indices of hadron tracks, no repetitions, pairwise combinations
                neutral_pairs_idxs = temp[np.where(np.sum(Charge[temp], axis=1)==0)] 
                  
                
                for pair in neutral_pairs_idxs:
                    idx1 = pair[0]
                    idx2 = pair[1]
                    buf_piplus = TLorentzVector()
                    buf_piminus = TLorentzVector()
                    buf_piplus.SetXYZM( bufpx[idx1], bufpy[idx1], bufpz[idx1] , m_pi ) 
                    buf_piminus.SetXYZM( bufpx[idx2], bufpy[idx2], bufpz[idx2] , m_pi ) 
                   
                    #Fill with JPsi and 2 hadron (opposite charge) tracks
                    #hist.Fill( (nul_mu+buf_piplus+buf_piminus).M())
                    
                    #Fill mu+mu+pi+pi_MM - mu+mu_MM + nominal mass value of J/Psi (=3097), should improve S/B separation
                    #mumu_array[0]=(nul_mu+buf_piplus+buf_piminus).M() - nul_mu.M() + 3097
                    #tdata.Fill()
                    hist.Fill( (nul_mu+buf_piplus+buf_piminus).M() - nul_mu.M() + 3097)
In my attempt to speed to the code I use only the branches I am interested in
#worry only about branches of interest, in the hope to speed the code
c.SetBranchStatus("*", 0)
c.SetBranchStatus("px",1)
c.SetBranchStatus("py",1)
c.SetBranchStatus("pz",1)
c.SetBranchStatus("muon",1)
c.SetBranchStatus("charge",1)
Can anyone suggest a way to speed up the code either in reading the entries or perform the pairwise combinations?
I have been thinking of writing a Root script to perform this task, and I do not seem to manage to perform the pairwise combinations successfully. Can anyone suggest a way to do this?
Thanks,
bldelane
