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 641 of file muonDis.py.

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

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

◆ boundaries()

muonDis.boundaries ( )

Definition at line 984 of file muonDis.py.

984def boundaries():
985 h['snd_1'] = ROOT.TLine(-49,15,-8,15)
986 h['snd_2'] = ROOT.TLine(-49,57,-8,57)
987 h['snd_3'] = ROOT.TLine(-49,57,-49,15)
988 h['snd_4'] = ROOT.TLine(-8,57,-8,15)
989 h['fas_1'] = ROOT.TLine(-12,12,12,12)
990 h['fas_2'] = ROOT.TLine(-12,-12,12,-12)
991 h['fas_3'] = ROOT.TLine(-12,-12,-12,12)
992 h['fas_4'] = ROOT.TLine(12,-12,12,12)
993 for n in range(1,5):
994 h['snd_'+str(n)] .SetLineColor(ROOT.kBlack)
995 h['snd_'+str(n)] .SetLineWidth(5)
996 h['snd_'+str(n)] .SetLineStyle(1)
997 h['fas_'+str(n)] .SetLineColor(ROOT.kCyan)
998 h['fas_'+str(n)] .SetLineWidth(5)
999 h['fas_'+str(n)] .SetLineStyle(2)
1000

◆ checkProdofMuDIS()

muonDis.checkProdofMuDIS ( )

Definition at line 529 of file muonDis.py.

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

◆ compMuDIS_P6withP8()

muonDis.compMuDIS_P6withP8 ( )

Definition at line 550 of file muonDis.py.

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

◆ 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 2223 of file muonDis.py.

2223def count_python_processes(macroName):
2224 username = pwd.getpwuid(os.getuid()).pw_name
2225 callstring = "ps -f -u " + username
2226# only works if screen is wide enough to print full name!
2227 status = subprocess.check_output(callstring,shell=True)
2228 n=0
2229 for x in status.decode().split('\n'):
2230 if not x.find(macroName)<0 and not x.find('python') <0: n+=1
2231 return n
2232

◆ drawMuon3D()

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

Definition at line 1493 of file muonDis.py.

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

◆ energyOfDISevents()

muonDis.energyOfDISevents ( )

Definition at line 2243 of file muonDis.py.

2243def energyOfDISevents():
2244 for cycle in range(4):
2245 for run in range(1,11):
2246 for k in [0,1000]:
2247 prod = str(run+cycle*100+k)
2248 print('opening','muonDis_'+str(prod)+'.root')
2249 f = ROOT.TFile('muonDis_'+str(prod)+'.root')
2250 if not f: continue
2251 if not f.FindKey('DIS'): continue
2252 sTree = f.DIS
2253 for sTree in f.DIS:
2254 for P in sTree.Particles:
2255 hname = 'E_'+str(P.GetPdgCode())
2256 if not hname in h:
2257 ut.bookHist(h,hname,'Energy '+PDG.GetParticle(P.GetPdgCode()).GetName(),1000,0.,5000.)
2258 rc = h[hname].Fill(P.Energy())
2259 ut.writeHists(h,'energyOfDISevents.root')
2260# maybe better to make a first selection of events with n or KL above 100GeV
2261# cautious, n or KL can also be produced by other particles. Specially problematic,
2262# 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 1001 of file muonDis.py.

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

