SND@LHC Software
Loading...
Searching...
No Matches
Functions | Variables
muonDis Namespace Reference

Functions

 photoabsorb (eps)
 
 G (x)
 
 nucl_cross (Ebeam, eps, A)
 
 SigmaAnalytic (Ebeam=5000., A=1, nsteps=10000)
 
 SigmaAnalyticVsEnergy (A=1)
 
 SigmaAnalyticVsA (Ebeam=500)
 
 myPrint (tc, tname, pathToPlots="/mnt/hgfs/microDisk/CERNBOX/SND@LHC/MuonDis/")
 
 convertAscii2Root (fname, version=2)
 
 muonPreTransport ()
 
 getMasssq (pid)
 
 rotate (ctheta, stheta, cphi, sphi, px, py, pz)
 
 getPythia6CrossSec (nstat, pmom=[])
 
 getPythia8CrossSec (nstat, pmom=[])
 
 readXsec (p)
 
 makeMuDISEvents (withElossFunction=False, nucleon='p+')
 
 checkProdofMuDIS ()
 
 compMuDIS_P6withP8 ()
 
 analyze (inFile)
 
 muonRateAtSND (withFaser=False, withEff=False, version=1)
 
 boundaries ()
 
 flukaMuons (version=1, Plimit=False, withFaser=True)
 
 muInterGeant4 (version=2, njobs=100)
 
 muondEdX (version=2, njobs=100, path='', withFaser=False, plotOnly=True)
 
 drawMuon3D (fname='muondEdX.root', hname='3d', gDir='muMinusGeant4_0')
 
 plotMuDisCrossSection ()
 
 muonDISfull (cycle=0, sMin=0, sMax=200, rMin=1, rMax=11, path='/eos/experiment/ship/user/truf/SND/muonDis/', debug=0, version=1, pythia6=True)
 
 thermNeutron ()
 
 analyzeDIS (NsubJobs=0, delta=13, hists="../Muons Extended Scoring Plane/muonDISfull.root", runCoverage=6.)
 
 count_python_processes (macroName)
 
 muonDISProduction (cycle, ecut=1., strippedEvents=False)
 
 energyOfDISevents ()
 
 selectEvents (Ecut=10)
 
 signalNeutrinos ()
 
 missing3Dproj (hist, ymin, ymax)
 

Variables

int theSeed = 0
 
 PDG = ROOT.TDatabasePDG.Instance()
 
 rnr = ROOT.TRandom()
 
dict flukaDict = {10:-13,11:13}
 
list pids = [-3222, -2212, -321, -211,-13, -11, 11,13, 211, 321, 2212, 3112]
 
dict pidsDict = {}
 
dict pidsDictRev = {}
 
float muonMass = 0.105658
 
dict normalization = {}
 
dict h = {}
 
float SND_Z = 483.262*u.m
 
int Eloss = 30.
 
str function = "13.6/1000.*sqrt([0])/x*(1.+0.038*log([0]))"
 
int xOverX0 = 75*u.m/10.02*u.cm
 
list EnergyScan = [10,20,30,40,50,75,100,200,300,400,500,1000,2000,3000,4000,5000,6000,7000,8000,10000]
 
 parser = ArgumentParser()
 
 dest
 
 type
 
 int
 
 help
 
 default
 
 float
 
 str
 
 options = parser.parse_args()
 
dict masssq = {}
 
dict FASERNU = {}
 
tuple scaleFactor = (42.*42.)/(24.*24.)
 
 ncpus = multiprocessing.cpu_count()
 
 cycle
 
 nMult
 
 sMin
 
 nStart
 
 sMax
 
 nEvents
 
 rMin
 
 rMax
 
 path
 
 muonIn
 
 pythia6
 
 nucleon
 

Function Documentation

◆ analyze()

muonDis.analyze (   inFile)

Definition at line 652 of file muonDis.py.

652def analyze(inFile):
653 NinteractionLength = 3
654 fin = ROOT.TFile('muDIScrossSec.root')
655 ROOT.gROOT.cd()
656 for pid in ['13','-13']:
657 h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid)
658 ut.bookHist(h,'inMu_'+str(pid),'Energy',200,0.,5000.,20,0.0,19.5)
659 ut.bookHist(h,'xy_In'+str(pid), ' x/y muon In',120,-60.,60.,120,-60.,60.)
660 ut.bookCanvas(h,'sec','xsec',900,600,1,1)
661 ut.bookHist(h,'muDISXsec',';E [GeV];#sigma [mb]',100,0.,10000.)
662 h['sec'].cd(1)
663 h['muDISXsec'].SetMinimum(0.)
664 h['muDISXsec'].SetMaximum(30.E-3)
665 h['muDISXsec'].Draw()
666 h['muDISXsec'].SetStats(0)
667 h['g_13'] .SetLineColor(ROOT.kGreen)
668 h['g_13'] .SetLineWidth(4)
669 h['g_13'] .Draw('same')
670 h['g_-13'] .Draw('same')
671 h['sec'].Print('muonXsec.png')
672
673 h['sTree'] = ROOT.TChain('DIS')
674 sTree = h['sTree']
675 for cycle in range(options.nMult):
676 for run in range(1,11):
677 for k in [0,1000]:
678 f = ROOT.TFile('muonDis_'+str(run+cycle*100+k)+'.root')
679 if not f.Get('DIS'):
680 print("file corrupted ",run+cycle*100+k)
681 continue
682 sTree.AddFile('muonDis_'+str(run+cycle*100+k)+'.root')
683 newFile = ROOT.TFile('muonDIS_SND.root','RECREATE')
684 newTree = sTree.CloneTree(0)
685 ROOT.gROOT.cd()
686 for nt in sTree:
687 muon = nt.InMuon[0]
688 muonEnergy = muon.Energy()
689 pid = muon.GetPdgCode()
690 wLHC = 5.83388*muon.GetWeight()/options.nMult/10. # ( options.nMult = number of cycles, 10 events per incoming muon in each cycle)
691 wDis = NinteractionLength * 97.5 / 1.67E-24 * h['g_'+str(pid)].Eval(muonEnergy ) * 1E-27
692 out = nt.Particles
693 rc = h['inMu_'+str(pid)].Fill(muonEnergy,out.GetEntries(),wLHC*wDis)
694 # place interaction point 5m in front of SND
695 z_int = -500.
696 lam = ( (SND_Z-z_int)-muon.Vz() ) / muon.Pz()
697 mxex = muon.Vx()+lam*muon.Px()
698 myex = muon.Vy()+lam*muon.Py()
699# missing update of time of flight
700 rc = h['xy_In'+str(muon.GetPdgCode())].Fill(mxex,myex,muon.GetWeight())
701 neutralParticles = []
702 chargedParticles = []
703 veto = False
704 for p in out:
705 pid = p.GetPdgCode()
706 E = p.Energy()
707 hname = "E_"+str(pid)
708 if not hname in h:
709 ut.bookHist(h,hname,'Energy',1000,0.,5000.)
710 ut.bookHist(h,'2dEsnd_'+str(pid), 'E in SND',1000,0.,5000.,100,0.,5000.)
711 ut.bookHist(h,'test1_'+str(pid), 'E in SND',100,0.,5000.)
712 ut.bookHist(h,'test2_'+str(pid), 'E in SND',100,0.,5000.)
713 ut.bookHist(h,'origin_xy'+str(pid), 'E in SND',100,-100.,100.,100,-100.,100.)
714 ut.bookHist(h,'2dEsnd_Veto_'+str(pid), 'E in SND',1000,0.,5000.,100,0.,5000.)
715 ut.bookHist(h,'Veto_mult_'+str(pid), 'veto multiplicity',100,0.,100.)
716 ut.bookHist(h,"xy_"+str(pid), 'x/y ',120,-60.,60.,120,-60.,60.)
717 rc = h[hname].Fill(E,p.GetWeight())
718 lamP = (SND_Z-z_int)/p.Pz()
719 xex = mxex + lamP*p.Px()
720 yex = myex + lamP*p.Py()
721 rc = h["xy_"+str(pid)].Fill(xex,yex,p.GetWeight())
722 if abs(pid) in [211,321,13,2212]: chargedParticles.append([pid,E])
723# counting in 15<y<57, 8<x<49
724# weight including interaction length, concrete interaction length = 42.4cm, 97.5 g/cm^2
725# Mproton = 1.67 10^-24 g, 1mb = 10^-27 cm^2
726# muon weight: add up the statistical weights
727# divide by the number of simulated p-p collisions (137,130,000) and
728# multiply by the actual collision rate (i.e. 8E8 at the nominal lumi of 1E34 cm-2 s-1): 5.83388
729 if -8 > xex and xex >-49 and 15 < yex and yex < 57:
730 if pid in [130,2112,-2112]: neutralParticles.append([pid,E,muonEnergy,wDis*wLHC])
731 rc = h["2dEsnd_"+str(pid)].Fill(E,muonEnergy ,wDis*wLHC)
732 rc = h["test1_"+str(pid)].Fill(muonEnergy ,wLHC)
733 rc = h["test2_"+str(pid)].Fill(muonEnergy)
734 if abs(pid) in [211,321,13,2212]:
735 if E > 1*u.GeV: veto = True
736# end loop over out particles
737 for x in neutralParticles:
738 rc = h["Veto_mult_"+str(x[0])].Fill(len(chargedParticles))
739 if len(chargedParticles)>0: rc = h['origin_xy'+str(x[0])].Fill(mxex,myex,x[3])
740 if not veto and len(neutralParticles)==1:
741 rc = h["2dEsnd_Veto_"+str(x[0])].Fill(x[1],x[2],x[3])
742 if len(neutralParticles)>0: rc = newTree.Fill()
743
744#
745 newFile.Write()
746 ut.bookCanvas(h,'muDIS_SND2','incoming muon energy',1500,900,2,1)
747 c=h['muDIS_SND2'].cd(1)
748 hname = "inMuEsnd_2112" # example
749 h[hname] = h["2dEsnd_2112" ].ProjectionY(hname)
750 ut.makeIntegralDistrib(h,hname)
751 h['I-'+hname].SetTitle('incoming muon energy;> E [GeV/c];N arbitrary units')
752 h['I-'+hname].SetMinimum(1.E-6)
753 # h['I-'+hname].SetMaximum(90)
754 h['I-'+hname].SetStats(0)
755 h['I-'+hname].SetLineWidth(3)
756 h['I-'+hname].Draw('hist')
757 c=h['muDIS_SND2'].cd(2)
758 h[hname].SetTitle('incoming muon energy; E [GeV/c];N arbitrary units')
759 h[hname].SetStats(0)
760 h[hname].SetLineWidth(3)
761 h[hname].Draw('hist')
762 h['muDIS_SND2'].Print('inMuEnergy.png')
763 parts = [130,2112,-2112]
764 for pid in parts:
765 hname = "Esnd_"+str(pid)
766 h[hname] = h["2dEsnd_"+str(pid)].ProjectionX(hname)
767 h["Esnd_Veto_"+str(pid)] = h["2dEsnd_Veto_"+str(pid)].ProjectionX("Esnd_Veto_"+str(pid))
768 ut.makeIntegralDistrib(h,hname)
769 h['I-'+hname].GetXaxis().SetRangeUser(0.,1000.)
770 h['I-'+hname].SetMinimum(1E-6)
771 ut.makeIntegralDistrib(h,"Esnd_Veto_"+str(pid))
772 ut.bookCanvas(h,'muDIS_SND','Kl and neutrons arriving at SND',2500,900,3,1)
773 for k in range(len(parts)):
774# 1E34 cm-2 s-1, 1fb = 1e-39 cm2, means 1fb requires 1E5 sec
775 fbScale = 1E5
776 tc=h['muDIS_SND'].cd(k+1)
777 tc.SetLogy(1)
778 pname = PDG.GetParticle(parts[k]).GetName()
779 hname = "Esnd_"+str(parts[k])
780 hnameVeto = "Esnd_Veto_"+str(parts[k])
781 for X in [hname,hnameVeto]:
782 h['IFB-'+X] = h['I-'+X].Clone('IFB-'+X)
783 h['IFB-'+X].Scale(150.*fbScale)
784 h['IFB-'+X].SetTitle(pname+';> E [GeV/c];N /150 fb^{-1}')
785 h['IFB-'+X]. GetYaxis().SetTitleOffset(1.4)
786 h['IFB-'+X].SetMaximum(90*150.)
787 h['IFB-'+X].SetStats(0)
788 h['IFB-'+X].SetLineWidth(3)
789 h['IFB-'+X].SetMinimum(0.1)
790 h['IFB-'+hname].Draw('hist')
791 h['IFB-'+hnameVeto].SetLineColor(ROOT.kRed)
792 h['IFB-'+hnameVeto].Draw('histsame')
793 hist = h['I-'+hname]
794 histVeto = h['I-Esnd_Veto_'+str(parts[k])]
795 R = hist.GetBinContent(1)
796 RVeto = histVeto.GetBinContent(1)
797 print(" %s: %5.2F N/sec %5.2F N fb with Veto: %5.2F N fb "%(pname,R,R*fbScale,RVeto*fbScale))
798 R = hist.GetBinContent(3)
799 RVeto = histVeto.GetBinContent(3)
800 print("E>10GeV, %s: %5.2F N/sec %5.2F N fb with Veto: %5.2F N fb "%(pname,R,R*fbScale,RVeto*fbScale))
801 R = hist.GetBinContent(21)
802 RVeto= histVeto.GetBinContent(21)
803 print("E>100GeV, %s: %5.2E N/sec %5.2F N fb with Veto: %5.2F N fb "%(pname,R,R*fbScale,RVeto*fbScale))
804 R = hist.GetBinContent(41)
805 RVeto= histVeto.GetBinContent(41)
806 print("E>200GeV, %s: %5.2E N/sec %5.2F N fb with Veto: %5.2F N fb "%(pname,R,R*fbScale,RVeto*fbScale))
807 R = hist.GetBinContent(101)
808 RVeto = histVeto.GetBinContent(101)
809 print("E>500GeV, %s: %5.2E N/sec %5.2F N fb with Veto: %5.2F N fb "%(pname,R,R*fbScale,RVeto*fbScale))
810 h['muDIS_SND'].Print('muDIS_SND.png')
811 h['muDIS_SND'].Print('muDIS_SND.pdf')

◆ analyzeDIS()

muonDis.analyzeDIS (   NsubJobs = 0,
  delta = 13,
  hists = "../Muons Extended Scoring Plane/muonDISfull.root",
  runCoverage = 6. 
)

Definition at line 1814 of file muonDis.py.

