Skip to content
Snippets Groups Projects
Commit f3175788 authored by Zhan Li's avatar Zhan Li Committed by lintao@ihep.ac.cn
Browse files

beam bkg

parent 1951a679
No related branches found
No related tags found
No related merge requests found
......@@ -2,4 +2,4 @@ build.*
build
spack*
./Generator/output/
./Generator/options/
./Generator/options/
\ No newline at end of file
<?xml version="1.0" encoding="UTF-8"?>
<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0"
xmlns:xs="http://www.w3.org/2001/XMLSchema"
xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd">
<info name="TDR_o1_v01"
title="CepC reference detctor for TDR"
author="Zhan Li"
url="http://cepc.ihep.ac.cn"
status="developing"
version="v01">
<comment>CepC reference detector simulation models used for TDR </comment>
</info>
<includes>
<gdmlFile ref="${DD4hepINSTALL}/DDDetectors/compact/elements.xml"/>
<gdmlFile ref="../CRD_common_v02/materials.xml"/>
</includes>
<define>
<constant name="world_size" value="10*m"/>
<constant name="world_x" value="world_size"/>
<constant name="world_y" value="world_size"/>
<constant name="world_z" value="world_size"/>
<include ref="${DD4hepINSTALL}/DDDetectors/compact/detector_types.xml"/>
</define>
<include ref="./TDR_Dimensions_v01_01.xml"/>
<!--old version, should be check/-->
<include ref="../CRD_common_v01/Beampipe_v01_02.xml"/>
<fields>
<field name="InnerSolenoid" type="solenoid"
inner_field="Field_nominal_value"
outer_field="0"
zmax="SolenoidCoil_half_length"
inner_radius="SolenoidCoil_center_radius"
outer_radius="Solenoid_outer_radius">
</field>
<field name="OuterSolenoid" type="solenoid"
inner_field="0"
outer_field="Field_outer_nominal_value"
zmax="SolenoidCoil_half_length"
inner_radius="Solenoid_outer_radius"
outer_radius="Yoke_barrel_inner_radius">
</field>
</fields>
</lccdd>
......@@ -34,7 +34,7 @@ rndmgensvc.Engine = rndmengine.name()
dsvc = k4DataSvc("EventDataSvc")
#geometry_option = "CepC_v4-onlyVXD.xml"
geometry_option = "CepC_v4_onlyTracker.xml"
#geometry_option = "CepC_v4.xml"
# geometry_option = "CepC_v4.xml"
geometry_path = os.path.join(os.getenv("DETCEPCV4ROOT"), "compact", geometry_option)
geosvc = GeomSvc("GeomSvc")
......
#!/usr/bin/env python
#Author: Zhan Li <lizhan@ihep.ac.cn>
#Created [2024-03-07 Thu 14:53]
import ROOT
import math
#read
#path1="/scratchfs/atlas/lizhan/cepc/CEPCSW/test-onlyVXD-single.root"
#path1="/scratchfs/atlas/lizhan/cepc/CEPCSW/test-detsim10.root"
def GetVXDHitMapFromFile(infile_name,out_file_name):
file = ROOT.TFile(infile_name)
tree = file.Get("events")
posx_VXD=[]
posy_VXD=[]
posz_VXD=[]
posr_VXD=[]
for entry in tree:
for elem in entry.VXDCollection:
x=elem.position.x
y=elem.position.y
z=elem.position.z
r = math.sqrt(x*x + y*y)
posx_VXD.append(x)
posy_VXD.append(y)
posz_VXD.append(z)
posr_VXD.append(r)
file.Close()
VXDHitMap = ROOT.TH2F("VXDHitMap","Hit Map of VXD",40,-125,125,20,0,80)
for i in range(0,len(posz_VXD)):
VXDHitMap.Fill(posz_VXD[i],posr_VXD[i])
ofile=ROOT.TFile(out_file_name,"Recreate")
VXDHitMap.Write()
ofile.Close()
def GetSITHitMapFromFile(infile_path,out_file_name):
file = ROOT.TFile(infile_path)
tree = file.Get("events")
posx_SIT=[]
posy_SIT=[]
posz_SIT=[]
posr_SIT=[]
for entry in tree:
for elem in entry.SITCollection:
x=elem.position.x
y=elem.position.y
z=elem.position.z
r = math.sqrt(x*x + y*y)
posx_SIT.append(x)
posy_SIT.append(y)
posz_SIT.append(z)
posr_SIT.append(r)
file.Close()
SITHitMap = ROOT.TH2F("SITHitMap","Hit Map of SIT",200,-2400,2400,20,120,320)
for i in range(0,len(posz_SIT)):
SITHitMap.Fill(posz_SIT[i],posr_SIT[i])
ofile=ROOT.TFile(out_file_name,"Recreate")
SITHitMap.Write()
ofile.Close()
print("test")
def GetSETHitMapFromFile(infile_path,out_file_name):
file = ROOT.TFile(infile_path)
tree = file.Get("events")
posx_SET=[]
posy_SET=[]
posz_SET=[]
posr_SET=[]
for entry in tree:
for elem in entry.SETCollection:
x=elem.position.x
y=elem.position.y
z=elem.position.z
r = math.sqrt(x*x + y*y)
posx_SET.append(x)
posy_SET.append(y)
posz_SET.append(z)
posr_SET.append(r)
file.Close()
SETHitMap = ROOT.TH2F("SETHitMap","Hit Map of SET",200,-2400,2400,20,1700,2000)
for i in range(0,len(posz_SET)):
SETHitMap.Fill(posz_SET[i],posr_SET[i])
ofile=ROOT.TFile(out_file_name,"Recreate")
SETHitMap.Write()
ofile.Close()
#draw
def DrawHitMap(in_file_name,TH2name,out_file_name):
file = ROOT.TFile(in_file_name)
TH2_HitMap=file.Get(TH2name)
c = ROOT.TCanvas("c", "c", 800, 800)
ROOT.gStyle.SetOptStat(0)
TH2_HitMap.GetXaxis().SetTitle("Z [mm]")
TH2_HitMap.GetXaxis().SetTitleOffset(0.77)
TH2_HitMap.GetXaxis().SetTitleSize(0.05)
TH2_HitMap.GetYaxis().SetTitle("R [mm]")
TH2_HitMap.GetYaxis().SetTitleOffset(0.8)
TH2_HitMap.GetYaxis().SetTitleSize(0.045)
TH2_HitMap.Draw("COL Z CJUST SAME")
c.SaveAs(out_file_name+".root")
c.SaveAs(out_file_name+".pdf")
path1="/scratchfs/atlas/lizhan/cepc/CEPCSW/test-SIT.root"
VXDOutFile="/scratchfs/atlas/lizhan/cepc/CEPCSW/VXDHM.root"
SITOutFile="/scratchfs/atlas/lizhan/cepc/CEPCSW/SITHM.root"
SETOutFile="/scratchfs/atlas/lizhan/cepc/CEPCSW/SETHM.root"
GetVXDHitMapFromFile(path1,VXDOutFile)
GetSITHitMapFromFile(path1,SITOutFile)
GetSETHitMapFromFile(path1,SETOutFile)
DrawHitMap(VXDOutFile,"VXDHitMap","Full_VXDHM_1evt")
DrawHitMap(SITOutFile,"SITHitMap","Full_SITHM_1evt")
DrawHitMap(SETOutFile,"SETHitMap","Full_SETHM_1evt")
\ No newline at end of file
......@@ -27,6 +27,10 @@ BeamBackgroundFileParserV1::BeamBackgroundFileParserV1(const std::string& filena
m_readTree->SetBranchAddress("cosy", &cosy);
m_readTree->SetBranchAddress("dp", &dp);
m_readTree->SetBranchAddress("dz", &dz);
m_readTree->SetBranchAddress("pid", &pid);
if(m_readTree->GetBranch("cosz")){
m_readTree->SetBranchAddress("cosz", &cosz);
}
m_rate = rate;
m_timewindow = timewindow;
......@@ -46,14 +50,17 @@ bool BeamBackgroundFileParserV1::load(IBeamBackgroundFileParser::BeamBackgroundD
// Now, we get a almost valid data
const double m2mm = 1e3; // convert from m to mm
result.pdgid = 11;
result.pdgid = pid;
result.x = x * m2mm;
result.y = y * m2mm;
result.z = (z+dz) * m2mm;
result.px = p * cosx;
result.py = p * cosy;
result.pz = p * std::sqrt(1-cosx*cosx-cosy*cosy);
if(m_readTree->GetBranch("cosz")){
result.pz = p * cosz;}
else{
result.pz = p * std::sqrt(1-cosx*cosx-cosy*cosy);}
result.mass = 0.000511; // assume e-/e+, mass is 0.511 MeV
......
......@@ -22,7 +22,8 @@ private:
double m_rate;
double m_timewindow;
double x, y, z, cosx, cosy, dz, dp;
double x, y, z, cosx, cosy, dz, dp, cosz;
int pid;
};
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment