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

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

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

◆ boundaries()

muonDis.boundaries ( )

Definition at line 987 of file muonDis.py.

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

◆ checkProdofMuDIS()

muonDis.checkProdofMuDIS ( )

Definition at line 532 of file muonDis.py.

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

◆ compMuDIS_P6withP8()

muonDis.compMuDIS_P6withP8 ( )

Definition at line 553 of file muonDis.py.

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

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

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

◆ drawMuon3D()

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

Definition at line 1496 of file muonDis.py.

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

◆ energyOfDISevents()

muonDis.energyOfDISevents ( )

Definition at line 2246 of file muonDis.py.

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

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

◆ 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.Get("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 muPart_replica = ROOT.TParticle(muPart)
508 muPart_TCA = iMuon.ConstructedAt(0)
509 ROOT.std.swap(muPart_replica, muPart_TCA)
510 myPythia.GenerateEvent()
511# remove all unnecessary stuff
512 myPythia.Pyedit(2)
513 for itrk in range(1,myPythia.GetN()+1):
514 did = myPythia.GetK(itrk,2)
515 dpx,dpy,dpz = rotate(ctheta,stheta,cphi,sphi,myPythia.GetP(itrk,1),myPythia.GetP(itrk,2),myPythia.GetP(itrk,3))
516 psq = dpx**2+dpy**2+dpz**2
517 E = ROOT.TMath.Sqrt(getMasssq(did)+psq)
518 part = ROOT.TParticle(did,0,0,0,0,0,dpx,dpy,dpz,E,x,y,z,t)
519 part.SetWeight(w/float(options.nMult))
520# copy to branch
521 nPart = dPart.GetEntries()
522 if dPart.GetSize() == nPart: dPart.Expand(nPart+10)
523 part_TCA = dPart.ConstructedAt(nPart)
524 ROOT.std.swap(part, part_TCA)
525 nMade+=1
526 if nMade%options.heartbeat==0: print('made so far ',options.run,' :',nMade,time.ctime())
527 dTree.Fill()
528 fout.cd()
529 dTree.Write()
530 myPythia.SetMSTU(11, 6)
531 print("finished ",options.run)

◆ missing3Dproj()

muonDis.missing3Dproj (   hist,
  ymin,
  ymax 
)

Definition at line 2301 of file muonDis.py.

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

◆ muInterGeant4()

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

Definition at line 1134 of file muonDis.py.

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

◆ muondEdX()

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

Definition at line 1177 of file muonDis.py.

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

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

◆ muonDISProduction()

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

Definition at line 2236 of file muonDis.py.

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

◆ muonPreTransport()

muonDis.muonPreTransport ( )

Definition at line 182 of file muonDis.py.

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

◆ muonRateAtSND()

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

Definition at line 804 of file muonDis.py.

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

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

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

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

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

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

◆ thermNeutron()

muonDis.thermNeutron ( )

Definition at line 1781 of file muonDis.py.

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

Variable Documentation

◆ cycle

muonDis.cycle

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

◆ nEvents

muonDis.nEvents

Definition at line 2316 of file muonDis.py.

◆ nMult

muonDis.nMult

Definition at line 2314 of file muonDis.py.

◆ normalization

dict muonDis.normalization = {}

Definition at line 23 of file muonDis.py.

◆ nStart

muonDis.nStart

Definition at line 2314 of file muonDis.py.

◆ nucleon

muonDis.nucleon

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

◆ rMax

muonDis.rMax

Definition at line 2316 of file muonDis.py.

◆ rMin

muonDis.rMin

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

◆ sMax

muonDis.sMax

Definition at line 2314 of file muonDis.py.

◆ sMin

muonDis.sMin

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