1814def analyzeDIS(NsubJobs=0,delta=13,hists="../Muons Extended Scoring Plane/muonDISfull.root",runCoverage=6.):
1815# pathToPlots = "/mnt/hgfs/microDisk/CERNBOX/SND@LHC/MuonDis/"
1816 pathToPlots = "/mnt/hgfs/microDisk/SND@LHC/MuonDis/reMake2022"
1817 latex = 'latex.txt'
1818 if options.pythia6<0:
1819 pathToPlots = "/mnt/hgfs/microDisk/CERNBOX/SND@LHC/MuonGeant4/"
1820 latex = 'latex-geant4.tex'
1821 if NsubJobs==0:
1822 # one cycle = 10x muon statistics
1823 ut.readHists(h,hists)
1824 else:
1825 s = 0
1826 while s < NsubJobs:
1827 sMin,sMax = s,s+delta
1828 print("reading muonDISfull-"+str(sMin)+"_"+str(sMax)+".root" )
1829 ut.readHists(h,"muonDISfull-"+str(sMin)+"_"+str(sMax)+".root")
1830 s+=delta
1831
1832 ut.bookCanvas(h,'muDIS_SND2','incoming muon energy',1500,900,2,1)
1833 c=h['muDIS_SND2'].cd(1)
1834 hname = "inMu"
1835 h[hname] = h["inMu_013" ].ProjectionX(hname)
1836 h[hname].Add(h["inMu_0-13" ].ProjectionX())
1837 ut.makeIntegralDistrib(h,hname)
1838 h['I-'+hname].SetTitle('incoming muon energy;> E [GeV/c];N arbitrary units')
1839 h['I-'+hname].SetMinimum(1.E-6)
1840 h['I-'+hname].SetStats(0)
1841 h['I-'+hname].SetLineWidth(3)
1842 h['I-'+hname].Draw('hist')
1843 c=h['muDIS_SND2'].cd(2)
1844 h[hname].SetTitle('incoming muon energy; E [GeV/c];N arbitrary units')
1845 h[hname].SetStats(0)
1846 h[hname].SetLineWidth(3)
1847 h[hname].Draw('hist')
1848 myPrint(h['muDIS_SND2'],'inMuEnergy.png',pathToPlots=pathToPlots)
1849#
1850 c=h['muDIS_SND2'].cd(1)
1851 hname = "xyz_mu_0"
1852 for x in ['13','-13']:
1853 h["xy_mu_0"+x] = h["xyz_mu_0"+x].Project3D('yx')
1854 h["xy_mu_0"+x].SetName("xy_mu_0"+x)
1855 h["xy_mu_0"+x].SetStats(0)
1856 h["xy_mu_0"+x].SetTitle(PDG.GetParticle(int(x)).GetName()+' ;x [cm]; y [cm]')
1857 h["xy_mu_0"+x].Draw('colz')
1858 c=h['muDIS_SND2'].cd(2)
1859 myPrint(h['muDIS_SND2'],'inMu_0XY',pathToPlots=pathToPlots)
1860 h["xyz_mu_0"]=h["xyz_mu_013"].Clone("xyz_mu_0")
1861 h["xyz_mu_0"].Add(h["xyz_mu_0-13"])
1862 h["yz_mu_0"] = h["xyz_mu_0"].Project3D('yz')
1863 h["yz_mu_0"].SetName("yz_mu_0")
1864 h["yz_mu_0"].SetStats(0)
1865 h["yz_mu_0"].SetTitle("; z [cm]; y [cm]")
1866 ut.bookCanvas(h,'muDIS_SNDyz','incoming muon',1200,900,1,1)
1867 c=h['muDIS_SNDyz'].cd()
1868 h["yz_mu_0"].Draw('colz')
1869 myPrint(h['muDIS_SNDyz'],'inMu_0YZ',pathToPlots=pathToPlots)
1870
1871 parts = [130,2112,-2112,310,3122,-3122,3322,-3322,22]
1872# 1E34 cm-2 s-1, 1fb = 1e-39 cm2, means 1fb requires 1E5 sec
1873 fbScale = 1E5
1874 ut.bookCanvas(h,'muDIS_N0In','primary neutrals',1200,1200,3,3)
1875 for k in range(len(parts)):
1876 pid = parts[k]
1877 pname = PDG.GetParticle(pid).GetName()
1878 tc=h['muDIS_N0In'].cd(k+1)
1879 tc.SetLogy(1)
1880 hname = "inE_"+str(pid) # all primary particles
1881 hConcr = "inE_Concrete_"+str(pid) # origin from concrete
1882 hout = "inE_out_"+str(pid) # going out
1883 h[hout] = h[hname].Clone(hout)
1884 h[hout].Add(h[hConcr],-1.)
1885 h[hout+'_projx'] = h[hout].ProjectionX(hout+'_projx')
1886 h[hname+'_projx'] = h[hname].ProjectionX(hname+'_projx')
1887 h[hout+'_projx'].SetMarkerStyle(20)
1888 h[hname+'_projx'].SetMarkerStyle(27)
1889 h[hname+'_projx'].SetStats(0)
1890 binw = int(h[hname+'_projx'].GetBinWidth(1))
1891 h[hname+'_projx'].GetYaxis().SetTitle('N [Hz/'+str(binw)+'GeV]')
1892 h[hout+'_projx'].SetMarkerStyle(29)
1893 h[hname+'_projx'].GetXaxis().SetRangeUser(0.,2900.)
1894 h[hname+'_projx'].GetXaxis().SetTitle('Energy [GeV] ')
1895 h[hname+'_projx'].SetTitle(pname)
1896 h[hname+'_projx'].GetYaxis().SetTitleOffset(1.2)
1897 h[hname+'_projx'].SetMaximum(20.)
1898 h[hname+'_projx'].Draw('PHIST')
1899 h[hname+'_projx'].Draw('histsame')
1900 h[hout+'_projx'].SetLineColor(ROOT.kRed)
1901 h[hout+'_projx'].Draw('SAMEPHIST')
1902 h[hout+'_projx'].Draw('histsame')
1903 myPrint(h['muDIS_N0In'],'EofPrimNeutrals',pathToPlots=pathToPlots)
1904
1905 tname = 'muDIS_inSND'
1906 ut.bookCanvas(h,tname,'neutrals not in concrete and z>-25cm',1200,1200,3,3)
1907 for k in range(len(parts)):
1908 pid = parts[k]
1909 pname = PDG.GetParticle(pid).GetName()
1910 tc=h[tname].cd(k+1)
1911 tc.SetLogy(1)
1912 h["inE_prim_"+str(pid)] = h["inE_"+str(pid)].Clone("inE_prim_"+str(pid))
1913 h["inE_prim_"+str(pid)].Add(h["inE_Concrete_"+str(pid)],-1)
1914 h["inE_seco_"+str(pid)] = h["inE_"+str(pid)+'_secondary'].Clone("inE_"+str(pid)+'_secondary_out')
1915 h["inE_seco_"+str(pid)].Add(h["inE_Concrete_"+str(pid)+'_secondary'],-1)
1916 for x in ['prim_','seco_']:
1917 zmin = h["inE_"+x+str(pid)].GetYaxis().FindBin(-25.)
1918 zmax = h["inE_"+x+str(pid)].GetYaxis().FindBin(+25.)
1919 hname = "inE_"+x+str(pid)+'_SND'
1920 h[hname] = h["inE_"+x+str(pid)].ProjectionX(hname,zmin,zmax)
1921 h[hname].SetStats(0)
1922 binw = int(h[hname].GetBinWidth(1))
1923 h[hname].GetYaxis().SetTitle('N [Hz/'+str(binw)+'GeV]')
1924 h[hname].GetXaxis().SetRangeUser(0.,2900.)
1925 h[hname].GetXaxis().SetTitle('Energy [GeV] ')
1926 h[hname].SetTitle(pname)
1927 h[hname].GetYaxis().SetTitleOffset(1.2)
1928 h[hname].SetMaximum(20.)
1929 h[hname].SetMinimum(1.E-6)
1930 hname = "inE_prim_"+str(pid)+'_SND'
1931 h[hname].SetLineColor(ROOT.kRed)
1932 h[hname].Draw('hist')
1933 hname = "inE_seco_"+str(pid)+'_SND'
1934 h[hname].SetLineColor(ROOT.kBlue)
1935 h[hname].Draw('histsame')
1936 myPrint(h[tname],'EofNeutralsInSND',pathToPlots=pathToPlots)
1937
1938 colour = {}
1939 colour[''] = ROOT.kRed
1940 colour['_secondary'] = ROOT.kBlue
1941 minMax = {'out_':[9,1E-6],'Concrete_':[30.,1E-5]}
1942 for k in range(len(parts)):
1943 pid = parts[k]
1944 for origin in ['_secondary','']:
1945 hname = "inE_"+str(pid)+origin # all primary particles
1946 hConcr = "inE_Concrete_"+str(pid) +origin # end in concrete
1947 hout = "inE_out_"+str(pid)+origin # end somewhere else
1948 h[hout] = h[hname].Clone(hout)
1949 h[hout].Add(h[hConcr],-1.)
1950
1951 for c in ['out_','Concrete_']:
1952 tname = 'muDIS_neuZend'+c
1953 ut.bookCanvas(h,tname,'end vertex z',1200,1200,3,3)
1954 for k in range(len(parts)):
1955 pid = parts[k]
1956 pname = PDG.GetParticle(pid).GetName()
1957 tc=h[tname].cd(k+1)
1958 tc.SetLogy(1)
1959 for origin in ['_secondary','']:
1960# only primaries
1961 hname = "inE_"+c+str(pid)+origin
1962 h[hname+'_projz'] = h[hname].ProjectionY(hname+'_projz')
1963 h[hname+'_projz'].SetLineColor(colour[origin])
1964 h[hname+'_projz'].SetStats(0)
1965 binw = int(h[hname+'_projz'].GetBinWidth(1))
1966 h[hname+'_projz'].GetYaxis().SetTitle('N [Hz/'+str(binw)+'GeV]')
1967 h[hname+'_projz'].GetXaxis().SetTitle('z [cm] ')
1968 h[hname+'_projz'].SetTitle(pname)
1969 h[hname+'_projz'].GetYaxis().SetTitleOffset(1.2)
1970 h[hname+'_projz'].SetMaximum(minMax[c][0])
1971 h[hname+'_projz'].SetMinimum(minMax[c][1])
1972 if origin=='_secondary': h[hname+'_projz'].Draw('HIST')
1973 else: h[hname+'_projz'].Draw('HISTsame')
1974 myPrint(h[tname],c+'zEndOfNeutrals',pathToPlots=pathToPlots)
1975
1976
1977# vertices, note: 'xyz_muInter','xyz_Inter','xyz_origin' are for events with hits
1978 for hist in ['xyz_muInter','xyz_Inter','xyzE100_Inter','xyz_origin','inE','inE_Concrete','xyzVeto_Inter','z_veto']:
1979 first = True
1980 for pid in parts:
1981 if pid==22: continue
1982 if first:
1983 h[hist] = h[hist+'_'+str(pid)].Clone(hist)
1984 first = False
1985 else: h[hist].Add(h[hist+'_'+str(pid)])
1986 h['inE-noConc']=h['inE'].Clone('inE-noConc')
1987 h['inE-noConc'].Add(h['inE_Concrete'],-1.)
1988 h['inE-noConc_22']=h['inE_22'].Clone('inE-noConc_22')
1989 h['inE-noConc_22'].Add(h['inE_Concrete_22'],-1.)
1990 for hist in ['xyz_muInter','xyz_Inter','xyzE100_Inter','xyz_origin','inE','inE-noConc','xyzVeto_Inter','z_veto']:
1991 if h[hist].ClassName() == 'TH2D':
1992 h[hist+'_z']=h[hist].ProjectionY(hist+'_z')
1993 h[hist+'_22_z']=h[hist+'_22'].ProjectionY(hist+'_22_z')
1994 h[hist+'_22_z'].SetTitle(';z [cm]')
1995 h[hist+'_z'].SetTitle(';z [cm]')
1996 else:
1997 for g in ['','_22']:
1998 h[hist+g+'_z'] =h[hist+g].ProjectionZ(hist+g+'_z')
1999 h[hist+g+'_z'].SetTitle(';z [cm]')
2000 h[hist+g+'_xy']=h[hist+g].Project3D('yx')
2001 h[hist+g+'_xy'].SetName(hist+g+'_xy')
2002 h[hist+g+'_xy'].SetTitle(';x [cm]; y [cm]')
2003 h[hist+g+'_xy'].SetStats(0)
2004 h[hist+g+'_xz']=h[hist+g].Project3D('xz')
2005 h[hist+g+'_xz'].SetName(hist+g+'_xz')
2006 h[hist+g+'_xz'].SetTitle(';z [cm]; x [cm]')
2007 h[hist+g+'_xz'].SetStats(0)
2008 h[hist+g+'_yz']=h[hist+g].Project3D('yz')
2009 h[hist+g+'_yz'].SetName(hist+g+'_yz')
2010 h[hist+g+'_yz'].SetTitle(';z [cm]; y [cm]')
2011 h[hist+g+'_yz'].SetStats(0)
2012 h[hist+'_z'].SetStats(0)
2013 h[hist+'_z'].SetTitle('; z [cm]')
2014 ut.bookCanvas(h,'vertices'+hist+'_z','vertices '+hist,1200,900,1,1)
2015 tc = h['vertices'+hist+'_z'].cd()
2016 if hist=='z_veto': tc.SetLogy()
2017 if hist=='xyzVeto_Inter': tc.SetLogy()
2018 if hist in ['xyz_muInter','xyz_origin']:
2019 h[hist+'_z'].Draw()
2020 else:
2021 h[hist+'_22_z'].SetStats(0)
2022 h[hist+'_22_z'].SetLineColor(ROOT.kRed)
2023 h[hist+'_22_z'].Draw()
2024 h[hist+'_z'].Draw('same')
2025 if hist=='xyzE100_Inter':
2026 txtlam = 'empty'
2027 txtrad = 'empty'
2028 if h[hist+'_z'].GetEntries()>10:
2029 rc = h[hist+'_z'].Fit('expo','S','',-25.,30.)
2030 fitresult = rc.Get()
2031 txtlam = '#lambda =%4.1Fcm'%(-1./fitresult.GetParams()[1])
2032 rc = h[hist+'_22_z'].Fit('expo','SL','',-30.,10.)
2033 fitresult = rc.Get()
2034 funRad =h[hist+'_22_z'].GetFunction('expo')
2035 funRad.SetBit(ROOT.TF1.kNotDraw)
2036 funRad.SetLineColor(h[hist+'_22_z'].GetLineColor())
2037 txtrad = 'X_{0} =%4.1Fcm'%(-1./fitresult.GetParams()[1])
2038 h[hist+'_z'].GetFunction('expo').SetLineColor(h[hist+'_z'].GetLineColor())
2039 h[hist+'_z'].Draw()
2040 h[hist+'_22_z'].Draw('same')
2041 T = ROOT.TLatex()
2042 T.DrawLatex(-10,0.005,txtlam)
2043 #T.DrawLatex(-10,0.0045,txtrad)
2044 myPrint(h['vertices'+hist+'_z'],hist+'_z',pathToPlots=pathToPlots)
2045 if h[hist].ClassName() == 'TH3D':
2046 for proj in ['xy','xz','yz']:
2047 ut.bookCanvas(h,'vertices'+hist+'_'+proj,'vertices '+hist,1200,900,1,1)
2048 tc = h['vertices'+hist+'_'+proj].cd()
2049 h[hist+'_'+proj].Draw('colz')
2050 myPrint(h['vertices'+hist+'_'+proj],hist+'_'+proj,pathToPlots=pathToPlots)
2051
2052 ut.bookCanvas(h,'verticesinE_z','vertices inE_z',1200,900,1,1)
2053 h['verticesinE_z'].cd()
2054 h["inE-noConc_z"].SetLineColor(ROOT.kRed)
2055 h["inE_z"].Draw('hist')
2056 h["inE-noConc_z"].Draw('samehist')
2057 myPrint(h['verticesinE_z'],"inE_z",pathToPlots=pathToPlots)
2058 ut.bookCanvas(h,'verticesxyz_Inter_z','vertices xyz_Inter_z',1200,900,1,1)
2059 h['verticesxyz_Inter_z'].cd()
2060 h['xyz_Inter_z'].GetXaxis().SetRangeUser(-50.,100.)
2061 h["xyz_Inter_z"].Draw('hist')
2062 myPrint(h['verticesxyz_Inter_z'],"verticesxyz_Inter_z",pathToPlots=pathToPlots)
2063# multiplicities
2064 withMult=False
2065 if withMult:
2066 ut.bookCanvas(h,'mult',' ',1200,900,1,1)
2067 tc = h['mult'].cd()
2068 tc.SetLogy(1)
2069 mults = {22:[19,ROOT.kMagenta],130:[20,ROOT.kRed],2112:[22,ROOT.kBlue],-2112:[23,ROOT.kCyan],
2070 310:[21,ROOT.kOrange],3122:[24,ROOT.kGreen-1],-3122:[25,ROOT.kGreen+1],
2071 3322:[26,ROOT.kGreen-2],-3322:[27,ROOT.kGreen+2]}
2072 h['legmult']=ROOT.TLegend(0.6,0.6,0.82,0.75)
2073 for m in mults:
2074 hist = 'mult_'+str(m)
2075 h[hist].SetStats(0)
2076 h[hist].SetTitle('; N')
2077 h[hist].SetMarkerStyle(mults[m][0])
2078 h[hist].SetMarkerColor(mults[m][1])
2079 h[hist].SetLineColor(mults[m][1])
2080 rc = h['legmult'].AddEntry(h[hist],PDG.GetParticle(m).GetName(),'PL')
2081 h['mult_22'].Draw()
2082 for m in mults: h['mult_'+str(m)].Draw('same')
2083 h['legmult'].Draw('same')
2084 myPrint(h['mult'],"neutralMultiplicities",pathToPlots=pathToPlots)
2085 for pid in parts:
2086 ut.makeIntegralDistrib(h,"inE_"+str(pid))
2087 ut.makeIntegralDistrib(h,"inE_Concrete_"+str(pid))
2088 for o in ["","prim"]:
2089 hname = "Esnd"+o+"_"+str(pid)
2090 h[hname] = h["E"+o+"_"+str(pid)].ProjectionX(hname)
2091 h["Esnd"+o+"_all_"+str(pid)] = h["E"+o+"_veto_"+str(pid)].ProjectionX("Esnd"+o+"_all_"+str(pid))
2092 h["Esnd"+o+"_all_"+str(pid)] .Add(h[hname] )
2093 ut.makeIntegralDistrib(h,hname)
2094 h['I-'+hname].GetXaxis().SetRangeUser(0.,1000.)
2095 h['I-'+hname].SetMinimum(1E-6)
2096 ut.makeIntegralDistrib(h,"Esnd"+o+"_all_"+str(pid))
2097
2098 ut.bookCanvas(h,'muDIS_SNDwithVeto','neutrals arriving at SND',1800,1200,3,3)
2099 for k in range(len(parts)):
2100 tc=h['muDIS_SNDwithVeto'].cd(k+1)
2101 tc.SetLogy(1)
2102 pname = PDG.GetParticle(parts[k]).GetName()
2103 hnameScaled = "Esnd_"+str(parts[k])+"_fb150"
2104 h[hnameScaled] = h["Esnd_"+str(parts[k])].Clone(hnameScaled)
2105 h[hnameScaled].Scale(150.*fbScale/runCoverage)
2106 h[hnameScaled].GetXaxis().SetRangeUser(0,500)
2107 h[hnameScaled].SetStats(0)
2108 h[hnameScaled].SetMaximum(4E6)
2109 h[hnameScaled].SetTitle(pname+'; E [GeV/c];N/10GeV /150 fb^{-1}')
2110 h[hnameScaled].SetLineWidth(3)
2111 h[hnameScaled].Draw('hist')
2112 myPrint(h['muDIS_SNDwithVeto'],'EofNeutralsInSND_vetoApplied',pathToPlots=pathToPlots)
2113
2114 ut.bookCanvas(h,'muDIS_SNDnoVetoRequired','neutrals arriving at SND',1800,1200,3,3)
2115 for k in range(len(parts)):
2116 tc=h['muDIS_SNDnoVetoRequired'].cd(k+1)
2117 tc.SetLogy(1)
2118 pname = PDG.GetParticle(parts[k]).GetName()
2119 hnameScaled = "Esnd_all_"+str(parts[k])+"_fb150"
2120 h[hnameScaled] = h["Esnd_all_"+str(parts[k])].Clone(hnameScaled)
2121 h[hnameScaled].Scale(150.*fbScale/runCoverage)
2122 h[hnameScaled].GetXaxis().SetRangeUser(0,500)
2123 h[hnameScaled].SetStats(0)
2124 h[hnameScaled].SetMaximum(4E6)
2125 h[hnameScaled].SetTitle(pname+'; E [GeV/c];N/10GeV /150 fb^{-1}')
2126 h[hnameScaled].SetLineWidth(3)
2127 h[hnameScaled].Draw('hist')
2128 myPrint(h['muDIS_SNDnoVetoRequired'],'EofNeutralsInSND_NoVetoApplied',pathToPlots=pathToPlots)
2129
2130 for o in ["","prim"]:
2131 ut.bookCanvas(h,'muDIS_SND'+o,'neutrals arriving at SND',1800,1200,3,3)
2132 for k in range(len(parts)):
2133 tc=h['muDIS_SND'+o].cd(k+1)
2134 tc.SetLogy(1)
2135 tc.SetLogx(0)
2136 pname = PDG.GetParticle(parts[k]).GetName()
2137 hname = "Esnd"+o+"_"+str(parts[k])
2138 hnameAll = "Esnd"+o+"_all_"+str(parts[k])
2139 for X in [hname,hnameAll]:
2140 h['IFB-'+X] = h['I-'+X].Clone('IFB-'+X)
2141 h['IFB-'+X].Scale(150.*fbScale)
2142 h['IFB-'+X].SetTitle(pname+';> E [GeV/c];N /150 fb^{-1}')
2143 h['IFB-'+X]. GetYaxis().SetTitleOffset(1.4)
2144 h['IFB-'+X]. GetXaxis().SetRangeUser(0.1,1200.)
2145 h['IFB-'+X].SetMaximum(2E7)
2146 h['IFB-'+X].SetStats(0)
2147 h['IFB-'+X].SetLineWidth(3)
2148 h['IFB-'+X].SetMinimum(1.0)
2149 h['IFB-'+hnameAll].Draw('hist')
2150 h['IFB-'+hname].SetLineColor(ROOT.kRed)
2151 h['IFB-'+hname].Draw('histsame')
2152 hist = h['I-'+hname]
2153 histAll = h['I-'+hnameAll]
2154 print("----> results for "+o)
2155 for E in [0,10,100,200,300,500,1000]:
2156 n = hist.FindBin(E)
2157 RVeto = hist.GetBinContent(n) / runCoverage
2158 R = histAll.GetBinContent(n) / runCoverage
2159 faser = ''
2160 if parts[k] in FASERNU:
2161 if E in FASERNU[parts[k]]:
2162 faser = "%5.2G"%( FASERNU[parts[k]][E]*scaleFactor )
2163 print("E>%i GeV, %s: %5.2F N/sec %5.2F N fb with Veto: %5.2F N fb | %5.2G veto %5.2G (%s) "%(E,pname,R,R*fbScale,RVeto*fbScale,R*150.*fbScale,RVeto*150.*fbScale,faser))
2164 myPrint(h['muDIS_SND'+o],'muDIS_SND'+o,pathToPlots=pathToPlots)
2165# make latex output
2166 fout = open(latex,'w')
2167 platex = {2112:'$n$' ,-2112: '$\overline{n}$' ,130:'$K_L$',310:'$K_S$',3122:'$\Lambda$',-3122:'$\overline{\Lambda}$',3322:'$\Xi + \overline{\Xi}$',22:'$\gamma$'}
2168
2169 for x in [2112,-2112,130,310,3122,-3122,3322,22]:
2170 line = platex[x]
2171 for E in [10,100,200,300,500,1000]:
2172 pname = PDG.GetParticle(x).GetName()
2173 hnameAll = "Esnd"+"_all_"+str(x)
2174 histAll = h['I-'+hnameAll]
2175 n = histAll.FindBin(E)
2176 R = histAll.GetBinContent(n) / runCoverage
2177 if x == 3322: R+=h['I-'+hnameAll.replace('3322','-3322')].GetBinContent(n) / runCoverage
2178 R150 = R*150.*fbScale
2179 if R150<1E-8: R150=0
2180 if R150<1000: line += "& $%5.1F$"%(R150)
2181 else: line += "& $%6.2G$"%(R150)
2182 line+= "\\\\ \n"
2183 #print(line)
2184 xline = line.replace('E+0','\,10^')
2185 fout.write(xline)
2186 fout.write("$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$\n")
2187 for x in [2112,-2112,130,310,3122,-3122,3322,22]:
2188 line = platex[x]
2189 for E in [10,100,200,300,500,1000]:
2190 pname = PDG.GetParticle(x).GetName()
2191 hname = "Esnd_"+str(x)
2192 hist= h['I-'+hname]
2193 n = hist.FindBin(E)
2194 R = hist.GetBinContent(n) / runCoverage
2195 if x == 3322: R+=h['I-'+hname.replace('3322','-3322')].GetBinContent(n) / runCoverage
2196 R150 = R*150.*fbScale
2197 if R150<1E-8: R150=0
2198 if R150<1000: line += "& $%5.1F$"%(R150)
2199 else: line += "& $%6.2G$"%(R150)
2200 line+= "\\\\ \n"
2201 #print(line)
2202 xline = line.replace('E+0','\,10^')
2203 fout.write(xline)
2204 fout.write("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n")
2205 for x in [2112,-2112,130,310,3122,-3122,3322,22]:
2206 line = platex[x]
2207 for E in [10,100,300,1000]:
2208 faser = FASERNU[x][E]*scaleFactor
2209 if x == 3322: faser+=FASERNU[-x][E]*scaleFactor
2210 if faser<1000: line += "& $%5.1F$"%(faser)
2211 else: line += "& $%6.2G$"%(faser)
2212 line+= "\\\\ \n"
2213 #print(line)
2214 xline = line.replace('E+0','\,10^')
2215 fout.write(xline)
2216 fout.close()
2217