◆ 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# incoming muon, id:px:py:pz:x:y:z:w
465# outgoing particles, id:px:py:pz
466 fout = ROOT.TFile('muonDis_'+str(options.run)+'.root','recreate')
467 dTree = ROOT.TTree('DIS','muon DIS')
468 iMuon = ROOT.TClonesArray("TParticle")
469 iMuonBranch = dTree.Branch("InMuon",iMuon,32000,-1)
470 dPart = ROOT.TClonesArray("TParticle")
471 dPartBranch = dTree.Branch("Particles",dPart,32000,-1)
472# read file with muons hitting concrete wall
473 fin = ROOT.TFile(options.muonIn) # id:px:py:pz:x:y:z:w
474 sTree = fin.nt
475 nTOT = sTree.GetEntries()
476 nEnd = min(nTOT,options.nStart + options.nEvents)
477# stop pythia printout during loop
478 myPythia.SetMSTU(11, 11)
479 print("start production ",options.nStart,nEnd)
480 nMade = 0
481 for k in range(options.nStart,nEnd):
482 rc = sTree.GetEvent(k)
483# make n events / muon
484 px,py,pz = sTree.px,sTree.py,sTree.pz
485 x,y,z = sTree.x,sTree.y,sTree.z-SND_Z
486 pid,w = int(sTree.id),sTree.w
487 if withElossFunction:
488 E = sTree.E - eLoss.Eval(sTree.E) # energy loss in 75m of rock
489 else:
490 E = sTree.E -27.
491 p = ROOT.TMath.Sqrt(px*px+py*py+pz*pz)
492 Peloss = ROOT.TMath.Sqrt(E*E-muonMass*muonMass)
493 if E < 5. : continue
494 t = sTree.t
495 # px=p*sin(theta)cos(phi),py=p*sin(theta)sin(phi),pz=p*cos(theta)
496 theta = ROOT.TMath.ACos(pz/p)
497 phi = ROOT.TMath.ATan2(py,px)
498 ctheta,stheta = ROOT.TMath.Cos(theta),ROOT.TMath.Sin(theta)
499 cphi,sphi = ROOT.TMath.Cos(phi),ROOT.TMath.Sin(phi)
500 escale = Peloss/p
501 muPart = ROOT.TParticle(pid,0,0,0,0,0,px*escale,py*escale,pz*escale,E,x,y,z,t)
502 muPart.SetWeight(w)
503 myPythia.Initialize('FIXT',mutype[pid],nucleon,E)
504 for n in range(options.nMult):
505 dPart.Clear()
506 iMuon.Clear()
507 iMuon[0] = muPart
508 myPythia.GenerateEvent()
509# remove all unnecessary stuff
510 myPythia.Pyedit(2)
511 for itrk in range(1,myPythia.GetN()+1):
512 did = myPythia.GetK(itrk,2)
513 dpx,dpy,dpz = rotate(ctheta,stheta,cphi,sphi,myPythia.GetP(itrk,1),myPythia.GetP(itrk,2),myPythia.GetP(itrk,3))
514 psq = dpx**2+dpy**2+dpz**2
515 E = ROOT.TMath.Sqrt(getMasssq(did)+psq)
516 part = ROOT.TParticle(did,0,0,0,0,0,dpx,dpy,dpz,E,x,y,z,t)
517 part.SetWeight(w/float(options.nMult))
518# copy to branch
519 nPart = dPart.GetEntries()
520 if dPart.GetSize() == nPart: dPart.Expand(nPart+10)
521 dPart[nPart] = part
522 nMade+=1
523 if nMade%options.heartbeat==0: print('made so far ',options.run,' :',nMade,time.ctime())
524 dTree.Fill()
525 fout.cd()
526 dTree.Write()
527 myPythia.SetMSTU(11, 6)
528 print("finished ",options.run)

◆ missing3Dproj()

muonDis.missing3Dproj (   hist,
  ymin,
  ymax 
)

Definition at line 2298 of file muonDis.py.

2298def missing3Dproj(hist,ymin,ymax):
2299 h[hist+'XZ']=h[hist].Project3D('xz')
2300 h[hist+'XZ'].SetName(hist+'XZ')
2301 h[hist+'XZ'].Reset()
2302 for ix in range(1,h[hist].GetNbinsX()+1):
2303 for iz in range(1,h[hist].GetNbinsZ()+1):
2304 N = 0
2305 for iy in range(ymin,ymax): N+=h[hist].GetBinContent(ix,iy,iz)
2306 h[hist+'XZ'].SetBinContent(iz,ix,N)
2307
2308

◆ muInterGeant4()

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

Definition at line 1131 of file muonDis.py.

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

◆ muondEdX()

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

Definition at line 1174 of file muonDis.py.

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

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

◆ muonDISProduction()

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

