# # To be used with pyRoot # # This macro demonstrates a simulation of pi0 -> gamma gamma decays # within ROOT. # For given pi0 energy with momentum along Z the decay is isotrop in the # center of mass system (cms). The cms decay angle thetaCM is drawn from # [0,pi]. No phi dependence is simulated. The photon energies in the # laboratory system are calculated using thetaCM, beta and gamma. # Units: c = 1 , Energy in GeV # The generated quantities E and theta are smeared according to the # resolution of a measurement. # WS 2016 J.M. # import ROOT import math import time from array import array InputFileName = "pi0DeacyCaloRooTree.root" ; Pi = 3.1416 # number of events nEvent = 1000 # setup, choose the energy of the pi0, units are GeV pi0Energy = 2. # Pi0 Mass pi0Mass = 0.1349766 # Calorimeter resolution parameters EnergyPar = array( 'd',[ 0.05,0.07,0.01]) # Angular resolution parameters ThetaPar = 0.004 # Get relativistic quantities pi0Mass2 = pi0Mass * pi0Mass beta = math.sqrt(pi0Energy*pi0Energy - pi0Mass*pi0Mass)/pi0Energy gamma = pi0Energy / pi0Mass print("Pi0 energy is set to ", pi0Energy) print("beta = ",beta," gamma= ",gamma ) # minimum opening angle alphaMin = 2. * math.asin(1./gamma) # minimum and maximum values of the energy EMax = gamma * pi0Mass/2. * (1. + beta ) EMin = gamma * pi0Mass/2. * (1. - beta ) print("Minimum opening angle : ",alphaMin) print("Maximum photon energy : ",EMax) print("Minimum photon energy : ",EMin) print("Generate ",nEvent," pi0 decays") # open root file which holds the generated data f = ROOT.TFile.Open(InputFileName,"RECREATE") # Initialize random number generator R = ROOT.TRandom3() # R.SetSeed(7892344) # fixed seed for testing R.SetSeed() # set seed to machine clock # need to set tree variables as array elements thetaCM = array('d',[0]) EPi0 = array('d',[pi0Energy]) E1 = array('d',[0]) E2 = array('d',[0]) Theta1 = array('d',[0]) Theta2 = array('d',[0]) alpha = array('d',[0]) mpi0Calc = array('d',[0]) E1Measure = array('d',[0]) E2Measure = array('d',[0]) Theta1Measure = array('d',[0]) Theta2Measure = array('d',[0]) alphaMeasure = array('d',[0]) mpi0CalcMeasure = array('d',[0]) photon1 = ROOT.TLorentzVector() photon2 = ROOT.TLorentzVector() photon1Measure = ROOT.TLorentzVector() photon2Measure = ROOT.TLorentzVector() # TTree Object tree = ROOT.TTree("pi0GenCalo","pi0 decay tree") tree.Branch("EPi0",EPi0,'EPi0/D') tree.Branch("thetaCM",thetaCM,'thetaCM/D') tree.Branch("E1",E1,"E1/D"); tree.Branch("E2",E2,"E2/D"); tree.Branch("Theta1",Theta1,"Theta1/D") tree.Branch("Theta2",Theta2,"Theta2/D") tree.Branch("alpha",alpha,"alpha/D") tree.Branch("mpi0Calc",mpi0Calc,"mpi0Calc/D") tree.Branch("E1Measure",E1Measure,"E1Measure/D") tree.Branch("E2Measure",E2Measure,"E2Measure/D") tree.Branch("Theta1Measure",Theta1Measure,"Theta1Measure/D") tree.Branch("Theta2Measure",Theta2Measure,"Theta2Measure/D") tree.Branch("alphaMeasure",alphaMeasure,"alphaMeasure/D") tree.Branch("mpi0CalcMeasure",mpi0CalcMeasure,"mpi0CalcMeasure/D") tree.Branch("photon1",photon1) # write TLorentzVector tree.Branch("photon2",photon2) # define the calo energy resolution function def sEnergy(arEnergy, ePar): return math.sqrt(ePar[0]*math.sqrt(arEnergy[0])*ePar[0]*math.sqrt(arEnergy[0]) + ePar[1]*ePar[1] + ePar[2]*arEnergy[0]*ePar[2]*arEnergy[0]) # define the calo angular resolution function def sTheta(arEn, tPar): return ( tPar[0] / math.sqrt(arEn[0]) ) # create a function, need to set the number for the resolution in theta # and energy sigmaEnergy = ROOT.TF1('sigmaEnergy', sEnergy, 0. , 10. , 3 ) sigmaEnergy.SetParameters(EnergyPar) sigmaTheta = ROOT.TF1('sigmaTheta', sTheta, 0. , 3.2 , 1 ) sigmaTheta.SetParameter(0,ThetaPar) # Event loop for i in range(nEvent) : # generate decay angle in the center of mass system, decay should # be isotropic thetaCM[0] = R.Uniform(0.,Pi) # Photon energies in the Lab system are E1[0] = gamma * pi0Mass * (1. + (beta * math.cos(thetaCM[0]))) / 2. E2[0] = gamma * pi0Mass * (1. - (beta * math.cos(thetaCM[0]))) / 2. # Scattering angle of the first Photon Theta1[0] = math.atan(math.sin(thetaCM[0]) / (gamma * ( beta + math.cos(thetaCM[0]) ))) Theta2[0] = math.atan(math.sin(thetaCM[0] + Pi ) / (gamma * ( beta + math.cos(thetaCM[0] + Pi ) ))) # Opening angle cosalpha = (1.0 - (pi0Mass2 /(2.*E1[0]*E2[0]))) alpha[0] = math.acos(cosalpha); # pi0mass calc mpi0Calc[0] = math.sqrt(2.* E1[0] * E2[0] * (1. - cosalpha)) # Get momenta of photons Px1 = E1[0] * math.cos(Theta1[0]) Py1 = 0.; Pz1 = E1[0] * math.sin(Theta1[0]) Px2 = E2[0] * math.cos(Theta2[0]) Py2 = 0.; Pz2 = E2[0] * math.sin(Theta2[0]) # Set FourVector of the photons and get the mass using the ROOT function # of TLorentzVector photon1.SetPxPyPzE(Px1,Py1,Pz1,E1[0]) photon2.SetPxPyPzE(Px2,Py2,Pz2,E2[0]) mpi0FourVec = (photon1+photon2).M() # # Transform to measured energies and angle # # Energy and Theta resolution sigmaE1 = sigmaEnergy.Eval(E1[0]) sigmaE2 = sigmaEnergy.Eval(E2[0]) E1Measure[0] = R.Gaus( E1[0] , sigmaE1) E2Measure[0] = R.Gaus( E2[0] , sigmaE2) # Theta resolution scales with 1/sqrt(E) sigmaTheta1 = sigmaTheta.Eval(E1[0]) sigmaTheta2 = sigmaTheta.Eval(E2[0]) Theta1Measure[0] = R.Gaus( Theta1[0] , sigmaTheta1) Theta2Measure[0] = R.Gaus( Theta2[0] , sigmaTheta2) alphaMeasure[0] = Theta2Measure[0] - Theta1Measure[0] # Get momenta of photons including energy resolution Px1Measure = E1Measure[0] * math.cos(Theta1Measure[0]) Py1Measure = 0.; Pz1Measure = E1Measure[0] * math.sin(Theta1Measure[0]) Px2Measure = E2Measure[0] * math.cos(Theta2Measure[0]) Py2Measure = 0.; Pz2Measure = E2Measure[0] * math.sin(Theta2Measure[0]) # Set FourVector of the photons, invariant mass photon1Measure.SetPxPyPzE(Px1Measure,Py1Measure,Pz1Measure,E1Measure[0]) photon2Measure.SetPxPyPzE(Px2Measure,Py2Measure,Pz2Measure,E2Measure[0]) mpi0CalcMeasure[0] = (photon1Measure+photon2Measure).M() # Monitor print out if i == 1 : print("Generated theta in CM: ",thetaCM[0]) print("Photon1: E = ",E1[0]," theta = ", Theta1[0], "Photon2 E = ", E2[0], " theta = ",Theta2[0]) print ("measured 4 Vector Photon 1: ", E1Measure[0] , " " , Px1Measure, " " , Py1Measure, " " , Pz1Measure , "Photon 2: " ,E2Measure[0] , " " , Px2Measure , " " , Py2Measure , " " , Pz2Measure ) print("Calculated Mass pi0 = ",mpi0Calc[0]," Input pi0 Mass ",pi0Mass," using FourVectors ", mpi0FourVec ," measured mpi0 " , mpi0CalcMeasure[0] ) # Fill data tree tree.Fill() # Write Tree f.Write()