◆ boundaries()

muonDis.boundaries ( )

Definition at line 995 of file muonDis.py.

995def boundaries():
996 h['snd_1'] = ROOT.TLine(-49,15,-8,15)
997 h['snd_2'] = ROOT.TLine(-49,57,-8,57)
998 h['snd_3'] = ROOT.TLine(-49,57,-49,15)
999 h['snd_4'] = ROOT.TLine(-8,57,-8,15)
1000 h['fas_1'] = ROOT.TLine(-12,12,12,12)
1001 h['fas_2'] = ROOT.TLine(-12,-12,12,-12)
1002 h['fas_3'] = ROOT.TLine(-12,-12,-12,12)
1003 h['fas_4'] = ROOT.TLine(12,-12,12,12)
1004 for n in range(1,5):
1005 h['snd_'+str(n)] .SetLineColor(ROOT.kBlack)
1006 h['snd_'+str(n)] .SetLineWidth(5)
1007 h['snd_'+str(n)] .SetLineStyle(1)
1008 h['fas_'+str(n)] .SetLineColor(ROOT.kCyan)
1009 h['fas_'+str(n)] .SetLineWidth(5)
1010 h['fas_'+str(n)] .SetLineStyle(2)
1011

◆ checkProdofMuDIS()

muonDis.checkProdofMuDIS ( )

Definition at line 540 of file muonDis.py.