Definition at line 2233 of file muonDis.py.

2233def muonDISProduction(cycle,ecut=1.,strippedEvents=False):
2234 for run in range(1,11):
2235 for k in [0,1000]:
2236 prod = str(run+cycle*100+k)
2237 if strippedEvents: fn = "muonDis_"+prod+"_stripped.root"
2238 else: fn = "muonDis_"+prod+".root"
2239 command = "python $SNDSW_ROOT/macro/run_simScript.py --shiplhc"
2240 os.system(command+" -f "+fn+" --MuDIS --output run_"+prod+" -n 99999999 --eMin "+str(ecut) + " &")
2241 while count_python_processes('run_simScript')>ncpus: time.sleep(200)
2242

◆ muonPreTransport()

muonDis.muonPreTransport ( )

Definition at line 182 of file muonDis.py.

182def muonPreTransport():
183 f = ROOT.TFile(options.muonIn)
184 nt = f.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 801 of file muonDis.py.

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

◆ 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 1525 of file muonDis.py.

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

◆ 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 2263 of file muonDis.py.

2263def selectEvents(Ecut=10):
2264 for cycle in range(10):
2265 for run in range(1,11):
2266 for k in [0,1000]:
2267 prod = str(run+cycle*100+k)
2268 f = ROOT.TFile('muonDis_'+str(prod)+'.root')
2269 sTree = f.DIS
2270 fstripped = ROOT.TFile('muonDis_'+str(prod)+'_stripped.root','recreate')
2271 newTree = sTree.CloneTree(0)
2272 for sTree in f.DIS:
2273 flag = False
2274 for P in sTree.Particles:
2275 if P.GetPdgCode() in [130,2112,-2112]:
2276 if P.Energy() > Ecut: rc = newTree.Fill()
2277 newTree.AutoSave()
2278 f.Close()
2279 print('f closed',prod)
2280 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 2281 of file muonDis.py.

2281def signalNeutrinos():
2282 nudict = {'e':'$\\nu_e$','ae':'$\overline{\\nu}_e$','mu':'$\\nu_\mu$','amu':'$\overline{\\nu}_\mu$','tau':'$\\nu_\\tau$','atau':'$\overline{\\nu}_\\tau$'}
2283 ut.readHists(h,'/mnt/hgfs/microDisk/CERNBOX/SND@LHC/Neutrinos/SND_Neutrinos_Interacting_CCDIS.root')
2284 for p in ['','a']:
2285 for nu in [ 'e','mu','tau']:
2286 line = ''
2287 hname = 'h'+p+'nu'+nu
2288 ut.makeIntegralDistrib(h,hname)
2289 hist = h['I-'+hname]
2290 line += nudict[p+nu]
2291 for E in [10,100,200,300,500,1000]:
2292 n = hist.FindBin(E)
2293 N = hist.GetBinContent(n)
2294 line += ' & '+'$%5.1F$'%(N)
2295 line +=' \\\\'
2296 print(line)
2297

◆ thermNeutron()

muonDis.thermNeutron ( )

Definition at line 1778 of file muonDis.py.

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

Variable Documentation

◆ cycle

muonDis.cycle

Definition at line 2311 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 2207 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 2313 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 2221 of file muonDis.py.

◆ nEvents

muonDis.nEvents

Definition at line 2313 of file muonDis.py.

◆ nMult

muonDis.nMult

Definition at line 2311 of file muonDis.py.

◆ normalization

dict muonDis.normalization = {}

Definition at line 23 of file muonDis.py.

◆ nStart

muonDis.nStart

Definition at line 2311 of file muonDis.py.

◆ nucleon

muonDis.nucleon

Definition at line 2316 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 2313 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 2313 of file muonDis.py.

◆ rMax

muonDis.rMax

Definition at line 2313 of file muonDis.py.

◆ rMin

muonDis.rMin

Definition at line 2313 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 2217 of file muonDis.py.

◆ sMax

muonDis.sMax

Definition at line 2311 of file muonDis.py.

◆ sMin

muonDis.sMin

Definition at line 2311 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.