540def checkProdofMuDIS():
541 redfac3d=100
542 for pid in [13,-13]:
543 ut.bookHist(h,"xyz_mu_"+str(pid), 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
544 for r in range(1,11):
545 for z in [0,1000]:
546 f=ROOT.TFile('/home/truf/ubuntu-1710/ship-ubuntu-1710-64/SND/MuonDis/Muons Extended Scoring Plane/muonDis_'+str(z+r)+'.root')
547 for sTree in f.DIS:
548 for m in sTree.InMuon:
549 rc = h['xyz_mu_'+str(m.GetPdgCode())].Fill(m.Vx(),m.Vy(),m.Vz())
550 ut.bookCanvas(h,'muDIS_inMu','incoming muon',1500,900,2,1)
551 c=h['muDIS_inMu'].cd(1)
552 for x in ['13','-13']:
553 h["xy_mu_"+x] = h["xyz_mu_"+x].Project3D('yx')
554 h["xy_mu_"+x].SetName("xy_mu_"+x)
555 h["xy_mu_"+x].SetStats(0)
556 h["xy_mu_"+x].SetTitle(PDG.GetParticle(int(x)).GetName()+' ;x [cm]; y [cm]')
557 h["xy_mu_"+x].Draw('colz')
558 c=h['muDIS_inMu'].cd(2)
559 myPrint(h['muDIS_SND2'],'inMu_XY')
560

◆ compMuDIS_P6withP8()

muonDis.compMuDIS_P6withP8 ( )

Definition at line 561 of file muonDis.py.

561def compMuDIS_P6withP8():
562 f = {}
563 for x in ['6','8']:
564 h['p'+x] = {}
565 ut.readHists(h['p'+x],'muDIScrossSec_Pythia'+x+'.root')
566 for aHist in h['p'+x]:
567 for t in ["nMult","Qhist","pTehist","xhist","yhist","pTrhist","pTdhist"]:
568 if aHist.find(t)<0:continue
569 h['p'+x][aHist].Scale(1./h['p'+x][aHist].GetEntries())
570 f['p6'] = ROOT.TFile('muDIScrossSec_Pythia6.root')
571 f['p8'] = ROOT.TFile('muDIScrossSec_Pythia8.root')
572 ROOT.gROOT.cd()
573 for pid in ['13','-13']:
574 for t in ['p','n']:
575 for g in ['p6','p8']:
576 tx = t
577 if t=='p' and g=='p6': tx='p+'
578 h['g_'+g+pid+t] = f[g].Get('g_'+pid+tx).Clone('g_'+g+pid+t)
579 h['g_'+g+pid+t].SetLineWidth(4)
580 if g=='p6':
581 h['g_'+g+pid+t].SetLineColor(ROOT.kBlue)
582 h['g_'+g+pid+t].SetMarkerColor(ROOT.kBlue)
583 if g=='p8':
584 h['g_'+g+pid+t].SetLineColor(ROOT.kGreen)
585 h['g_'+g+pid+t].SetMarkerColor(ROOT.kGreen)
586 if pid=='13': h['g_'+g+pid+t].SetMarkerStyle(23)
587 if pid=='-13': h['g_'+g+pid+t].SetMarkerStyle(22)
588 ut.bookCanvas(h,'sec','xsec',900,600,1,1)
589 ut.bookHist(h,'muDISXsec',';E [GeV];#sigma [mb]',100,0.,10000.)
590 h['sec'].cd(1)
591 h['muDISXsec'].SetMinimum(5E-5)
592 h['muDISXsec'].SetMaximum(30.E-3)
593 h['muDISXsec'].SetStats(0)
594 h['muDISXsec'].Draw()
595 for pid in ['13','-13']:
596 for t in ['p','n']:
597 for g in ['p6','p8']:
598 h['g_'+g+pid+t].Draw('same')
599 T=ROOT.TLatex()
600 T.DrawLatex(2000,0.007,'Pythia6')
601 T.SetLineColor(ROOT.kBlue)
602 T.DrawLatex(4000,0.0003,'Pythia8')
603 T.SetLineColor(ROOT.kBlue)
604 h['sec'].Print('muonXsecP6P8.png')
605
606 E = str(500.0)
607 ut.bookCanvas(h,'mul',E,1200,600,3,1)
608 h['mul'].cd(1)
609 h['p6']['nMult13p+'+E].SetStats(0)
610 h['p6']['nMult13p+'+E].SetLineColor(ROOT.kBlue)
611 h['p6']['nMult13p+'+E].Draw()
612 h['p8']['nMult13p'+E].SetStats(0)
613 h['p8']['nMult13p'+E].SetLineColor(ROOT.kGreen)
614 h['p8']['nMult13p'+E].Draw('same')
615 T.SetLineColor(ROOT.kRed)
616 T.DrawLatexNDC(0.5,0.45,'mean P6: %5.2F'%(h['p6']['nMult13p+'+E].GetMean()))
617 T.DrawLatexNDC(0.5,0.5,'mean P8: %5.2F'%(h['p8']['nMult13p'+E].GetMean()))
618 h['mul'].cd(2)
619 h['p6']['nMultcha13p+'+E].SetStats(0)
620 h['p6']['nMultcha13p+'+E].SetLineColor(ROOT.kBlue)
621 h['p6']['nMultcha13p+'+E].GetXaxis().SetRangeUser(-0.5,30.5)
622 h['p6']['nMultcha13p+'+E].Draw()
623 h['p8']['nMultneu13p'+E].SetStats(0)
624 h['p8']['nMultcha13p'+E].SetLineColor(ROOT.kGreen)
625 h['p8']['nMultcha13p'+E].Draw('same')
626 T.DrawLatexNDC(0.5,0.45,'mean P6: %5.2F'%(h['p6']['nMultcha13p+'+E].GetMean()))
627 T.DrawLatexNDC(0.5,0.5,'mean P8: %5.2F'%(h['p8']['nMultcha13p'+E].GetMean()))
628 h['mul'].cd(3)
629 h['p6']['nMultneu13p+'+E].SetStats(0)
630 h['p6']['nMultneu13p+'+E].SetLineColor(ROOT.kBlue)
631 h['p6']['nMultneu13p+'+E].GetXaxis().SetRangeUser(-0.5,8.5)
632 h['p6']['nMultneu13p+'+E].Draw()
633 h['p8']['nMultneu13p'+E].SetStats(0)
634 h['p8']['nMultneu13p'+E].SetLineColor(ROOT.kGreen)
635 h['p8']['nMultneu13p'+E].Draw('same')
636 T.DrawLatexNDC(0.5,0.45,'mean P6: %5.2F'%(h['p6']['nMultneu13p+'+E].GetMean()))
637 T.DrawLatexNDC(0.5,0.5,'mean P8: %5.2F'%(h['p8']['nMultneu13p'+E].GetMean()))
638 h['mul'].Print('muonXsecP6P8_nMult'+E+'.png')
639 ut.bookCanvas(h,'Q2',E,900,600,1,1)
640 rc = h['Q2'].cd()
641 rc.SetLogy(1)
642 h['p6']['Qhist13p+500.0'].SetLineColor(ROOT.kBlue)
643 h['p6']['Qhist13p+500.0'].GetXaxis().SetRangeUser(0.,12.)
644 h['p6']['Qhist13p+500.0'].SetMaximum(1)
645 h['p6']['Qhist13p+500.0'].SetStats(0)
646 h['p6']['Qhist13p+500.0'].Draw()
647 h['p8']['Qhist13p500.0'].SetLineColor(ROOT.kGreen)
648 h['p8']['Qhist13p500.0'].SetStats(0)
649 h['p8']['Qhist13p500.0'].Draw('same')
650 h['Q2'].Print('muonXsecP6P8_Q2'+E+'.png')
651

◆ convertAscii2Root()

muonDis.convertAscii2Root (   fname,
  version = 2 
)

Definition at line 117 of file muonDis.py.

117def convertAscii2Root(fname,version=2):
118# convert fluka ascii table into root tntuple
119# /mnt/hgfs/microDisk/SND/MuonDis/unit30_Nm unit30_Pm
120# version 0: col 0 is the run number
121# col 1 is the event number
122# col 2 is the FLUKA particle type (10=muon+,11=muon-)
123# col 3 is the generation number (1 is a direct p-p collision product)
124# col 4 is the kinetic energy in GeV
125# col 5 is the statistical weight
126# col 6 is the x coordinate in cm (x=0 is as IP1, with the x-axis pointing outside the ring on the hor plane)
127# col 7 is the y coordinate in cm (y=0 is as IP1, with the y-axis pointing opposite to the gravity)
128# col 8 is the direction cosine wrt x-axis
129# col 9 is the direction cosine wrt y-axis
130# col 10 is the particle age in s (0 is the p-p collision time)
131# col 11 is useless
132
133# version 2: no change for columns 0 to 9, column 10 contains the z-coordinate (in cm) at the scoring plane
134 # Col 10: z coord (cm)
135 # Col 11: Last decay x cooord (cm)
136 # Col 12: Last decay y cooord (cm)
137 # Col 13: Last decay z cooord (cm)
138 # Col 14: Last decay ID
139 # Col 15: Last interaction x cooord (cm)
140 # Col 16: Last interaction y cooord (cm)
141 # Col 17: Last interaction z cooord (cm)
142 # Col 18: Last interaction ID
143
144 f = open(fname)
145 variables = "run:event:id:generation:E:w:x:y:px:py:t:z:pz"
146 # 0 1 2 3 4 5 6 7 8 9 10 11 12
147 fntuple = ROOT.TFile.Open(fname+'.root', 'RECREATE')
148 nt = ROOT.TNtupleD("nt","muon",variables)
149 for l in f.readlines():
150 if not l.find('#')<0: continue
151 tmp = l.replace('\n','').split(' ')
152 line = []
153 column = []
154 for x in tmp:
155 if x!='': line.append(float(x))
156# convert to more standard variables
157 column.append(line[0])
158 column.append(line[1])
159 column.append(flukaDict[line[2]])
160 column.append(line[3])
161 E = line[4] + muonMass
162 P = ROOT.TMath.Sqrt(E*E-muonMass*muonMass)
163 column.append(E)
164 column.append(line[5])
165 column.append(line[6])
166 column.append(line[7])
167 column.append(P*line[8])
168 column.append(P*line[9])
169 if version == 2:
170 column.append(0) # time had been dropped unfortunately
171 column.append(line[10])
172 else:
173 column.append(line[10]*u.s)
174 if version == 0: column.append(409.*u.m)
175 if version == 1: column.append(418.684*u.m + line[6]*ROOT.TMath.Tan(2.338/180.*ROOT.TMath.Pi()))
176 column.append(P*ROOT.TMath.Sqrt(1-line[8]**2-line[9]**2))
177 theTuple = array('d',column)
178 nt.Fill(theTuple)
179 fntuple.cd()
180 nt.Write()
181

◆ count_python_processes()

muonDis.count_python_processes (   macroName)

Definition at line 2234 of file muonDis.py.

2234def count_python_processes(macroName):
2235 username = pwd.getpwuid(os.getuid()).pw_name
2236 callstring = "ps -f -u " + username
2237# only works if screen is wide enough to print full name!
2238 status = subprocess.check_output(callstring,shell=True)
2239 n=0
2240 for x in status.decode().split('\n'):
2241 if not x.find(macroName)<0 and not x.find('python') <0: n+=1
2242 return n
2243

◆ drawMuon3D()

muonDis.drawMuon3D (   fname = 'muondEdX.root',
  hname = '3d',
  gDir = 'muMinusGeant4_0' 
)

Definition at line 1504 of file muonDis.py.

1504def drawMuon3D(fname='muondEdX.root',hname='3d',gDir='muMinusGeant4_0'):
1505 # muonDISfull.root h['xyz_Inter_'+str(pid)]
1506 ut.readHists(h,fname)
1507 ROOT.gStyle.SetPalette(ROOT.kRainBow)
1508 pal = ROOT.TColor.GetPalette()
1509 h['g'] = ROOT.TFile(gDir+'/geofile_full.conical.Ntuple-TGeant4.root')
1510 geoManager = h['g'] .FAIRGeom
1511 geoManager .SetNsegments(12)
1512 material = ROOT.TGeoMaterial("dummy")
1513 medium = ROOT.TGeoMedium('dummy',1,material)
1514 h['top'] = geoManager.GetTopVolume()
1515 top = h['top']
1516 nx,ny,nz = h[hname].GetXaxis().GetNbins(),h[hname].GetYaxis().GetNbins(),h[hname].GetZaxis().GetNbins()
1517 boxes = []
1518 hmax = pal.GetSize()/h[hname].GetMaximum()
1519 for ix in range(1,nx):
1520 for iy in range(1,ny):
1521 for iz in range(1,nz):
1522 C = h[hname].GetBinContent(ix,iy,iz)
1523 if not C>0: continue
1524 x,y,z = h[hname].GetXaxis().GetBinCenter(ix),h[hname].GetYaxis().GetBinCenter(iy),h[hname].GetZaxis().GetBinCenter(iz)
1525 bn = str(ix)+':'+str(iy)+':'+str(iz)
1526 box = geoManager.MakeBox(bn,medium,0.5,0.5,0.5)
1527 kc = int(pal.GetAt(int(C*hmax)))
1528 # print("color",int(C*hmax),kc)
1529 box.SetLineColor(kc)
1530 boxes.append(box)
1531 top.AddNode(box, 1 , ROOT.TGeoTranslation(x,y, z))
1532 for n in top.GetNodes():
1533 if n.GetName().find('TI')>0: n.SetVisibility(0)
1534 top.Draw('ogl')
1535

◆ energyOfDISevents()

muonDis.energyOfDISevents ( )

Definition at line 2254 of file muonDis.py.

2254def energyOfDISevents():
2255 for cycle in range(4):
2256 for run in range(1,11):
2257 for k in [0,1000]:
2258 prod = str(run+cycle*100+k)
2259 print('opening','muonDis_'+str(prod)+'.root')
2260 f = ROOT.TFile('muonDis_'+str(prod)+'.root')
2261 if not f: continue
2262 if not f.FindKey('DIS'): continue
2263 sTree = f.DIS
2264 for sTree in f.DIS:
2265 for P in sTree.Particles:
2266 hname = 'E_'+str(P.GetPdgCode())
2267 if not hname in h:
2268 ut.bookHist(h,hname,'Energy '+PDG.GetParticle(P.GetPdgCode()).GetName(),1000,0.,5000.)
2269 rc = h[hname].Fill(P.Energy())
2270 ut.writeHists(h,'energyOfDISevents.root')
2271# maybe better to make a first selection of events with n or KL above 100GeV
2272# cautious, n or KL can also be produced by other particles. Specially problematic,
2273# DIS events with no n or KL above 100GeV, but pion or protons make n or KL in transport!

◆ flukaMuons()

muonDis.flukaMuons (   version = 1,
  Plimit = False,
  withFaser = True 
)

Definition at line 1012 of file muonDis.py.

1012def flukaMuons(version=1,Plimit=False,withFaser=True):
1013 norm = normalization[version]
1014 path="/mnt/hgfs/microDisk/CERNBOX/SND@LHC/FLUKA/"
1015 pMin = 0
1016 if Plimit: pMin = 30.
1017 if version==0:
1018 fnames = {13:path+"muons_up/version0/unit30_Nm.root",-13:path+"muons_up/version0/unit30_Pm.root"}
1019 if version==1:
1020 if Plimit: pMin = 27.
1021 fnames = {13:path+"muons_up/version1/unit30_Nm.root",-13:path+"muons_up/version1/unit30_Pm.root"}
1022 if version==3:
1023 if Plimit: pMin = 27.
1024 fnames = {13:path+"muons_down/muons_VCdown_IR1-LHC.root"}
1025 boundaries()
1026 ut.bookCanvas(h,'TXY','muons',1800,650,2,1)
1027 ut.bookCanvas(h,'TE','muons',1800,650,2,1)
1028 ut.bookCanvas(h,'TW','muons',1800,650,2,1)
1029
1030 for mu in [13,-13]:
1031 ut.bookHist(h,'xySCO_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
1032 ut.bookHist(h,'xySND_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
1033 ut.bookHist(h,'xySNDX_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
1034 ut.bookHist(h,'E_'+str(mu), ';E [GeV]; N / 5 GeV',1000,0.,5000.)
1035 ut.bookHist(h,'W_'+str(mu), 'w',100,0.,0.15)
1036 for mu in fnames:
1037 fin = ROOT.TFile(fnames[mu])
1038 for sTree in fin.Get("nt"):
1039 m2cm = 1.
1040 if version == 0: m2cm = 100.
1041 P = ROOT.TVector3(sTree.px,sTree.py,sTree.pz)
1042 if P.Mag()<pMin: continue
1043 X = ROOT.TVector3(sTree.x,sTree.y,sTree.z*m2cm) # how could have this happened? used u.m for conversion!
1044 pid,w = int(sTree.id),sTree.w
1045 wLHC = norm*w
1046 rc = h['xySCO_'+str(pid)].Fill(X[0],X[1],wLHC)
1047 rc = h["E_"+str(pid)].Fill(P.Mag(),wLHC)
1048 rc = h['W_'+str(pid)].Fill(w)
1049 lam = ( SND_Z-X[2])/P[2]
1050 Xex = X+lam*P
1051 rc = h['xySND_'+str(pid)].Fill(Xex[0],Xex[1],wLHC)
1052 if abs(sTree.x)>50. or abs(sTree.y)>50.: rc = h['xySNDX_'+str(pid)].Fill(Xex[0],Xex[1],wLHC)
1053#
1054 for Y in ['xySCO_','xySND_','xySNDX_']:
1055 if version==0 and Y=='xySNDX_': continue
1056 k = 0
1057 for mu in [13,-13]:
1058 k+=1
1059 h[Y+str(mu)].SetStats(0)
1060 h[Y+str(mu)].SetMaximum(5.)
1061 if mu<0: h[Y+str(mu)].SetTitle("#mu^{+}")
1062 if mu>0: h[Y+str(mu)].SetTitle("#mu^{-}")
1063 tc = h['TXY'].cd(k)
1064 tc.SetRightMargin(0.2)
1065 h[Y+str(mu)].GetZaxis().SetTitleOffset(1.3)
1066 h[Y+str(mu)].Draw('colz')
1067 tc.Update()
1068 pal = h[Y+str(mu)].GetListOfFunctions()[0]
1069 pal.SetX1NDC(0.81)
1070 pal.SetX2NDC(0.85)
1071 if not Y.find('xySND')<0:
1072 for n in range(1,5):
1073 h['snd_'+str(n)].Draw('same')
1074 if withFaser: h['fas_'+str(n)].Draw('same')
1075 myPrint(h['TXY'],'FlukaMuons'+Y+'v'+str(version))
1076 k = 0
1077 for mu in [13,-13]:
1078 k+=1
1079 tc = h['TE'].cd(k)
1080 tc.SetLogy(1)
1081 if mu<0: h["E_"+str(mu)].SetTitle("#mu^{+}")
1082 if mu>0: h["E_"+str(mu)].SetTitle("#mu^{-}")
1083 h["E_"+str(mu)].SetStats(0)
1084 h["E_"+str(mu)].Draw()
1085 myPrint(h['TE'],'FlukaMuonsE_v'+str(version))
1086 k = 0
1087 for mu in [13,-13]:
1088 k+=1
1089 tc = h['TW'].cd(k)
1090 tc.SetLogy(1)
1091 if mu<0: h["W_"+str(mu)].SetTitle("#mu^{+}")
1092 if mu>0: h["W_"+str(mu)].SetTitle("#mu^{-}")
1093 h["W_"+str(mu)].SetStats(0)
1094 h["W_"+str(mu)].Draw()
1095 myPrint(h['TW'],'FlukaMuonsW_v'+str(version))
1096
1097 if not 'stats' in h: h['stats'] = {}
1098 h['stats'][version] = {}
1099 stats = h['stats'][version]
1100 for mu in fnames:
1101 stats[mu] = {}
1102 for z in ['snd','fas']:
1103 stats[mu][z] = []
1104 for Y in ['xySND_','xySNDX_']:
1105 hist = h[Y+str(mu)]
1106 xMin = hist.GetXaxis().FindBin(h[z+'_1'].GetX1())
1107 xMax = hist.GetXaxis().FindBin(h[z+'_1'].GetX2())
1108 yMin = hist.GetYaxis().FindBin(h[z+'_1'].GetY1())
1109 yMax = hist.GetYaxis().FindBin(h[z+'_2'].GetY1())
1110 if yMin>yMax:
1111 tmp = yMin
1112 yMin = yMax
1113 yMax = tmp
1114 sqcm = (xMax-xMin)*(yMax-yMin)
1115 stats[mu][z].append(hist.Integral(xMin,xMax,yMin,yMax))
1116 stats[mu][z].append(sqcm)
1117 stats[mu]['total'] = h['xySND_'+str(mu)].GetSumOfWeights()
1118 stats[mu]['1x1'] = h['xySND_'+str(mu)].GetSumOfWeights() - h['xySNDX_'+str(mu)].GetSumOfWeights()
1119 for mu in fnames:
1120 print("-------- "+PDG.GetParticle(mu).GetName()+" ----------")
1121 if version==1 and 0 in h['stats']:
1122 print('increae of rates, total: %5.2F absolute(1x1) %5.2F new/old=%5.2F'%(stats[mu]['1x1']/stats[mu]['total'],stats[mu]['1x1'],stats[mu]['1x1']/h['stats'][0][mu]['1x1']))
1123 newOvOld = (stats[mu]['snd'][0] -stats[mu]['snd'][1])/(h['stats'][0][mu]['snd'][0] -h['stats'][0][mu]['snd'][1])
1124 newOvOldTOT = stats[mu]['snd'][0]/h['stats'][0][mu]['snd'][0]
1125 print('increae of rates, SND: %5.2F absolute(1x1) %5.2F new/old=%5.2F new/old(tot)=%5.2F'%(stats[mu]['snd'][1]/stats[mu]['snd'][0],stats[mu]['snd'][0] -stats[mu]['snd'][1],newOvOld,newOvOldTOT ))
1126 newOvOld = (stats[mu]['fas'][0] -stats[mu]['fas'][1])/(h['stats'][0][mu]['fas'][0] -h['stats'][0][mu]['fas'][1])
1127 newOvOldTOT = stats[mu]['fas'][0]/h['stats'][0][mu]['fas'][0]
1128 print('increae of rates, FAS: %5.2F absolute(1x1) %5.2F new/old=%5.2F new/old(tot)=%5.2F'%(stats[mu]['fas'][1]/stats[mu]['fas'][0] ,stats[mu]['fas'][0] -stats[mu]['fas'][1],newOvOld,newOvOldTOT ))
1129 else:
1130 print('increae of rates, total: %5.2F absolute(1x1) %5.2F'%(h['xySNDX_'+str(mu)].GetSumOfWeights()/h['xySND_'+str(mu)].GetSumOfWeights(),
1131 h['xySND_'+str(mu)].GetSumOfWeights()-h['xySNDX_'+str(mu)].GetSumOfWeights()) )
1132 print('increae of rates, SND: %5.2F absolute(1x1) %5.2F'%(stats[mu]['snd'][1]/stats[mu]['snd'][0],stats[mu]['snd'][0] -stats[mu]['snd'][1] ))
1133 print('increae of rates, FAS: %5.2F absolute(1x1) %5.2F'%(stats[mu]['fas'][1]/stats[mu]['fas'][0] ,stats[mu]['fas'][0] -stats[mu]['fas'][1] ))
1134# sum of muons
1135 if version==1 and 0 in h['stats']:
1136 for z in ['snd','fas']:
1137 sumMu = h['stats'][1][13][z][0]+h['stats'][1][-13][z][0]
1138 increase = sumMu/(h['stats'][0][13][z][0]+h['stats'][0][-13][z][0])
1139 print('%s: %5.2F /cm2/s increase: %5.2F '%(z,sumMu/h['stats'][1][13][z][2],increase))
1140
1141

◆ G()

muonDis.G (   x)

Definition at line 45 of file muonDis.py.

45def G(x):
46 return 3./x**3*(x**2/2.-1.+ROOT.TMath.Exp(-x)*(1+x))

◆ getMasssq()

muonDis.getMasssq (   pid)

Definition at line 212 of file muonDis.py.

212def getMasssq(pid):
213 apid = abs(int(pid))
214 if not apid in masssq:
215 masssq[apid] = PDG.GetParticle(apid).Mass()**2
216 return masssq[apid]

◆ getPythia6CrossSec()

muonDis.getPythia6CrossSec (   nstat,
  pmom = [] 
)

Definition at line 227 of file muonDis.py.

227def getPythia6CrossSec(nstat,pmom=[]):
228 if len(pmom)==0:
229 pmom=[5.,10.,15.,20.,25.,50.,75.,100.,150.,200.,250,300.,400.,500.,750.,1000.,1500.,2500.,5000.,7500.,10000.]
230 mutype = {-13:'gamma/mu+',13:'gamma/mu-'}
231 target = ['p+','n']
232 for pid in mutype:
233 for t in target:
234 h['g_'+str(pid)+t] = ROOT.TGraph()
235 h['g_'+str(pid)+t].SetName('g_'+str(pid)+t)
236 np = 0
237 for P in pmom:
238 print('run pythia6 for ',pid,' on ',t,'with p=',P)
239 # Histograms.
240 tag = str(pid)+t+str(P)
241 ut.bookHist(h,"Qhist"+tag,"Q;[GeV]", 100, 0., 50.)
242 ut.bookHist(h,"pTehist"+tag,"pT of scattered muon;[GeV]", 100, 0., 50.)
243 ut.bookHist(h,"xhist"+tag,"x", 100, 0., 1.)
244 ut.bookHist(h,"yhist"+tag,"y", 100, 0., 1.)
245 ut.bookHist(h,"pTrhist"+tag,"pT of radiated parton; [GeV]", 100, 0., 50.)
246 ut.bookHist(h,"pTdhist"+tag,"ratio pT_parton/pT_muon", 100, 0., 5.)
247 ut.bookHist(h,"nMult"+tag,"particle multiplicity", 50, -0.5, 49.5)
248 ut.bookHist(h,"nMultcha"+tag,"charged particle multiplicity", 50, -0.5, 49.5)
249 ut.bookHist(h,"nMultneu"+tag,"neutral particle multiplicity", 50, -0.5, 49.5)
250#
251 myPythia = ROOT.TPythia6()
252 myPythia.SetMSEL(2) # msel 2 includes diffractive parts
253 myPythia.SetPARP(2,2) # To get below 10 GeV, you have to change PARP(2)
254 R = int(time.time()%900000000)
255 myPythia.SetMRPY(1,R)
256# stop pythia printout during loop
257 myPythia.SetMSTU(11, 11)
258 myPythia.Initialize('FIXT',mutype[pid],t,P)
259 for n in range(nstat):
260 myPythia.GenerateEvent()
261 if 0>1:
262 for i in range(1,myPythia.GetN()+1):
263 tmp = ROOT.Pythia8.Vec4(myPythia.GetP(i,1),myPythia.GetP(i,2),myPythia.GetP(i,3),myPythia.GetP(i,4))
264 print(i,myPythia.GetK(i,2),tmp.pAbs())
265 pProton = ROOT.Pythia8.Vec4(myPythia.GetP(2,1),myPythia.GetP(2,2),myPythia.GetP(2,3),myPythia.GetP(2,4))
266 peIn = ROOT.Pythia8.Vec4(myPythia.GetP(1,1),myPythia.GetP(1,2),myPythia.GetP(1,3),myPythia.GetP(1,4))
267 peOut = ROOT.Pythia8.Vec4(myPythia.GetP(3,1),myPythia.GetP(3,2),myPythia.GetP(3,3),myPythia.GetP(3,4))
268 pPhoton = peIn - peOut
269 # Q2, W2, Bjorken x, y.
270 Q2 = - pPhoton.m2Calc()
271 W2 = (pProton + pPhoton).m2Calc()
272 x = Q2 / (2. * (pProton * pPhoton))
273 y = (pProton * pPhoton) / (pProton * peIn)
274 h['Qhist'+tag].Fill( ROOT.TMath.Sqrt(Q2) )
275 h['xhist'+tag].Fill( x )
276 h['yhist'+tag].Fill( y )
277 h['pTehist'+tag].Fill( peOut.pT() )
278 # pT spectrum of partons being radiated in shower.
279 nMult = [0,0,0]
280 myPythia.Pyedit(2)
281 for i in range(1,myPythia.GetN()+1):
282 nMult[0]+=1
283 pidCode = myPythia.GetK(i,2) # gives the particle code
284 ch = myPythia.Pychge(pidCode)
285 if abs( ch )>0 : nMult[1]+=1
286 elif pidCode != 22 : nMult[2]+=1
287 tmp = ROOT.Pythia8.Vec4(myPythia.GetP(i,1),myPythia.GetP(i,2),myPythia.GetP(i,3),myPythia.GetP(i,4))
288 h['pTrhist'+tag].Fill( tmp.pT() )
289 h['pTdhist'+tag].Fill( tmp.pT() / peOut.pT() )
290 h['nMult'+tag].Fill(nMult[0])
291 h['nMultcha'+tag].Fill(nMult[1])
292 h['nMultneu'+tag].Fill(nMult[2])
293
294# very ugly procedure, but myPythia.GetPyint5() does not work!
295 myPythia.Pystat(0)
296 myPythia.Pystat(4)
297 xsec = readXsec(myPythia)
298 # print('sigma',xsec,ROOT.fixXsec(myPythia))
299 h['g_'+str(pid)+t].SetPoint(np,P,xsec)
300 np+=1
301 fout = ROOT.TFile('muDIScrossSec_Pythia6.root','RECREATE')
302 for pid in mutype:
303 for t in target:
304 h['g_'+str(pid)+t].Write()
305 for P in pmom:
306 tag = str(pid)+t+str(P)
307 for x in [ "nMult"+tag,"nMultcha"+tag,"nMultneu"+tag,"Qhist"+tag,"pTehist"+tag,
308 "xhist"+tag,"yhist"+tag,"pTrhist"+tag,"pTdhist"+tag]:
309 h[x].Write()
310

◆ getPythia8CrossSec()

muonDis.getPythia8CrossSec (   nstat,
  pmom = [] 
)

Definition at line 331 of file muonDis.py.

331def getPythia8CrossSec(nstat,pmom=[]):
332 if len(pmom)==0:
333 # pmom=[15.,20.,25.,50.,75.,100.,150.,200.,250,300.,400.,500.,750.,1000.,1500.,2500.,5000.,7500.,10000.]
334 pmom=[50.,75.,100.,150.,200.,250,300.,400.,500.,750.,1000.,1500.,2500.,5000.,7500.,10000.]
335 mutype = {-13:'mu+',13:'mu-'}
336 generators = {}
337 Q2min = 0. # 25.
338 for pid in mutype:
339 for g in ['p','n']:
340 h['g_'+str(pid)+g] = ROOT.TGraph()
341 h['g_'+str(pid)+g].SetName('g_'+str(pid)+g)
342 np = 0
343 for P in pmom:
344 generators[g]=ROOT.Pythia8.Pythia()
345 if g=='p': generators[g].settings.mode("Beams:idB", 2212)
346 else: generators[g].settings.mode("Beams:idB", 2112)
347 generators[g].settings.mode("Next:numberCount",options.heartbeat)
348 generators[g].settings.mode("Beams:frameType", 2)
349 generators[g].settings.parm("Beams:eA",P)
350 generators[g].settings.mode("Beams:idA", pid)
351 generators[g].settings.parm("Beams:eB",0.)
352 generators[g].readString("PDF:pSet = "+options.PDFpSet)
353#
354 # Neutral current (with gamma/Z interference).
355 generators[g].readString("WeakBosonExchange:ff2ff(t:gmZ) = on")
356#
357 # charged current.
358 generators[g].readString("WeakBosonExchange:ff2ff(t:W) = on")
359#
360 # Phase-space cut: minimal Q2 of process.
361 generators[g].settings.parm("PhaseSpace:Q2Min", Q2min)
362#
363 # Set dipole recoil on. Necessary for DIS + shower.
364 generators[g].readString("SpaceShower:dipoleRecoil = on")
365#
366 # Allow emissions up to the kinematical limit,
367 # since rate known to match well to matrix elements everywhere.
368 generators[g].readString("SpaceShower:pTmaxMatch = 2")
369#
370 # QED radiation off lepton not handled yet by the new procedure.
371 generators[g].readString("PDF:lepton = off")
372 generators[g].readString("TimeShower:QEDshowerByL = off")
373#
374 generators[g].init()
375 # Histograms.
376 tag = str(pid)+g+str(P)
377 ut.bookHist(h,"Qhist"+tag,"Q;[GeV]", 100, 0., 50.)
378 ut.bookHist(h,"pTehist"+tag,"pT of scattered muon;[GeV]", 100, 0., 50.)
379 ut.bookHist(h,"xhist"+tag,"x", 100, 0., 1.)
380 ut.bookHist(h,"yhist"+tag,"y", 100, 0., 1.)
381 ut.bookHist(h,"pTrhist"+tag,"pT of radiated parton; [GeV]", 100, 0., 50.)
382 ut.bookHist(h,"pTdhist"+tag,"ratio pT_parton/pT_muon", 100, 0., 5.)
383 ut.bookHist(h,"nMult"+tag,"particle multiplicity", 50, -0.5, 49.5)
384 ut.bookHist(h,"nMultcha"+tag,"charged particle multiplicity", 50, -0.5, 49.5)
385 ut.bookHist(h,"nMultneu"+tag,"neutral particle multiplicity", 50, -0.5, 49.5)
386#
387 for n in range(nstat):
388 rc = generators[g].next()
389 if not rc: continue
390 event = generators[g].event
391 # Four-momenta of proton, electron, virtual photon/Z^0/W^+-.
392 pProton = event[2].p()
393 peIn = event[1].p()
394 peOut = event[3].p()
395 pPhoton = peIn - peOut
396 # Q2, W2, Bjorken x, y.
397 Q2 = - pPhoton.m2Calc()
398 W2 = (pProton + pPhoton).m2Calc()
399 x = Q2 / (2. * (pProton * pPhoton))
400 y = (pProton * pPhoton) / (pProton * peIn)
401 h['Qhist'+tag].Fill( ROOT.TMath.Sqrt(Q2) )
402 h['xhist'+tag].Fill( x )
403 h['yhist'+tag].Fill( y )
404 h['pTehist'+tag].Fill( event[6].pT() )
405 # pT spectrum of partons being radiated in shower.
406 nMult = [0,0,0]
407 for i in range(event.size()):
408 # print(i,event[i].statusAbs(),event[i].id(), event[i].status())
409 if (event[i].status()>0):
410 nMult[0]+=1
411 if abs( event[i].chargeType() )>0 : nMult[1]+=1
412 elif event[i].id() != 22 : nMult[2]+=1
413 if (event[i].statusAbs() == 43):
414 h['pTrhist'+tag].Fill( event[i].pT() )
415 h['pTdhist'+tag].Fill( event[i].pT() / event[6].pT() )
416 h['nMult'+tag].Fill(nMult[0])
417 h['nMultcha'+tag].Fill(nMult[1])
418 h['nMultneu'+tag].Fill(nMult[2])
419 generators[g].stat()
420 xsec = ROOT.fixInfo(generators[g])
421 h['g_'+str(pid)+g].SetPoint(np,P,xsec)
422 np+=1
423 fout = ROOT.TFile('muDIScrossSec_Pythia8.root','RECREATE')
424 for pid in mutype:
425 for g in ['p','n']:
426 h['g_'+str(pid)+g].Write()
427 for P in pmom:
428 tag = str(pid)+g+str(P)
429 for x in ["nMult"+tag,"nMultcha"+tag,"nMultneu"+tag,"Qhist"+tag,"pTehist"+tag,
430 "xhist"+tag,"yhist"+tag,"pTrhist"+tag,"pTdhist"+tag]:
431 h[x].Write()
432
433

◆ makeMuDISEvents()

muonDis.makeMuDISEvents (   withElossFunction = False,
  nucleon = 'p+' 
)

Definition at line 448 of file muonDis.py.

448def makeMuDISEvents(withElossFunction=False,nucleon='p+'):
449# for energy loss:
450 if withElossFunction:
451 ut.readHists(h,"meanEloss.root")
452 eLoss = h['TCeloss'].FindObject('pol3').Clone('eLoss')
453#
454 myPythia = ROOT.TPythia6()
455 myPythia.SetMSEL(2) # msel 2 includes diffractive parts
456 myPythia.SetPARP(2,2) # To get below 10 GeV, you have to change PARP(2)
457 for kf in [211,321,130,310,3112,3122,3222,3312,3322,3334]:
458 kc = myPythia.Pycomp(kf)
459 myPythia.SetMDCY(kc,1,0)
460 R = int(time.time()%900000000)
461 myPythia.SetMRPY(1,R)
462 mutype = {-13:'gamma/mu+',13:'gamma/mu-'}
463# DIS event
464# run number
465# event number
466# incoming muon, id:px:py:pz:x:y:z:w
467# outgoing particles, id:px:py:pz
468 run = array('i', [0])
469 event = array('i', [0])
470 fout = ROOT.TFile('muonDis_'+str(options.run)+'.root','recreate')
471 dTree = ROOT.TTree('DIS','muon DIS')
472 dTree.Branch("run", run, "run/I")
473 dTree.Branch("event", event, "event/I")
474 iMuon = ROOT.TClonesArray("TParticle")
475 iMuonBranch = dTree.Branch("InMuon",iMuon,32000,-1)
476 dPart = ROOT.TClonesArray("TParticle")
477 dPartBranch = dTree.Branch("Particles",dPart,32000,-1)
478# read file with muons hitting concrete wall
479 fin = ROOT.TFile(options.muonIn) # id:px:py:pz:x:y:z:w
480 sTree = fin.Get("nt")
481 nTOT = sTree.GetEntries()
482 nEnd = min(nTOT,options.nStart + options.nEvents)
483# stop pythia printout during loop
484 myPythia.SetMSTU(11, 11)
485 print("start production ",options.nStart,nEnd)
486 nMade = 0
487 for k in range(options.nStart,nEnd):
488 rc = sTree.GetEvent(k)
489 run[0] = int(sTree.run)
490 event[0] = int(sTree.event)
491# make n events / muon
492 px,py,pz = sTree.px,sTree.py,sTree.pz
493 x,y,z = sTree.x,sTree.y,sTree.z-SND_Z
494 pid,w = int(sTree.id),sTree.w
495 if withElossFunction:
496 E = sTree.E - eLoss.Eval(sTree.E) # energy loss in 75m of rock
497 else:
498 E = sTree.E -27.
499 p = ROOT.TMath.Sqrt(px*px+py*py+pz*pz)
500 Peloss = ROOT.TMath.Sqrt(E*E-muonMass*muonMass)
501 if E < 5. : continue
502 t = sTree.t
503 # px=p*sin(theta)cos(phi),py=p*sin(theta)sin(phi),pz=p*cos(theta)
504 theta = ROOT.TMath.ACos(pz/p)
505 phi = ROOT.TMath.ATan2(py,px)
506 ctheta,stheta = ROOT.TMath.Cos(theta),ROOT.TMath.Sin(theta)
507 cphi,sphi = ROOT.TMath.Cos(phi),ROOT.TMath.Sin(phi)
508 escale = Peloss/p
509 muPart = ROOT.TParticle(pid,0,0,0,0,0,px*escale,py*escale,pz*escale,E,x,y,z,t)
510 muPart.SetWeight(w)
511 myPythia.Initialize('FIXT',mutype[pid],nucleon,E)
512 for n in range(options.nMult):
513 dPart.Clear()
514 iMuon.Clear()
515 muPart_replica = ROOT.TParticle(muPart)
516 muPart_TCA = iMuon.ConstructedAt(0)
517 ROOT.std.swap(muPart_replica, muPart_TCA)
518 myPythia.GenerateEvent()
519# remove all unnecessary stuff
520 myPythia.Pyedit(2)
521 for itrk in range(1,myPythia.GetN()+1):
522 did = myPythia.GetK(itrk,2)
523 dpx,dpy,dpz = rotate(ctheta,stheta,cphi,sphi,myPythia.GetP(itrk,1),myPythia.GetP(itrk,2),myPythia.GetP(itrk,3))
524 psq = dpx**2+dpy**2+dpz**2
525 E = ROOT.TMath.Sqrt(getMasssq(did)+psq)
526 part = ROOT.TParticle(did,0,0,0,0,0,dpx,dpy,dpz,E,x,y,z,t)
527 part.SetWeight(w/float(options.nMult))
528# copy to branch
529 nPart = dPart.GetEntries()
530 if dPart.GetSize() == nPart: dPart.Expand(nPart+10)
531 part_TCA = dPart.ConstructedAt(nPart)
532 ROOT.std.swap(part, part_TCA)
533 nMade+=1
534 if nMade%options.heartbeat==0: print('made so far ',options.run,' :',nMade,time.ctime())
535 dTree.Fill()
536 fout.cd()
537 dTree.Write()
538 myPythia.SetMSTU(11, 6)
539 print("finished ",options.run)

◆ missing3Dproj()

muonDis.missing3Dproj (   hist,
  ymin,
  ymax 
)

Definition at line 2309 of file muonDis.py.

2309def missing3Dproj(hist,ymin,ymax):
2310 h[hist+'XZ']=h[hist].Project3D('xz')
2311 h[hist+'XZ'].SetName(hist+'XZ')
2312 h[hist+'XZ'].Reset()
2313 for ix in range(1,h[hist].GetNbinsX()+1):
2314 for iz in range(1,h[hist].GetNbinsZ()+1):
2315 N = 0
2316 for iy in range(ymin,ymax): N+=h[hist].GetBinContent(ix,iy,iz)
2317 h[hist+'XZ'].SetBinContent(iz,ix,N)
2318
2319

◆ muInterGeant4()

muonDis.muInterGeant4 (   version = 2,
  njobs = 100 
)

Definition at line 1142 of file muonDis.py.

1142def muInterGeant4(version=2,njobs=100):
1143 redfac3d = 10
1144 norm = normalization[version]
1145 uninterestingParticles = [11,-11,-12,12,14,-14,13,-13]
1146 ut.bookHist(h,"xyz_muInter_", 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
1147 ut.bookHist(h,"xyz_origin_", 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
1148
1149 if version == 2: g = ROOT.TFile('muGeant4_0/geofile_full.conical.Ntuple-TGeant4.root')
1150 else: g = ROOT.TFile('muMinusGeant4_0/geofile_full.conical.Ntuple-TGeant4.root')
1151 geoManager = g.FAIRGeom
1152 for subjob in range(njobs):
1153 if version == 2: muons = ['muGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root']
1154 else: muons = ['muMinusGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root','muPlusGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root']
1155 pids = {}
1156 for fname in muons:
1157 f = ROOT.TFile(fname)
1158 ROOT.gROOT.cd()
1159 nEv = -1
1160 for sTree in f.Get("cbmsim"):
1161 nEv+=1
1162 muon = sTree.MCTrack[0]
1163 w = norm*muon.GetWeight()
1164 muID = muon.GetPdgCode()
1165 for p in sTree.vetoPoint:
1166 pid = p.PdgCode()
1167 if not pid in pids: pids[pid]=0
1168 pids[pid]+=1
1169 if pid in uninterestingParticles: continue
1170 if not "xyz_origin_"+str(pid) in h:
1171 for x in ["xyz_muInter_","xyz_origin_"]:
1172 h[x+str(pid)]=h[x].Clone(x+str(pid))
1173 track = p.GetTrackID()
1174 exitPoint = p.LastPoint()
1175 rc = h["xyz_origin_"+str(pid)].Fill(exitPoint[0],exitPoint[1],exitPoint[2],w)
1176 while track > 0:
1177 mother = sTree.MCTrack[track].GetMotherId()
1178 if mother == 0: break
1179 track = mother
1180 if track >0:
1181 m = sTree.MCTrack[track]
1182 rc = h["xyz_muInter_"+str(pid)].Fill(m.GetStartX(),m.GetStartY(),m.GetStartZ(),w)
1183
1184

◆ muondEdX()

muonDis.muondEdX (   version = 2,
  njobs = 100,
  path = '',
  withFaser = False,
  plotOnly = True 
)

Definition at line 1185 of file muonDis.py.

1185def muondEdX(version=2,njobs=100,path='',withFaser=False, plotOnly=True):
1186# version 1 /home/truf/ubuntu-1710/ship-ubuntu-1710-48/SND/MuonDis/
1187
1188 # python -i $SNDBUILD/sndsw/shipLHC/run_simSND.py --PG --Estart 10 --Eend 5000 --EVz -7100 --EVx -30 --EVy 40 --pID 13 -n 100000 --FastMuon
1189 # python -i $SNDBUILD/sndsw/shipLHC/run_simSND.py --PG --Estart 10 --Eend 5000 --EVz -7100 --EVx -30 --EVy 40 --pID -13 -n 100000 --FastMuon
1190 # python $SNDSW_ROOT/shipLHC/run_simSND.py -f unit30_Nm.root --Ntuple --FastMuon --output muMinusGeant4 -n 999999999
1191 # python $SNDSW_ROOT/shipLHC/run_simSND.py -f unit30_Pm.root --Ntuple --FastMuon --output muPlusGeant4 -n 999999999
1192 neuPartList = [22,130,310,2112,-2112,3122,-3122,3322,-3322]
1193 if not plotOnly:
1194 norm = normalization[version]
1195
1196 ut.bookHist(h,'eloss','energyLoss', 500,0., 5000.,500,0.,500.)
1197 ut.bookHist(h,'dx','multiple scattering x', 500,0., 5000.,1000,-0.02,0.02)
1198 ut.bookHist(h,'dy','multiple scattering y', 500,0., 5000.,1000,-0.02,0.02)
1199 ut.bookHist(h,'oE','original energy', 500,0.,5000.)
1200 ut.bookHist(h,'oL','flight path', 100,0.,10000.)
1201 ut.bookHist(h,'3d','exit points',200,-100.,100.,200,-100.,100.,600,-500.,100.)
1202 ut.bookHist(h,'z','last z ',1000,-1000.,1000.)
1203 ut.bookHist(h,'SNDmuP', 'SND entry points',200,-100.,100.,200,-100.,100.)
1204 ut.bookHist(h,'SNDmuM','SND entry points',200,-100.,100.,200,-100.,100.)
1205 ut.bookHist(h,'SNDmuP_E', 'SND muon enery',200,0.,5000.)
1206 ut.bookHist(h,'SNDmuM_E','SND muon energy',200,0.,5000.)
1207 ut.bookHist(h,'allmuE_13', 'all muon enery',200,0.,5000.)
1208 ut.bookHist(h,'allmuE_-13','all muon energy',200,0.,5000.)
1209 ut.bookHist(h,'SNDMuFiltermuP', 'SND entry points',200,-100.,100.,200,-100.,100.)
1210 ut.bookHist(h,'SNDMuFiltermuM','SND entry points',200,-100.,100.,200,-100.,100.)
1211 ut.bookHist(h,'xy_13', 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
1212 ut.bookHist(h,'xy_-13', 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
1213
1214 for p in neuPartList: ut.bookHist(h,"E_"+str(p),'Energy', 500,0.,5000.)
1215
1216 if version == 3: g = ROOT.TFile(path+'muGeant4_VCdown_0/geofile_full.conical.Ntuple-TGeant4.root')
1217 elif version == 2: g = ROOT.TFile(path+'muGeant4_0/geofile_full.conical.Ntuple-TGeant4.root')
1218 else: g = ROOT.TFile(path+'muMinusGeant4_0/geofile_full.conical.Ntuple-TGeant4.root')
1219 geoManager = g.FAIRGeom
1220 for subjob in range(njobs):
1221 if version == 3: muons = [path+'muGeant4_VCdown_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root']
1222 elif version == 2: muons = [path+'muGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root']
1223 else: muons = [path+'muMinusGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root',path+'muPlusGeant4_'+str(subjob)+'/ship.conical.Ntuple-TGeant4.root']
1224 for fname in muons:
1225 h[fname] = ROOT.TFile(fname)
1226 ROOT.gROOT.cd()
1227 nEv = -1
1228 for sTree in h[fname].Get("cbmsim"):
1229 nEv+=1
1230 muon = sTree.MCTrack[0]
1231 w = norm*muon.GetWeight()
1232 muID = muon.GetPdgCode()
1233 if muon.GetPdgCode()==0:
1234 print('something strange here PdgCode=0',fname)
1235 rc = h['oE'].Fill(muon.GetEnergy(),w)
1236 for x in sTree.MCTrack:
1237 pid = x.GetPdgCode()
1238 if pid in neuPartList: rc = h["E_"+str(pid)].Fill(x.GetEnergy())
1239 trajectories = {}
1240 for p in sTree.vetoPoint:
1241 track = p.GetTrackID()
1242 if track < 0: continue # trajectory not stored
1243 if track != 0: continue # only take original muon for the moment
1244
1245 if not track in trajectories: trajectories[track]={}
1246 traj = trajectories[track]
1247 end = p.LastPoint()
1248 middle = ROOT.TVector3(p.GetX(),p.GetY(),p.GetZ())
1249 first = 2*middle - end
1250 Pin = ROOT.Math.PxPyPzMVector(p.GetPx(),p.GetPy(),p.GetPz(),muonMass)
1251 Pout = ROOT.Math.PxPyPzMVector(p.LastMom()[0],p.LastMom()[1],p.LastMom()[2],muonMass)
1252 if not Pout.Z()>0:
1253 enterCave = False
1254 if sTree.ScifiPoint.GetEntries()>0 or sTree.MuFilterPoint.GetEntries()>0 and end.z() < -24.5: # otherwise muon stops in or after SND
1255 print('something strange here, pout.z<0',nEv,fname,sTree.ScifiPoint.GetEntries(),sTree.MuFilterPoint.GetEntries())
1256 p.Print()
1257 else:
1258 test = end + 1./Pout.Z()*p.LastMom()
1259 node = geoManager.FindNode(test[0],test[1],test[2])
1260 enterCave = node.GetName() .find('Vrock') <0 # particle leaving concrete and entering rock
1261 traj[first[2]]=[first[0],first[1],Pin,False]
1262 traj[end[2]]=[end[0],end[1],Pout,enterCave]
1263 for track in trajectories:
1264 traj = trajectories[track]
1265 zPos = list(traj.keys())
1266 zPos.sort()
1267 T = [ROOT.TGraph(),ROOT.TGraph()]
1268 k = 0
1269 for z in zPos:
1270 for j in range(2):
1271 T[j].SetPoint(k,z,traj[z][j])
1272 k+=1
1273 xex,yex = T[0].Eval(0),T[1].Eval(0)
1274 if xex==0 and yex==0:
1275 T[0].Print()
1276 T[1].Print()
1277 print(nEv,fname,track,traj,zPos)
1278 rc = h['xy_'+str(muID)].Fill(xex,yex,w)
1279 first = True
1280 for z in zPos:
1281 if not traj[z][3]: continue
1282 oEnergy = muon.GetEnergy()
1283 P = traj[z][2]
1284 if first:
1285 h['allmuE_'+str(muID)].Fill(P.E())
1286 first = False
1287 eloss = muon.GetEnergy() - P.E()
1288 rc = h['eloss'].Fill(oEnergy,eloss,w)
1289 start = ROOT.TVector3(sTree.MCTrack[0].GetStartX(),sTree.MCTrack[0].GetStartY(),sTree.MCTrack[0].GetStartZ())
1290 L = end-start
1291 rc = h['oL'].Fill(L.Mag(),w)
1292 if P.Z()>0:
1293 dAlpha = muon.GetPx()/muon.GetPz() - P.X()/P.Z()
1294 rc = h['dx'].Fill(oEnergy,dAlpha)
1295 dAlpha = muon.GetPy()/muon.GetPz()- P.Y()/P.Z()
1296 rc = h['dy'].Fill(oEnergy,dAlpha)
1297 rc = h['3d'].Fill(traj[z][0],traj[z][1],z,w)
1298 break
1299 hitScifi = False
1300 for p in sTree.ScifiPoint:
1301# rate of muons entering SND, count only one hit
1302 omu = p.GetTrackID()
1303 if omu != 0 or abs(p.PdgCode()) != 13: continue
1304 hname = 'SNDmuM'
1305 if p.PdgCode()<0: hname = 'SNDmuP'
1306 rc = h[hname].Fill(p.GetX(),p.GetY(),w)
1307 mom = ROOT.TVector3(p.GetPx(),p.GetPy(),p.GetPz())
1308 rc = h[hname+'_E'].Fill(mom.Mag(),w)
1309 hitScifi=True
1310 break
1311 hitMuFilter = False
1312 for p in sTree.MuFilterPoint:
1313# rate of muons entering SND, count only one hit
1314 omu = p.GetTrackID()
1315 if omu != 0 or abs(p.PdgCode()) != 13: continue
1316 hname = 'SNDMuFiltermuM'
1317 if p.PdgCode()<0: hname = 'SNDMuFiltermuP'
1318 rc = h[hname].Fill(p.GetX(),p.GetY(),w)
1319 hitMuFilter=True
1320 break
1321# go for first point, exit from concrete
1322 noHit = True
1323 for p in sTree.ScifiPoint:
1324 if abs(p.PdgCode()) != 13: continue
1325 if p.GetDetectorID()==47:
1326 noHit = False
1327 break
1328 if noHit: continue
1329#
1330 ut.writeHists(h,'muondEdX_'+str(version)+'.root')
1331 else:
1332 ut.readHists(h,'muondEdX_'+str(version)+'.root')
1333
1334 for det in ['SND','SNDMuFilter']:
1335 S = h[det+'muP'].Clone('S')
1336 S.Add(h[det+'muM'])
1337 hname = det+'muP'
1338 nx,ny = h[hname].GetXaxis().GetNbins(),h[hname].GetYaxis().GetNbins()
1339 projx = S.ProjectionX()
1340 projy = S.ProjectionY()
1341 for ix in range(1,nx+1):
1342 if projx.GetBinContent(ix)>0:
1343 xmin = projx.GetBinCenter(ix)-projx.GetBinWidth(ix)/2.
1344 break
1345 for ix in range(nx,0,-1):
1346 if projx.GetBinContent(ix)>0:
1347 xmax = projx.GetBinCenter(ix)+projx.GetBinWidth(ix)/2.
1348 break
1349 for iy in range(1,ny+1):
1350 if projy.GetBinContent(iy)>0:
1351 ymin = projy.GetBinCenter(iy)-projy.GetBinWidth(iy)/2.
1352 break
1353 for iy in range(ny,0,-1):
1354 if projy.GetBinContent(iy)>0:
1355 ymax = projy.GetBinCenter(iy)+projy.GetBinWidth(iy)/2.
1356 break
1357 sqcm = (xmax-xmin)*(ymax-ymin)
1358 muP = h[det+'muP'].GetSumOfWeights()
1359 muM = h[det+'muM'].GetSumOfWeights()
1360 print("%s square: %5.2F,%5.2F,%5.2F,%5.2F"%(det,xmin,xmax,ymin,ymax))
1361 print(" mu+ rate = %5.2F Hz %5.2F Hz/cm2"%(muP,muP/sqcm))
1362 print(" mu- rate = %5.2F Hz %5.2F Hz/cm2"%(muM,muM/sqcm))
1363 print(" sum rate = %5.2F Hz %5.2F Hz/cm2"%(muM+muP,(muM+muP)/sqcm))
1364
1365 ut.bookCanvas(h,'TCeloss','eloss',900,600,1,1)
1366 tc = h['TCeloss'].cd()
1367 h['eloss_mean'] = h['eloss'].ProfileX('mean',1,-1,'g')
1368 h['eloss_mean'].GetXaxis().SetRangeUser(10,5000)
1369 fname = 'pol3'
1370 h['eloss_mean'].Fit(fname,'S','',20.,5000.)
1371 tc.SetLogx(1)
1372 h['eloss_mean'].SetTitle('; incoming momentun [GeV/c] ; mean energy loss [GeV]')
1373 h['eloss_mean'].Draw('hist')
1374 fun = h['eloss_mean'].GetFunction(fname)
1375 fun.Draw('same')
1376 tc.Update()
1377 stats = h['eloss_mean'].FindObject('stats')
1378 stats.SetOptFit(111)
1379 stats.SetOptStat(0)
1380 stats.SetX1NDC(0.2)
1381 stats.SetY1NDC(0.5)
1382 stats.SetX2NDC(0.6)
1383 stats.SetY2NDC(0.75)
1384 tc.Update()
1385 myPrint(tc,'meanEloss_'+str(version))
1386#
1387 ut.bookCanvas(h,'eff','efficiency',900,600,1,1)
1388 tc = h['eff'].cd()
1389 eff = h['eloss'].ProjectionX()
1390 eff.Divide(h['oE'])
1391 eff.SetStats(0)
1392 eff.SetTitle("; incoming momentun [GeV/c] ; efficiency")
1393 eff.GetXaxis().SetRangeUser(0.,500.)
1394 eff.Draw()
1395 myPrint(tc,'efficiency_'+str(version))
1396#
1397 ut.bookCanvas(h,'TCMS','multiple scattering',900,600,1,1)
1398 tc = h['TCMS'].cd()
1399 msDx = h['dx'].ProfileX('msDx',1,-1,'s')
1400 N = msDx.GetNbinsX()
1401 Xmin = msDx.GetBinCenter(1)-msDx.GetBinWidth(1)
1402 Xmax = msDx.GetBinCenter(N)+msDx.GetBinWidth(N)
1403 ut.bookHist(h,'msDxE','msDxE',N,Xmin,Xmax)
1404 msDxE = h['msDxE']
1405 # angular distributions are distorted due to different paths and different integrated material
1406 for j in range(1,msDx.GetNbinsX()+1):
1407 tmp = h['dy'].ProjectionY('tmp',j,j)
1408 rms = 0
1409 Erms = 0
1410 if tmp.GetEntries()>20:
1411 rc = tmp.Fit('gaus','SQL')
1412 fitResult = rc.Get()
1413 rms = fitResult.Parameter(2)
1414 Erms = fitResult.ParError(2)
1415 msDxE.SetBinContent(j,rms)
1416 msDxE.SetBinError(j,Erms)
1417 msDxE.SetTitle("; incoming momentun [GeV/c] ; MS angle [rad]")
1418 msDxE.SetMaximum(0.02)
1419 msDxE.SetMinimum(0.0001)
1420 tc.SetLogy(1)
1421 msDxE.GetXaxis().SetRangeUser(0.,1000.)
1422 msDxE.Draw('hist')
1423 rc = msDxE.Fit(h['MS'],'S','',20.,400.)
1424 tc.Update()
1425 stats = msDxE.FindObject('stats')
1426 stats.SetOptFit(111)
1427 h['MS'].Draw('same')
1428 myPrint(h['TCMS'],'multipleScattering_'+str(version))
1429# rates
1430 boundaries()
1431 ut.bookCanvas(h,'TXY','muons',1800,650,2,1)
1432 k = 0
1433 for mu in [13,-13]:
1434 k+=1
1435 hist = h['xy_'+str(mu)]
1436 hist.SetStats(0)
1437 hist.SetMaximum(1.)
1438 if mu<0: hist.SetTitle("#mu^{+}")
1439 if mu>0: hist.SetTitle("#mu^{-}")
1440 tc = h['TXY'].cd(k)
1441 tc.SetRightMargin(0.2)
1442 hist.GetZaxis().SetTitleOffset(1.3)
1443 hist.Draw('colz')
1444 tc.Update()
1445 pal = hist.GetListOfFunctions()[0]
1446 pal.SetX1NDC(0.81)
1447 pal.SetX2NDC(0.85)
1448 for n in range(1,5):
1449 h['snd_'+str(n)].Draw('same')
1450 if withFaser: h['fas_'+str(n)].Draw('same')
1451 myPrint(h['TXY'],'Geant4Muonsxy_v'+str(version))
1452
1453 h['stats'] = {}
1454 stats = h['stats']
1455 for mu in [13,-13]:
1456 stats[mu] = {}
1457 for z in ['snd','fas']:
1458 stats[mu][z] = []
1459 hist = h['xy_'+str(mu)]
1460 xMin = hist.GetXaxis().FindBin(h[z+'_1'].GetX1())
1461 xMax = hist.GetXaxis().FindBin(h[z+'_1'].GetX2())
1462 yMin = hist.GetYaxis().FindBin(h[z+'_1'].GetY1())
1463 yMax = hist.GetYaxis().FindBin(h[z+'_2'].GetY1())
1464 if yMin>yMax:
1465 tmp = yMin
1466 yMin = yMax
1467 yMax = tmp
1468 sqcm = (xMax-xMin)*(yMax-yMin)
1469 stats[mu][z].append(hist.Integral(xMin,xMax,yMin,yMax))
1470 stats[mu][z].append(sqcm)
1471 print('%s, %s: total = %5.2F Hz %5.2F Hz/cm2 '%(mu,z,stats[mu][z][0],stats[mu][z][0]/stats[mu][z][1]))
1472
1473# interaction rate in target:
1474 ut.bookCanvas(h,'Tmuons','muons',1600,900,2,1)
1475 fin = ROOT.TFile('muDIScrossSec.root')
1476 ROOT.gROOT.cd()
1477 L = {'13':'SNDmuM','-13':'SNDmuP'}
1478 for pid in L:
1479 h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid)
1480 hname = L[pid]+'_inter'
1481 if pid=='13': h[L[pid]+'_E'].SetLineColor(ROOT.kRed)
1482 else: h[L[pid]+'_E'].SetLineColor(ROOT.kBlue)
1483 h[L[pid]+'_E'].SetTitle(';E [GeV]; N/s')
1484 h[L[pid]+'_E'].SetStats(0)
1485 h[hname]=h[L[pid]+'_E'].Clone(hname)
1486 for k in range(1,h[hname].GetNbinsX()+1):
1487 muonEnergy = h[L[pid]+'_E'].GetBinCenter(k)
1488 wDis = 5.9*19.3 / 1.67E-24 * h['g_'+str(pid)].Eval(muonEnergy) * 1E-27
1489 h[hname].SetBinContent(k,wDis*h[L[pid]+'_E'].GetBinContent(k))
1490 ut.makeIntegralDistrib(h,L[pid]+'_E')
1491 ut.makeIntegralDistrib(h,hname)
1492 tc = h['Tmuons'].cd(1)
1493 tc.SetLogy(1)
1494 h['I-SNDmuM_E'].Draw('hist')
1495 h['I-SNDmuP_E'].Draw('histsame')
1496 tc = h['Tmuons'].cd(2)
1497 tc.SetLogy(1)
1498 h['I-SNDmuM_inter'].Draw('hist')
1499 h['I-SNDmuP_inter'].Draw('histsame')
1500# 1fb-1 requires 1E5 sec, 25fb-1 = 2.5E6 seconds, -> 0.25 million mu interactions with E>200GeV.
1501 myPrint(h['Tmuons'],'Geant4MuonE_v'+str(version))
1502
1503# drawMuon3D(fname='muonDISfull.root',hname='3d')

◆ muonDISfull()

muonDis.muonDISfull (   cycle = 0,
  sMin = 0,
  sMax = 200,
  rMin = 1,
  rMax = 11,
  path = '/eos/experiment/ship/user/truf/SND/muonDis/',
  debug = 0,
  version = 1,
  pythia6 = True 
)

Definition at line 1571 of file muonDis.py.

1571def muonDISfull(cycle = 0, sMin=0,sMax=200,rMin=1,rMax=11,path = '/eos/experiment/ship/user/truf/SND/muonDis/',debug=0,version=1,pythia6=True):
1572#
1573 muRange = [0,1000]
1574 geofile = 'geofile_full.conical.muonDIS-TGeant4.root'
1575 geofileExample = 'run_1_1/'+geofile
1576 datafile = 'ship.conical.muonDIS-TGeant4.root'
1577 if not pythia6:
1578 muRange = ['muMinusGeant4','muPlusGeant4']
1579 geofile ='geofile_full.conical.Ntuple-TGeant4.root'
1580 geofileExample = 'muMinusGeant4_0/'+geofile
1581 datafile = '/ship.conical.Ntuple-TGeant4.root'
1582 if cycle ==100: datafile = '/ship.conical.Ntuple-TGeant4_boost100.0.root'
1583
1584 redfac3d = 10
1585 norm = normalization[version]
1586 fin = ROOT.TFile('muDIScrossSec.root')
1587 ROOT.gROOT.cd()
1588 for pid in ['13','-13']:
1589 h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid)
1590 ut.bookHist(h,'wDIS_'+pid,'muon DIS probability',100,0.,0.1)
1591 for z in ['0','inter']:
1592 ut.bookHist(h,'inMu_'+z+str(pid),'muon Energy vs #direct hits',200,0.,5000.,20,0.0,19.5)
1593 ut.bookHist(h,"xyz_mu_"+z+str(pid), 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
1594 neuPartList = [22,130,310,2112,-2112,3122,-3122,3322,-3322]
1595 ut.bookHist(h,'fullStats','statistics',2000,-0.5,1999.5)
1596 ut.bookHist(h,'pidsDict','pids dictionary',100,-0.5,99.5)
1597 n = 1
1598 for x in pidsDict:
1599 h['pidsDict'].SetBinContent(n,x)
1600 n+=1
1601
1602 for p in neuPartList:
1603 for o in ['','_secondary']:
1604 ut.bookHist(h,'inE_'+str(p)+o,'energy',500,0.,5000.,600,-500.,100.) # end vertex of primary neutral particle
1605 ut.bookHist(h,'inE_Concrete_'+str(p)+o,'energy',500,0.,5000.,600,-500.,100.) # end vertex of primary neutral particle if in concrete
1606 ut.bookHist(h,"E_"+str(p),'Energy', 500,0.,5000.,200,-30.,170.)
1607 ut.bookHist(h,"E_veto_"+str(p),'Energy', 500,0.,5000.,200,-30.,170.)
1608 ut.bookHist(h,"E_vetoParticle_"+str(p),'Energy', 500,0.,50.,20,-0.5,19.5)
1609 ut.bookHist(h,"Eprim_"+str(p),'Energy', 500,0.,5000.,200,-30.,170.)
1610 ut.bookHist(h,"Eprim_veto_"+str(p),'Energy', 500,0.,5000.,200,-30.,170.)
1611 ut.bookHist(h,"startZ_"+str(p),'z',200,-100.,300.)
1612 ut.bookHist(h,"xyz_Inter_"+str(p), 'x/y/z ',200,-100.,100.,200,-100.,100.,int(400/redfac3d),-100.,100.)
1613 ut.bookHist(h,"xyzE100_Inter_"+str(p), 'x/y/z ',200,-100.,100.,200,-100.,100.,int(400/redfac3d),-100.,100.)
1614 ut.bookHist(h,"xyzVeto_Inter_"+str(p), 'x/y/z ',200,-100.,100.,200,-100.,100.,int(400/redfac3d),-100.,100.)
1615 ut.bookHist(h,"xyz_muInter_"+str(p), 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
1616 ut.bookHist(h,"xyz_origin_"+str(p), 'x/y /z',200,-100.,100.,200,-100.,100.,int(600/redfac3d) ,-500.,100.)
1617 ut.bookHist(h,"z_veto_"+str(p), 'z veto vs z neutral',100,-50.,50.,100,-50.,50.)
1618 ut.bookHist(h,"xy_z-30_veto_"+str(p), 'xy at z=-30cm',200,-100.,100.,200,-100.,100.)
1619 ut.bookHist(h,"xy_z-30_"+str(p), 'xy at z=-30cm',200,-100.,100.,200,-100.,100.)
1620
1621# python -i $SNDBUILD/sndsw/macro/run_simScript.py --shiplhc -f muonDis_1108.root --MuDIS --output run_1108 -n 99999999
1622# prod = str(run+cycle*100+k), k=0 or 1000, mu+ or mu-, number of cycles: 10 events per incoming muon in each cycle)
1623 if path.find('eos')<0: g = ROOT.TFile.Open(path+geofileExample)
1624 else: g = ROOT.TFile.Open(os.environ['EOSSHIP']+path+geofileExample)
1625 geoManager = g.FAIRGeom
1626 if path.find('eos')<0:
1627 topDir = os.listdir(path)
1628 else:
1629 topDir = str( subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls -l "+path,shell=True) )
1630 # for cycle in [options.nMult ]:
1631 for run in range(rMin,rMax):
1632 for k in muRange:
1633 for subjob in range(sMin,sMax):
1634 if pythia6:
1635 if options.Emin==0: prod = 'run_'+str(run+cycle*100+k)+'_'+str(subjob)
1636 else: prod = "ecut"+str(options.Emin)+'_run_'+str(run+cycle*100+k)+'_'+str(subjob)
1637 else:
1638 prod = k+"_"+str(subjob)
1639 if path.find('eos')<0:
1640 if not prod in topDir:
1641 print('prod not found ',prod)
1642 continue
1643 if not geofile in os.listdir(path+prod):
1644 print('no geofile found ',path+prod)
1645 continue
1646 else:
1647 if not "/"+prod in topDir:
1648 print('prod not found on eos',prod)
1649 continue
1650 temp = subprocess.check_output("xrdfs "+os.environ['EOSSHIP']+" ls -l "+path+prod,shell=True)
1651 if not geofile in str(temp):
1652 print('no geofile found on eos',path+prod)
1653 continue
1654 if pythia6: rc = h['fullStats'].Fill(run+cycle*100+k)
1655 if path.find('eos')<0: h['f'] = ROOT.TFile.Open(path+prod+datafile)
1656 else: h['f'] = ROOT.TFile.Open(os.environ['EOSSHIP']+path+prod+datafile)
1657 nEv = -1
1658 for sTree in h['f'].Get("cbmsim"):
1659 nEv+=1
1660 if nEv%1000 == 0: print('N ',nEv,prod)
1661 muon = sTree.MCTrack[0]
1662 muonEnergy = muon.GetEnergy()
1663 mupid = muon.GetPdgCode()
1664 muInterVx = ROOT.TVector3(muon.GetStartX(),muon.GetStartY(),muon.GetStartZ())
1665 if pythia6:
1666 # weights on muonDIs files: incoming muon has original FLUKA weight, outgoing particles FLUKA weight/options.nMult
1667 # unfortunately, weight of outgoing particles is overwritten by MuDISGenerator with density along trajectory as weight, g/cm^2
1668 nMult = 10
1669 wLHC = norm*muon.GetWeight()/float(nMult)
1670 wInter = sTree.MCTrack[2].GetWeight() # put density along trajectory as weight, g/cm^2
1671 # fast MC used NinteractionLength * 97.5 g cm-2 interaction length of concrete --> wInter
1672 wDis = wInter / 1.67E-24 * h['g_'+str(mupid)].Eval(muonEnergy ) * 1E-27
1673 rc = h['wDIS_'+str(mupid)].Fill(wDis)
1674 W = wLHC*wDis
1675 else:
1676 W = norm*muon.GetWeight()
1677 rc = h['inMu_0'+str(mupid)].Fill(muonEnergy,0,W)
1678 rc = h['xyz_mu_0'+str(mupid)].Fill(muon.GetStartX(),muon.GetStartY(),muon.GetStartZ(),W)
1679
1680 motherOf = {}
1681 for d in sTree.MCTrack:
1682 motherOf[d] = d.GetMotherId()
1683 daughterOf = dict(zip(motherOf.values(), motherOf.keys()))
1684
1685 neutrals = {}
1686 tagged = -1
1687 n = -1
1688 for m in sTree.MCTrack:
1689 n+=1
1690 Apid = abs(m.GetPdgCode())
1691 if Apid in neuPartList:
1692# check origin of neutral particle, only count particles from outside
1693 nodeOrigin = geoManager.FindNode(m.GetStartX(),m.GetStartY(),m.GetStartZ())
1694 om = nodeOrigin.GetName()
1695 if not( om=='VTI18_1' or om=='Vrock_1' or om=='cave_1'): continue
1696# find end vertex:
1697 endZ = 999.
1698 concrete = False
1699 if n in daughterOf:
1700 d = daughterOf[n]
1701 endZ = d.GetStartZ()
1702 nodeInter = geoManager.FindNode(d.GetStartX(),d.GetStartY(),d.GetStartZ())
1703 if nodeInter.GetName()=='VTI18_1' or nodeInter.GetName()=='Vrock_1': concrete = True
1704 if not endZ<999.: continue
1705 origin = ''
1706 if m.GetProcID()!=0 : origin = '_secondary'
1707 rc = h['inE_'+str(m.GetPdgCode())+origin].Fill(m.GetEnergy(),endZ,W)
1708 if concrete: rc = h['inE_Concrete_'+str(m.GetPdgCode())+origin].Fill(m.GetEnergy(),endZ,W)
1709 elif debug>2 and m.GetEnergy()>75. and endZ>-25. and endZ<25. and Apid==130 :
1710 tagged = n
1711 print('tagged ',h['f'].GetName(),nEv,tagged)
1712 #elif endZ<-100:
1713 # print('what is this place?',nodeInter.GetName(),m.GetPdgCode(),m.GetP(),d.GetProcName())
1714 # decay or photonpair prod or hadronic inelastic
1715 nm = nodeInter.GetName()
1716 if nm!='VTI18_1' and nm!='Vrock_1' and nm!='cave_1' and endZ<25: # only interested in emulsion
1717 # if Apid==130: print('what is this place called',h['f'].GetName(),nEv,n,d,nm)
1718 neutrals[n] = [m.GetProcID(),nm,d.GetStartX(),d.GetStartY(),d.GetStartZ()]
1719 if sTree.ScifiPoint.GetEntries()<1: continue
1720 if len(neutrals)<1 : continue
1721# find neutral with highest energy
1722 eMax, neu = -1,0
1723 for x in neutrals:
1724 E = sTree.MCTrack[x].GetEnergy()
1725 if E>eMax:
1726 eMax,neu=E,x
1727 if debug>3: print('debug',len(neutrals),E,sTree.MCTrack[x])
1728# look for veto hit
1729# special case, low energy neutrons can also make a hit.
1730 vetoPoint, zMinCharged = -1,999.
1731 for det in [sTree.ScifiPoint,sTree.MuFilterPoint]: # also loop over MuFilterPoint to get hits in veto station
1732 nP = -1
1733 for p in det:
1734 nP+=1
1735 part = PDG.GetParticle(p.PdgCode())
1736 charge = -1
1737 if part: charge = part.Charge()
1738 if charge != 0 and p.GetZ() < zMinCharged :
1739 vetoPoint,zMinCharged = p,p.GetZ()
1740 if debug>2 and tagged >0:
1741 print("What happened",neu)
1742 Npart = sTree.MCTrack[neu]
1743 E = Npart.GetEnergy()
1744 pid = Npart.GetPdgCode()
1745 # neutral particle interaction point, neutral particle origin, muon interaction point
1746 neutInterVx = [neutrals[neu][2],neutrals[neu][3],neutrals[neu][4]]
1747 veto = zMinCharged < neutInterVx[2]
1748 rc = h['xyz_Inter_'+str(pid)].Fill(neutInterVx[0],neutInterVx[1],neutInterVx[2],W)
1749 if E>100.: rc = h['xyzE100_Inter_'+str(pid)].Fill(neutInterVx[0],neutInterVx[1],neutInterVx[2],W)
1750 rc = h['xyz_origin_'+str(pid)].Fill(Npart.GetStartX(),Npart.GetStartY(),Npart.GetStartZ(),W)
1751 rc = h['xyz_muInter_'+str(pid)].Fill(muInterVx.X(),muInterVx.Y(),muInterVx.Z(),W)
1752 if veto:
1753 rc = h['xyzVeto_Inter_'+str(pid)].Fill(neutInterVx[0],neutInterVx[1],neutInterVx[2],W)
1754 rc = h["E_veto_"+str(pid)].Fill(E,neutInterVx[2],W)
1755 p = vetoPoint
1756 P = ROOT.TVector3(p.GetPx(),p.GetPy(),p.GetPz())
1757 pidVeto = p.PdgCode()
1758 if not pidVeto in pidsDictRev:
1759 print("!!!! unknown pid",pidVeto)
1760 L = len(pidsDictRev)
1761 pidsDictRev[pidVeto] = L
1762 rc = h['pidsDict'].SetBinContent(L+1,pidVeto)
1763 rc = h["E_vetoParticle_"+str(pid)].Fill(P.Mag(),pidsDictRev[pidVeto],W)
1764 rc = h["z_veto_"+str(pid)].Fill(neutInterVx[2],zMinCharged)
1765 else:
1766 rc = h["E_"+str(pid)].Fill(E,neutInterVx[2],W)
1767 if Npart.GetProcID()==0:
1768 if veto: rc = h["Eprim_veto_"+str(pid)].Fill(E,neutInterVx[2],W)
1769 else: rc = h["Eprim_"+str(pid)].Fill(E,neutInterVx[2],W)
1770#
1771# check for particles in the first muon station to study about increased veto detector
1772 for p in sTree.MuFilterPoint:
1773 if p.GetZ()<50. and p.GetZ()>0:
1774 track = p.GetTrackID()
1775 if track<0: continue
1776 V = sTree.MCTrack[track]
1777 if V.GetStartZ()<-35:
1778 rc = h["xy_z-30_"+str(pid)].Fill(p.GetX(),p.GetY(),W)
1779 if veto: rc = h["xy_z-30_veto_"+str(pid)].Fill(p.GetX(),p.GetY(),W)
1780 h['f'].Close()
1781 h['3d']=h['xyz_muInter_130'].Clone('3d')
1782 h['3d'].Add(h['xyz_muInter_2112'])
1783 h['3d'].Add(h['xyz_muInter_-2112'])
1784 if pythia6:
1785 if options.Emin==0: ut.writeHists(h,"muonDISfull-"+str(sMin)+"_"+str(sMax)+"_"+str(cycle)+".root")
1786 else: ut.writeHists(h,"muonDISfull-"+str(options.Emin)+'_'+str(sMin)+"_"+str(sMax)+"_"+str(cycle)+".root")
1787 else:
1788 ut.writeHists(h,"muonGeant4full-"+str(sMin)+"_"+str(sMax)+".root")

◆ muonDISProduction()

muonDis.muonDISProduction (   cycle,
  ecut = 1.,
  strippedEvents = False 
)

Definition at line 2244 of file muonDis.py.

2244def muonDISProduction(cycle,ecut=1.,strippedEvents=False):
2245 for run in range(1,11):
2246 for k in [0,1000]:
2247 prod = str(run+cycle*100+k)
2248 if strippedEvents: fn = "muonDis_"+prod+"_stripped.root"
2249 else: fn = "muonDis_"+prod+".root"
2250 command = "python $SNDSW_ROOT/macro/run_simScript.py --shiplhc"
2251 os.system(command+" -f "+fn+" --MuDIS --output run_"+prod+" -n 99999999 --eMin "+str(ecut) + " &")
2252 while count_python_processes('run_simScript')>ncpus: time.sleep(200)
2253

◆ muonPreTransport()

muonDis.muonPreTransport ( )

Definition at line 182 of file muonDis.py.

182def muonPreTransport():
183 f = ROOT.TFile(options.muonIn)
184 nt = f.Get("nt")
185 foutName = options.muonIn.replace('.root','_z'+str(options.z)+'.root')
186 fout = ROOT.TFile(foutName,'recreate')
187 variables = ""
188 for n in range(nt.GetListOfLeaves().GetEntries()):
189 variables+=nt.GetListOfLeaves()[n].GetName()
190 if n < nt.GetListOfLeaves().GetEntries()-1: variables+=":"
191 sTree = ROOT.TNtupleD("nt","muon",variables)
192 for n in range(nt.GetEntries()):
193 rc = nt.GetEvent(n)
194 E = nt.E - 27.
195 if E>0:
196 column = []
197 lam = ((options.z+SND_Z)-nt.z)/nt.pz
198 for n in range(nt.GetListOfLeaves().GetEntries()):
199 val = nt.GetListOfLeaves()[n].GetValue()
200 if nt.GetListOfLeaves()[n].GetName()=='E': val = E
201 if nt.GetListOfLeaves()[n].GetName()=='x': val = nt.x+lam*nt.px
202 if nt.GetListOfLeaves()[n].GetName()=='y': val = nt.y+lam*nt.py
203 if nt.GetListOfLeaves()[n].GetName()=='z': val = options.z+SND_Z
204 column.append(val)
205 theTuple = array('d',column)
206 rc = sTree.Fill(theTuple)
207 sTree.AutoSave()
208 fout.Close()
209

◆ muonRateAtSND()

muonDis.muonRateAtSND (   withFaser = False,
  withEff = False,
  version = 1 
)

Definition at line 812 of file muonDis.py.

812def muonRateAtSND(withFaser=False,withEff=False,version=1):
813 norm = normalization[version]
814 if withEff:
815 ut.readHists(h,"efficiency.root")
816 h['efficiency'] = h['eff'].FindObject('eloss_px').Clone('efficiency')
817 fname = 'pol4'
818 rc = h['efficiency'].Fit(fname,'S','',20,300.)
819 effFun = h['efficiency'].GetFunction(fname)
820 effPMax = 200.
821 ut.readHists(h,"meanEloss.root")
822 eLoss = h['TCeloss'].FindObject('pol3').Clone('eLoss')
823 Rmult = 100 # if 1, MS scattering replaced by efficiency
824 fnames = {13:"unit30_Nm.root",-13:"unit30_Pm.root"}
825 ut.bookCanvas(h,'muDIS_SNDXY_1','muons',1200,900,1,1)
826 ut.bookCanvas(h,'muDIS_SNDXY_2','muons',1200,900,1,1)
827 ut.bookCanvas(h,'muDIS_SNDE','muons',1200,900,1,1)
828 for mu in fnames:
829 fin = ROOT.TFile(fnames[mu]) # id:px:py:pz:x:y:z:w
830 ut.bookHist(h,'xy_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
831 ut.bookHist(h,'xyMS_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ',200,-100.,100.,200,-100.,100.)
832 ut.bookHist(h,'xyW_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
833 ut.bookHist(h,'xysco_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ' , 200,-100.,100.,200,-100.,100.)
834 ut.bookHist(h,'xyMSW_'+str(mu), 'x/y;x [cm];y [cm];N/sec/cm^{2} ',200,-100.,100.,200,-100.,100.)
835 ut.bookHist(h,'Esco_'+str(mu), 'E in SND',1000,0.,5000.)
836 ut.bookHist(h,'Esnd_'+str(mu), 'E in SND',1000,0.,5000.)
837 ut.bookHist(h,'EsndMS_'+str(mu), 'E in SND',1000,0.,5000.)
838 ut.bookHist(h,'Efaser_'+str(mu), 'E in faser',1000,0.,5000.)
839 ut.bookHist(h,'EfaserMS_'+str(mu), 'E in faser',1000,0.,5000.)
840 ut.bookHist(h,'tanThetaXY_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01)
841 ut.bookHist(h,'tanThetaXY30_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01)
842 ut.bookHist(h,'tanThetaXYSND_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01)
843 ut.bookHist(h,'tanThetaXYMS_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01)
844 ut.bookHist(h,'tanThetaXYMSfaser_'+str(mu), 'tan theta X/Y',200,-0.01,0.01,200,-0.01,0.01)
845 ut.bookHist(h,'theta', 'mult scattering angle',200,-1.,1.)
846 for sTree in fin.Get("nt"):
847 px,py,pz = sTree.px,sTree.py,sTree.pz
848 m2cm = 1.
849 if version == 0: m2cm = 100
850 x,y,z = sTree.x,sTree.y,sTree.z*m2cm # how could have this happened? used u.m for conversion!
851 # if abs(x)<50. and abs(y)<50. : continue
852 pid,w = int(sTree.id),sTree.w
853 p = ROOT.TMath.Sqrt(px*px+py*py+pz*pz)
854 wLHC = norm*w
855 rc = h["Esco_"+str(mu)].Fill(p,wLHC)
856 rc = h['xysco_'+str(pid)].Fill(x,y,wLHC)
857# apply some multiple scattering and energy loss, and adjust weight for efficiency replacing MS angle
858 if Rmult==-1:
859 if p < effPMax: wLHC=wLHC*effFun.Eval(p)
860 Psnd = p - eLoss.Eval(p)
861 else:
862 Psnd = p - 30.
863 if version == 1: Psnd = p - 26.
864 lam = ( SND_Z-z)/pz
865 mxex = x+lam*px
866 myex = y+lam*py
867 # print('debug',SND_Z,z,lam,x,mxex,y,myex,px,py)
868 rc = h['xy_'+str(pid)].Fill(mxex,myex)
869 rc = h['xyW_'+str(pid)].Fill(mxex,myex,wLHC)
870 tanThetaX = ROOT.TMath.Tan(ROOT.TMath.ATan2(px,pz))
871 tanThetaY = ROOT.TMath.Tan(ROOT.TMath.ATan2(py,pz))
872 rc = h['tanThetaXY_'+str(mu)].Fill(tanThetaX,tanThetaY,wLHC)
873 if Psnd>0:
874 rc = h['tanThetaXY30_'+str(mu)].Fill(tanThetaX,tanThetaY,wLHC)
875 if -8 > mxex and mxex >-49 and 15 < myex and myex < 57:
876 rc = h["Esnd_"+str(mu)].Fill(Psnd,wLHC)
877 if -12 < mxex and mxex < 12 and -12 < myex and myex < 12:
878 rc = h["Efaser_"+str(mu)].Fill(Psnd,wLHC)
879# 13.6/P sqrt(L/x0)(1+0.038 log(L/x0)), L/X0 = 700
880 for r in range(Rmult):
881 mxexMS = mxex
882 myexMS = myex
883 if Rmult>1:
884 thetaX = rnr.Gaus(0.,h['MS'].Eval(p))
885 thetaY = rnr.Gaus(0.,h['MS'].Eval(p))
886 rc = h['theta'].Fill(thetaX)
887 rc = h['theta'].Fill(thetaY)
888 mxexMS = mxex + (SND_Z-z)*ROOT.TMath.Tan(thetaX)
889 myexMS = myex + (SND_Z-z)*ROOT.TMath.Tan(thetaY)
890 tanThetaX = ROOT.TMath.Tan(ROOT.TMath.ATan2(px,pz)+thetaX)
891 tanThetaY = ROOT.TMath.Tan(ROOT.TMath.ATan2(py,pz)+thetaY)
892 rc = h['tanThetaXYSND_'+str(mu)].Fill(tanThetaX,tanThetaY,wLHC/Rmult)
893 if Psnd >0:
894 rc = h['xyMS_'+str(pid)].Fill(mxexMS,myexMS)
895 rc = h['xyMSW_'+str(pid)].Fill(mxexMS,myexMS,wLHC/Rmult)
896 if -8 > mxexMS and mxexMS >-49 and 15 < myexMS and myexMS < 57:
897 rc = h["EsndMS_"+str(mu)].Fill(Psnd,wLHC/Rmult)
898 if Psnd >0: rc = h['tanThetaXYMS_'+str(mu)].Fill(tanThetaX,tanThetaY,wLHC/Rmult)
899 if -12 < mxexMS and mxexMS < 12 and -12 < myexMS and myexMS < 12:
900 rc = h["EfaserMS_"+str(mu)].Fill(Psnd,wLHC/Rmult)
901 if Psnd >0: rc = h['tanThetaXYMSfaser_'+str(mu)].Fill(tanThetaX,tanThetaY,wLHC/Rmult)
902 for hname in ["Esco_"+str(mu),"Esnd_"+str(mu),"EsndMS_"+str(mu),"Efaser_"+str(mu),"EfaserMS_"+str(mu)]:
903 if 'I-'+hname in h: A=h.pop('I-'+hname)
904 ROOT.gROOT.cd()
905 ut.makeIntegralDistrib(h,hname)
906 h['I-'+hname].GetXaxis().SetRangeUser(0.,10000.)
907 h['I-'+hname].SetMinimum(1E-6)
908 pname = PDG.GetParticle(mu).GetName()
909 if hname.find("Esco")==0:
910 print("at scoring plane %s, %s: %5.2F N/sec "%(pname,hname,h['I-'+hname].GetBinContent(1)))
911 else:
912 print("at SND dE=-30GeV %s, %s: %5.2F N/sec "%(pname,hname,h['I-'+hname].GetBinContent(1)))
913 h['snd_1'] = ROOT.TLine(-49,15,-8,15)
914 h['snd_2'] = ROOT.TLine(-49,57,-8,57)
915 h['snd_3'] = ROOT.TLine(-49,57,-49,15)
916 h['snd_4'] = ROOT.TLine(-8,57,-8,15)
917 h['fas_1'] = ROOT.TLine(-12,12,12,12)
918 h['fas_2'] = ROOT.TLine(-12,-12,12,-12)
919 h['fas_3'] = ROOT.TLine(-12,-12,-12,12)
920 h['fas_4'] = ROOT.TLine(12,-12,12,12)
921 for n in range(1,5):
922 h['snd_'+str(n)] .SetLineColor(ROOT.kBlack)
923 h['snd_'+str(n)] .SetLineWidth(5)
924 h['snd_'+str(n)] .SetLineStyle(1)
925 h['fas_'+str(n)] .SetLineColor(ROOT.kCyan)
926 h['fas_'+str(n)] .SetLineWidth(5)
927 h['fas_'+str(n)] .SetLineStyle(2)
928 for Y in ['xyMSW_','xysco_']:
929 for mu in fnames:
930 h[Y+str(mu)].SetStats(0)
931 if version==0:
932 h[Y+str(mu)].GetXaxis().SetRangeUser(-60.,60.)
933 h[Y+str(mu)].GetYaxis().SetRangeUser(-60.,60.)
934 else:
935 # h[Y+str(mu)].SetMaximum(10)
936 h[Y+str(mu)].GetXaxis().SetRangeUser(-80.,80.)
937 h[Y+str(mu)].GetYaxis().SetRangeUser(-80.,80.)
938 if mu<0: h[Y+str(mu)].SetTitle("#mu^{+}")
939 if mu>0: h[Y+str(mu)].SetTitle("#mu^{-}")
940 tc = h['muDIS_SNDXY_1'].cd()
941 tc.SetRightMargin(0.2)
942 h[Y+'13'].GetZaxis().SetTitleOffset(1.3)
943 h[Y+'13'].Draw('colz')
944 tc.Update()
945 pal = h[Y+'13'].GetListOfFunctions()[0]
946 pal.SetX1NDC(0.81)
947 pal.SetX2NDC(0.85)
948 if Y=='xyMSW_':
949 for n in range(1,5):
950 h['snd_'+str(n)].Draw('same')
951 if withFaser: h['fas_'+str(n)].Draw('same')
952 tc = h['muDIS_SNDXY_2'].cd()
953 tc.SetRightMargin(0.2)
954 h[Y+'-13'].GetZaxis().SetTitleOffset(1.5)
955 h[Y+'-13'].Draw('colz')
956 tc.Update()
957 pal = h[Y+'-13'].GetListOfFunctions()[0]
958 pal.SetX1NDC(0.8)
959 pal.SetX2NDC(0.85)
960 if Y=='xyMSW_':
961 for n in range(1,5):
962 h['snd_'+str(n)].Draw('same')
963 if withFaser: h['fas_'+str(n)].Draw('same')
964 if Y=='xyMSW_':
965 myPrint(h['muDIS_SNDXY_1'],'muDIS_SNDXY_muM')
966 myPrint(h['muDIS_SNDXY_2'],'muDIS_SNDXY_muP')
967 else:
968 myPrint(h['muDIS_SNDXY_1'],'muDIS_SCOXY_muM')
969 myPrint(h['muDIS_SNDXY_2'],'muDIS_SCOXY_muP')
970
971 R = h['I-EsndMS_13'].GetBinContent(1)+h['I-EsndMS_-13'].GetBinContent(1)
972 print('mu+ and mu- at SND: %5.2F N/sec %5.2F Hz/cm2'%(R,R/(41.*42.)))
973 R = h['I-EfaserMS_13'].GetBinContent(1)+h['I-EfaserMS_-13'].GetBinContent(1)
974 print('mu+ and mu- at Faser: %5.2F N/sec %5.2F Hz/cm2'%(R,R/(24.*24.)))
975 rc = h['muDIS_SNDE'].cd(1)
976 rc.SetLogy()
977 h['I-Esnd_-13'].SetLineColor(ROOT.kRed)
978 h['I-Esnd_13'].SetLineColor(ROOT.kBlue)
979 h['I-Esnd_13'].SetFillStyle(2)
980 h['I-Esnd_13'].SetStats(0)
981 h['I-Esnd_-13'].SetStats(0)
982 h['I-Esnd_-13'].SetFillStyle(2)
983 h['I-Esnd_13'].SetMinimum(1E-3)
984 h['I-Esnd_13'].SetTitle(';>P [GeV/c]; Hz')
985 h['I-Esnd_13'].Draw('BAR')
986 h['I-Esnd_-13'].Draw('sameBAR')
987 h['t-13' ]=ROOT.TLatex(500,0.1,PDG.GetParticle(-13).GetName())
988 h['t13' ]=ROOT.TLatex(2500,0.3,PDG.GetParticle(13).GetName())
989 h['t13' ].SetTextColor(ROOT.kWhite)
990 h['t13' ].Draw('same')
991 h['t-13'].Draw('same')
992 h['muDIS_SNDE'].Print('muDIS_ESND.png')
993

◆ myPrint()

muonDis.myPrint (   tc,
  tname,
  pathToPlots = "/mnt/hgfs/microDisk/CERNBOX/SND@LHC/MuonDis/" 
)

Definition at line 112 of file muonDis.py.

112def myPrint(tc,tname,pathToPlots = "/mnt/hgfs/microDisk/CERNBOX/SND@LHC/MuonDis/"):
113 for z in ['.png','.pdf','.root']:
114 tc.Print(tname+z)
115 os.system('mv '+tname+z+' '+pathToPlots)
116

◆ nucl_cross()

muonDis.nucl_cross (   Ebeam,
  eps,
  A 
)

Definition at line 47 of file muonDis.py.

47def nucl_cross(Ebeam,eps,A):
48 # returns xsec in microbarns
49 # A in g / mol ??
50 Mmu = muonMass
51 m1sq = 0.54
52 m2sq = 1.8
53 alpha = 1./137.03599976
54 nu = eps/Ebeam
55 t = Mmu**2*nu**2/(1.-nu)
56 k = 1.-2./nu+2./nu**2
57 x = 0.00282*ROOT.TMath.Power(A,1./3.)*photoabsorb(eps)
58 if eps<5: sigma = 0
59 # if eps<3: sigma = 0 # 3.8.2015, test, text (PDG) says, <5 GeV not reliable
60 else: sigma = alpha/(2.*ROOT.TMath.Pi())*A*photoabsorb(eps)*nu*(0.75*G(x)*(k*ROOT.TMath.Log(1+m1sq/t) - \
61 k*m1sq/(m1sq+t) - 2.*Mmu**2/t)+0.25*(k*ROOT.TMath.Log(1+m2sq/t)- 2.*Mmu**2/t) + \
62 Mmu**2/(2.*t)*(0.75*G(x)*m1sq/(m1sq+t)+0.25*m2sq/t*ROOT.TMath.Log(1.+t/m2sq)))
63 return sigma
64

◆ photoabsorb()

muonDis.photoabsorb (   eps)

Definition at line 43 of file muonDis.py.

43def photoabsorb(eps):
44 return 114.3 + 1.647*(ROOT.TMath.Log(0.0213*eps))**2

◆ plotMuDisCrossSection()

muonDis.plotMuDisCrossSection ( )

Definition at line 1536 of file muonDis.py.

1536def plotMuDisCrossSection():
1537 fin = ROOT.TFile('muDIScrossSec.root')
1538 for pid in ['13','-13']:
1539 h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid)
1540
1541 SigmaAnalyticVsEnergy(A=1)
1542 SigmaAnalyticVsA(Ebeam=500)
1543
1544 ut.bookCanvas(h,'sec','xsec',900,600,1,1)
1545 ut.bookHist(h,'muDISXsec',';E [GeV];#sigma [mb]',100,0.,10000.)
1546 ut.bookHist(h,'muDISXsecA',';A ; #sigma [mb]',100,0.,100.)
1547 h['sec'].cd(1)
1548 h['muDISXsec'].SetMinimum(0.)
1549 h['muDISXsec'].SetMaximum(30.E-3)
1550 h['muDISXsecA'].SetMinimum(0.)
1551 h['muDISXsecA'].SetMaximum(3000.E-3)
1552 h['muDISXsec'].Draw()
1553 h['muDISXsec'].SetStats(0)
1554 h['g_13'].SetLineColor(ROOT.kGreen)
1555 h['g_-13'].SetLineColor(ROOT.kRed)
1556 h['g_13'].SetLineWidth(4)
1557 h['g_-13'].SetLineWidth(4)
1558 h['g_13'].Draw('same')
1559 h['g_-13'].Draw('same')
1560 h['AnalyticCross'].SetLineWidth(4)
1561 h['AnalyticCross'].SetLineColor(ROOT.kMagenta)
1562 h['AnalyticCross'].Draw('same')
1563 h['lxsec']=ROOT.TLegend(0.14,0.71,0.64,0.84)
1564 h['lxsec'].AddEntry(h['AnalyticCross'],'analytic cross section (Bezrukov and Bugaev, A=1)','PL')
1565 h['lxsec'].AddEntry(h['g_-13'],'Pythia6 for #mu^{+} p','PL')
1566 h['lxsec'].AddEntry(h['g_13'],'Pythia6 for #mu^{-} p','PL')
1567 h['lxsec'].Draw('same')
1568
1569 myPrint(h['sec'],'muonXsec')
1570

◆ readXsec()

muonDis.readXsec (   p)

Definition at line 434 of file muonDis.py.

434def readXsec(p):
435 f = open("fort.11")
436 tmp = None
437 X = f.readlines()
438 for i in range(1,len(X)):
439 l = X[len(X)-i]
440 if l.find('All included')<0: continue
441 tmp = l.replace("\n",'').split('I')
442 break
443 if not tmp:
444 print("empty file")
445 return 0
446 return float(tmp[3].replace('D','E'))
447

◆ rotate()

muonDis.rotate (   ctheta,
  stheta,
  cphi,
  sphi,
  px,
  py,
  pz 
)

Definition at line 217 of file muonDis.py.

217def rotate(ctheta,stheta,cphi,sphi,px,py,pz):
218 #rotate around y-axis
219 px1=ctheta*px+stheta*pz
220 pzr=-stheta*px+ctheta*pz
221 #rotate around z-axis
222 pxr=cphi*px1-sphi*py
223 pyr=sphi*px1+cphi*py
224 return pxr,pyr,pzr
225

◆ selectEvents()

muonDis.selectEvents (   Ecut = 10)

Definition at line 2274 of file muonDis.py.

2274def selectEvents(Ecut=10):
2275 for cycle in range(10):
2276 for run in range(1,11):
2277 for k in [0,1000]:
2278 prod = str(run+cycle*100+k)
2279 f = ROOT.TFile('muonDis_'+str(prod)+'.root')
2280 sTree = f.DIS
2281 fstripped = ROOT.TFile('muonDis_'+str(prod)+'_stripped.root','recreate')
2282 newTree = sTree.CloneTree(0)
2283 for sTree in f.DIS:
2284 flag = False
2285 for P in sTree.Particles:
2286 if P.GetPdgCode() in [130,2112,-2112]:
2287 if P.Energy() > Ecut: rc = newTree.Fill()
2288 newTree.AutoSave()
2289 f.Close()
2290 print('f closed',prod)
2291 fstripped.Close()

◆ SigmaAnalytic()

muonDis.SigmaAnalytic (   Ebeam = 5000.,
  A = 1,
  nsteps = 10000 
)

Definition at line 65 of file muonDis.py.

65def SigmaAnalytic(Ebeam=5000.,A=1,nsteps = 10000):
66 gname = 'xsec_'+str(A)
67 h[gname] = ROOT.TGraph(nsteps)
68 deps = Ebeam/float(nsteps)
69 eps = deps/2.
70 for i in range(nsteps):
71 sigma = nucl_cross(Ebeam,eps,A)
72 h[gname].SetPoint(i,eps/Ebeam,sigma)
73 eps+=deps
74

◆ SigmaAnalyticVsA()

muonDis.SigmaAnalyticVsA (   Ebeam = 500)

Definition at line 86 of file muonDis.py.

86def SigmaAnalyticVsA(Ebeam=500):
87 nsteps = 100
88 gname = 'AnalyticCross_'+str(Ebeam)
89 h[gname] = ROOT.TGraph(nsteps)
90 h[gname].SetName(gname)
91 i=0
92 for A in range(1,nsteps+1):
93 SigmaAnalytic(Ebeam,A)
94 h[gname].SetPoint(i,A, h['xsec_'+str(A)].Integral()/1000.) # convert to mbarns
95 i+=1

◆ SigmaAnalyticVsEnergy()

muonDis.SigmaAnalyticVsEnergy (   A = 1)

Definition at line 75 of file muonDis.py.

75def SigmaAnalyticVsEnergy(A=1):
76 nsteps = len(EnergyScan)
77 gname = 'AnalyticCross_'+str(A)
78 h[gname] = ROOT.TGraph(nsteps)
79 h[gname].SetName(gname)
80 i=0
81 for Ebeam in EnergyScan:
82 SigmaAnalytic(Ebeam,A)
83 h[gname].SetPoint(i,Ebeam, h['xsec_'+str(A)].Integral()/A/1000.) # convert to mbarns
84 i+=1
85

◆ signalNeutrinos()

muonDis.signalNeutrinos ( )

Definition at line 2292 of file muonDis.py.

2292def signalNeutrinos():
2293 nudict = {'e':'$\\nu_e$','ae':'$\overline{\\nu}_e$','mu':'$\\nu_\mu$','amu':'$\overline{\\nu}_\mu$','tau':'$\\nu_\\tau$','atau':'$\overline{\\nu}_\\tau$'}
2294 ut.readHists(h,'/mnt/hgfs/microDisk/CERNBOX/SND@LHC/Neutrinos/SND_Neutrinos_Interacting_CCDIS.root')
2295 for p in ['','a']:
2296 for nu in [ 'e','mu','tau']:
2297 line = ''
2298 hname = 'h'+p+'nu'+nu
2299 ut.makeIntegralDistrib(h,hname)
2300 hist = h['I-'+hname]
2301 line += nudict[p+nu]
2302 for E in [10,100,200,300,500,1000]:
2303 n = hist.FindBin(E)
2304 N = hist.GetBinContent(n)
2305 line += ' & '+'$%5.1F$'%(N)
2306 line +=' \\\\'
2307 print(line)
2308

◆ thermNeutron()

muonDis.thermNeutron ( )

Definition at line 1789 of file muonDis.py.

1789def thermNeutron():
1790# under development, not used, low energy neutrons not stored.
1791 f=ROOT.TFile.Open('root://eospublic.cern.ch//eos/experiment/ship/user/truf/SND/muonDis/run_1_0/ship.conical.muonDIS-TGeant4.root')
1792 ut.bookHist(h,'n_2112','logE MeV',100,-10,2)
1793 ut.bookHist(h,'n_-2112','logE MeV',100,-10,2)
1794 L=ROOT.TLorentzVector()
1795 nMult = 10
1796 norm = 5.83388
1797 fin = ROOT.TFile('muDIScrossSec.root')
1798 ROOT.gROOT.cd()
1799 for pid in ['13','-13']:
1800 h['g_'+pid] = fin.Get('g_'+pid).Clone('g_'+pid)
1801 for sTree in f.Get("cbmsim"):
1802 muon = sTree.MCTrack[0]
1803 muonEnergy = muon.GetEnergy()
1804 mupid = muon.GetPdgCode()
1805 wLHC = norm*muon.GetWeight()/float(nMult)
1806 wInter = sTree.MCTrack[2].GetWeight()
1807 wDis = wInter / 1.67E-24 * h['g_'+str(mupid)].Eval(muonEnergy ) * 1E-27
1808 W = wLHC*wDis
1809 for m in sTree.MCTrack:
1810 if abs(m.GetPdgCode()) == 2112:
1811 m.Get4Momentum(L)
1812 logE = ROOT.TMath.Log10( 1000.*(L.Energy()-L.M()))
1813 rc = h['n_'+str(m.GetPdgCode())].Fill(logE,W)

Variable Documentation

◆ cycle

muonDis.cycle

Definition at line 2322 of file muonDis.py.

◆ default

muonDis.default

Definition at line 97 of file muonDis.py.

◆ dest

muonDis.dest

Definition at line 97 of file muonDis.py.

◆ Eloss

int muonDis.Eloss = 30.

Definition at line 32 of file muonDis.py.

◆ EnergyScan

list muonDis.EnergyScan = [10,20,30,40,50,75,100,200,300,400,500,1000,2000,3000,4000,5000,6000,7000,8000,10000]

Definition at line 41 of file muonDis.py.

◆ FASERNU

dict muonDis.FASERNU = {}

Definition at line 2218 of file muonDis.py.

◆ float

muonDis.float

Definition at line 101 of file muonDis.py.

◆ flukaDict

dict muonDis.flukaDict = {10:-13,11:13}

Definition at line 13 of file muonDis.py.

◆ function

str muonDis.function = "13.6/1000.*sqrt([0])/x*(1.+0.038*log([0]))"

Definition at line 34 of file muonDis.py.

◆ h

dict muonDis.h = {}

Definition at line 29 of file muonDis.py.

◆ help

muonDis.help

Definition at line 97 of file muonDis.py.

◆ int

muonDis.int

Definition at line 97 of file muonDis.py.

◆ masssq

dict muonDis.masssq = {}

Definition at line 210 of file muonDis.py.

◆ muonIn

muonDis.muonIn

Definition at line 2324 of file muonDis.py.

◆ muonMass

float muonDis.muonMass = 0.105658

Definition at line 21 of file muonDis.py.

◆ ncpus

muonDis.ncpus = multiprocessing.cpu_count()

Definition at line 2232 of file muonDis.py.

◆ nEvents

muonDis.nEvents

Definition at line 2324 of file muonDis.py.

◆ nMult

muonDis.nMult

Definition at line 2322 of file muonDis.py.

◆ normalization

dict muonDis.normalization = {}

Definition at line 23 of file muonDis.py.

◆ nStart

muonDis.nStart

Definition at line 2322 of file muonDis.py.

◆ nucleon

muonDis.nucleon

Definition at line 2327 of file muonDis.py.

◆ options

muonDis.options = parser.parse_args()

Definition at line 110 of file muonDis.py.

◆ parser

muonDis.parser = ArgumentParser()

Definition at line 96 of file muonDis.py.

◆ path

muonDis.path

Definition at line 2324 of file muonDis.py.

◆ PDG

muonDis.PDG = ROOT.TDatabasePDG.Instance()

Definition at line 10 of file muonDis.py.

◆ pids

list muonDis.pids = [-3222, -2212, -321, -211,-13, -11, 11,13, 211, 321, 2212, 3112]

Definition at line 15 of file muonDis.py.

◆ pidsDict

dict muonDis.pidsDict = {}

Definition at line 16 of file muonDis.py.

◆ pidsDictRev

dict muonDis.pidsDictRev = {}

Definition at line 17 of file muonDis.py.

◆ pythia6

muonDis.pythia6

Definition at line 2324 of file muonDis.py.

◆ rMax

muonDis.rMax

Definition at line 2324 of file muonDis.py.

◆ rMin

muonDis.rMin

Definition at line 2324 of file muonDis.py.

◆ rnr

muonDis.rnr = ROOT.TRandom()

Definition at line 11 of file muonDis.py.

◆ scaleFactor

tuple muonDis.scaleFactor = (42.*42.)/(24.*24.)

Definition at line 2228 of file muonDis.py.

◆ sMax

muonDis.sMax

Definition at line 2322 of file muonDis.py.

◆ sMin

muonDis.sMin

Definition at line 2322 of file muonDis.py.

◆ SND_Z

float muonDis.SND_Z = 483.262*u.m

Definition at line 31 of file muonDis.py.

◆ str

muonDis.str

Definition at line 105 of file muonDis.py.

◆ theSeed

int muonDis.theSeed = 0

Definition at line 8 of file muonDis.py.

◆ type

muonDis.type

Definition at line 97 of file muonDis.py.

◆ xOverX0

int muonDis.xOverX0 = 75*u.m/10.02*u.cm

Definition at line 38 of file muonDis.py.