Skip to content

Commit 2923b6b

Browse files
jrstevenjlabjrstevenjlab
committed
initial commit of DIRC macros, can be culled later as needed
1 parent 8cd15be commit 2923b6b

19 files changed

+2627
-0
lines changed

‎dirc/README.md‎

Lines changed: 87 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,87 @@
1+
########## Instructions to setup DIRC environment ###########
2+
3+
Use 2018-dirc-commissioning branch for:
4+
-HDDS (https://github.com/JeffersonLab/hdds)
5+
-sim-recon (https://github.com/JeffersonLab/sim-recon)
6+
-HDGeant4 (https://github.com/JeffersonLab/HDGeant4)
7+
8+
To setup your own builds of with these branches use the following commands:
9+
10+
> setenv BUILD_SCRIPTS /group/halld/Software/build_scripts
11+
> mkdir builds
12+
> cd builds
13+
> $BUILD_SCRIPTS/my_sim-recon_jlab.sh
14+
15+
This will pull the master of each repository and build it, after this you'll need to switch to the proper branch
16+
17+
> cd hdds/
18+
> git checkout dirc_commissioning_2018
19+
> scons -u install
20+
> cd ../sim-recon/src/
21+
> git checkout dirc_commissioning_2018
22+
> scons -u install -j8
23+
> cd ../../hdgeant4/
24+
> git checkout dirc_commissioning_2018
25+
> make
26+
27+
Finally, you're environment can be set using the builds/setup_gluex.csh file creted by the build scripts. In order to use the local xml geometry with the DIRC you'll need to add the line below to your setup_gluex.csh script so the local geometry definition is used instead of CCDB
28+
29+
setenv JANA_GEOMETRY_URL xmlfile://${HDDS_HOME}/main_HDDS.xml
30+
31+
Now your environment is set to run DIRC simulations!
32+
33+
-------------------------------------------------------------------------
34+
35+
####### Example DIRC simulation workflow #######
36+
37+
#1) Set particle generator/source by using control.in file and run.mac files
38+
# Set KINE and SCAP cards for single particle MC in control.in
39+
# KINE controls momentum (GeV/c) and angle (degrees) of paticle gun
40+
# SCAP controls vertex postion of particle gun: should be (0., 0., 65.)
41+
# Set number of events to generate using run.mac
42+
# /run/beamOn 100: will generate 100 events
43+
44+
#2) Generate single track events and smear hits
45+
hdgeant4 run.mac
46+
mcsmear out.hddm
47+
48+
#3) Analyze smeared .hddm output with several plugins (Note: use -PTRIG:BYPASS=1 to skip trigger simulation)
49+
# pid_dirc: convert DIRC hits into custom .root format for LUT algorithm
50+
# truth_dirc: make diagnostic truth histograms from DIRC hits
51+
# monitoring_hists: make diagnostic histograms for other detectors (eg. check tracking is working)
52+
hd_root -o drc.root -PPLUGINS=pid_dirc,truth_dirc,monitoring_hists -PTRIG:BYPASS=1 out_smeared.hddm
53+
54+
#4) Plot DIRC's hit pattern:
55+
root loadlib.C drawHP.C'("drc.root")'
56+
57+
For a GUI with 3D visualization one need to have OpenGl and Qt libraries installed.
58+
Also -DG4UI_USE_EXECUTIVE flag needs to be added in GNUmakefile
59+
60+
####### LUT creation #######
61+
62+
#1) Generate photons from the end of a single bar for (1M events in run.mac)
63+
# DIRCLUT card in control.in controls which bar the LUT will be generated for
64+
hdgeant4 run.mac
65+
66+
#2) Create ROOT tree with LUT from lut_dirc plugin for a single bar:
67+
hd_root -nthreads=10 -o lut.root --plugin=lut_dirc out.hddm
68+
69+
#3) Averege LUT directions for a single bar:
70+
root loadlib.C glxlut_avr.C+'("lut.root")'
71+
72+
#3.5) Add step for generating LUT for all bars via batch system
73+
74+
#4) Combine averaged LUT tables from many bars into single file:
75+
root loadlib.C glxlut_add.C+'("lut/loopLut/lut_avr_*.root")'
76+
77+
####### LUT reconstruction (standalone macro) #######
78+
79+
#1) Simulate pion and kaon samples with fixed momentum and anglues (see example simulation workflow above). Then combine into a single ROOT file with hadd.
80+
81+
#2) Run the LUT reconstruction macro to show the pi/K separation power
82+
root loadlib.C reco_lut.C'("drc_kaon_pion.root")'
83+
84+
####### LUT reconstruction (JANA factory) #######
85+
86+
Work in progress...
87+

‎dirc/control.in‎

Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
cINFILE 'bggen/bggen.hddm'
2+
c OUTFILE 'sim_hprime2600_100k_notr_nodc.hddm'
3+
4+
RUNG 40000
5+
OUTFILE 'out.hddm'
6+
7+
TRIG 100
8+
9+
c generation of the DIRC lookup table (LUT)
10+
c 0 - LUT for all bars simultaneously
11+
c {0..47} - LUT for a given radiator
12+
cDIRCLUT 0
13+
14+
c particle momentum theta phi delta_momentum delta_theta delta_phi
15+
c KINE 1050 0.000000003 90. 0. 0. 0. 0.
16+
17+
c KINE 108 4. 11. -90. 0. 0. 0.
18+
c KINE 108 4. 11. -90. 0. 0. 0.
19+
20+
c KINE 114 4. 11. -90. 0. 0. 0.
21+
KINE 111 4. 4.0 90. 0. 0. 0.
22+
23+
24+
25+
c The SCAP card determines the vertex position for the particle gun.
26+
c vertex_x vertex_y vertex_z (default 0,0,0)
27+
cSCAP 293.8 10.43 585.8
28+
cSCAP 0. -8.83 584.9
29+
cSCAP 0. 10.43 584.9
30+
SCAP 0. 0. 65.
31+
32+
c vertex_extent_r vertex_extent_z
33+
c TGTWIDTH 0.0 0.
34+
35+
POSTSMEAR 0
36+
DELETEUNSMEARED 0
37+
38+
HALO 5e-5
39+
TRAJECTORIES 0
40+
41+
BGGATE -200. 200.
42+
BGRATE 4.80
43+
BGTAGONLY 1
44+
45+
TREFSIGMA 10.
46+
RNDM 3
47+
TOFMAX 1e-5
48+
49+
HADR 1
50+
CKOV 1
51+
LABS 1
52+
53+
ABAN 0
54+
DEBU 1 10 1000
55+
NOSECONDARIES 0
56+
57+
cBFIELDMAP 'Magnets/Solenoid/solenoid_1200A_poisson_20140520'
58+
c BFIELDTYPE 'NoField'
59+
PSBFIELDMAP 'Magnets/PairSpectrometer/PS_1.8T_20150513_test'
60+
61+
SAVEHITS 0
62+
SHOWERSINCOL 0
63+
DRIFTCLUSTERS 0
64+
65+
END

‎dirc/drawHP.C‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,34 @@
1+
#define glx__sim
2+
#include "/work/halld2/home/jrsteven/2018-dirc/builds/sim-recon/src/plugins/Analysis/pid_dirc/DrcEvent.h"
3+
#include "glxtools.C"
4+
5+
void drawHP(TString infile="drc.root"){
6+
if(!glx_initc(infile,1,"data/drawHP")) return;
7+
8+
DrcHit hit;
9+
for (Int_t e=0; e<glx_ch->GetEntries(); e++){
10+
glx_ch->GetEntry(e);
11+
for (Int_t t=0; t<glx_events->GetEntriesFast(); t++){
12+
glx_nextEventc(e,t,10);
13+
if(e > glx_ch->GetEntries()-2) cout<<"particle is "<<glx_names[glx_findPdgId(glx_event->GetPdg())]<<endl;
14+
if(glx_event->GetParent()>0) continue;
15+
for(Int_t h=0; h<glx_event->GetHitSize(); h++){
16+
hit = glx_event->GetHit(h);
17+
Int_t pmt = hit.GetPmtId();
18+
Int_t pix = hit.GetPixelId();
19+
TVector3 gpos = hit.GetPosition();
20+
Double_t time = hit.GetLeadTime();
21+
22+
// if(pmt != 36) continue;
23+
//if(pix != 0) continue;
24+
//cout<<"channel - "<<hit.GetChannel()<<endl;
25+
if(pmt<108) glx_hdigi[pmt]->Fill(pix%8, 7-pix/8);
26+
//if(pmt>=108) glx_hdigi[pmt-108]->Fill(pix%8, 7-pix/8);
27+
}
28+
}
29+
}
30+
glx_drawDigi("m,p,v\n",0);
31+
glx_canvasAdd(glx_cdigi);
32+
glx_canvasSave(1,0);
33+
34+
}

‎dirc/drawHitPats1.C‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,38 @@
1+
#define glx__sim
2+
#include "../../../../sim-recon/master/src/plugins/Analysis/pid_dirc/DrcEvent.h"
3+
#include "glxtools.C"
4+
5+
void drawHitPats1(double xmin=-40., double xmax=45., double ymin=-40., double ymax=45.){
6+
TString infile = "data/kaons/kaons.root";
7+
if(!glx_initc(infile,1,"data/drawHP")) return;
8+
9+
DrcHit hit;
10+
TVector3 hpos;
11+
int Neve = 0;
12+
if(Neve == 0) Neve = glx_ch->GetEntries();
13+
for (Int_t e=0; e<Neve; e++){
14+
glx_ch->GetEntry(e);
15+
for (Int_t t=0; t<glx_events->GetEntriesFast(); t++){
16+
glx_nextEventc(e,t,10);
17+
if(glx_event->GetParent()>0) continue;
18+
hpos = glx_event->GetPosition();
19+
if(hpos.X() < xmin || hpos.X() > xmax || hpos.Y() < ymin || hpos.Y() > ymax) continue;
20+
for(Int_t h=0; h<glx_event->GetHitSize(); h++){
21+
hit = glx_event->GetHit(h);
22+
Int_t pmt = hit.GetPmtId();
23+
Int_t pix = hit.GetPixelId();
24+
TVector3 gpos = hit.GetPosition();
25+
Double_t time = hit.GetLeadTime();
26+
if(pmt<108) glx_hdigi[pmt]->Fill(pix%8, 7-pix/8);
27+
//if(pmt>=108) glx_hdigi[pmt-108]->Fill(pix%8, 7-pix/8);
28+
}
29+
}
30+
}
31+
32+
glx_drawDigi("m,p,v\n",0);
33+
glx_cdigi->SaveAs(Form("data/drawHitPat/hp_x%.1f_y%.1f.png",xmin + (xmax-xmin)*0.5, ymin + (ymax-ymin)*0.5));
34+
35+
cout<<"x = "<<xmin + (xmax-xmin)*0.5<<", y = "<<ymin + (ymax-ymin)*0.5<<endl;
36+
cout<<"xmin = "<<xmin<<", xmax = "<<xmax<<", ymin = "<<ymin<<", ymax = "<<ymax<<endl;
37+
38+
}

‎dirc/drawMult.C‎

Lines changed: 97 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,97 @@
1+
#define glx__sim
2+
#include "../../../../sim-recon/master/src/plugins/Analysis/pid_dirc/DrcEvent.h"
3+
#include "glxtools.C"
4+
5+
//void drawMult(TString infile="data/out_kaons_100k.root"){
6+
//void drawMult(TString infile="data/sim_etaprime2300.root"){
7+
void drawMult(TString infile="data/sim_hprime2600_100k_tr_nodc.root"){
8+
//void drawMult(TString infile="data/sim_hprime2600_100k_notr_nodc.root"){
9+
10+
//void drawMult(TString infile="data/sim_hprime2600.root"){
11+
//void drawMult(TString infile="data/sim_hprime2600.root"){
12+
if(!glx_initc(infile,1,"data/drawHP")) return;
13+
14+
// gStyle->SetOptStat(0);
15+
TH2F *hPoint = new TH2F("hPoint",";x [cm]; y [cm]",200,-120,120,200,-120,120);
16+
TH1F *hMult = new TH1F("hMult",";detected photons [#]; [#]",500,0,500);
17+
TH1F *hEnergy = new TH1F("hEnergy",";photon energy [GeV]; [#]",100,0,10);
18+
19+
const auto nmax(20);
20+
double minx=-100, maxx=100;
21+
TH1F *hMultX[nmax];
22+
for(auto i=0; i<nmax; i++){
23+
hMultX[i] = new TH1F(Form("hMultX_%d",i),Form("hMultX_%d;x [cm]; stat [#]",i),300,0,300);
24+
}
25+
26+
TVector3 hpos,gpos;
27+
DrcHit hit;
28+
for (auto e=0; e<glx_ch->GetEntries(); e++){
29+
glx_ch->GetEntry(e);
30+
for (auto t=0; t<glx_events->GetEntriesFast(); t++){
31+
glx_nextEventc(e,t,100);
32+
if(glx_event->GetParent()>0) continue;
33+
hpos = glx_event->GetPosition();
34+
double x(hpos.X()), y(hpos.Y());
35+
hPoint->Fill(x, y);
36+
37+
int nhits=glx_event->GetHitSize();
38+
for(auto h=0; h<nhits; h++){
39+
hit = glx_event->GetHit(h);
40+
Int_t pmt = hit.GetPmtId();
41+
Int_t pix = hit.GetPixelId();
42+
gpos = hit.GetPosition();
43+
Double_t time = hit.GetLeadTime();
44+
if(pmt<108) glx_hdigi[pmt]->Fill(pix%8, 7-pix/8);
45+
//if(pmt>=108) glx_hdigi[pmt-108]->Fill(pix%8, 7-pix/8);
46+
hEnergy->Fill(hit.GetEnergy()*1E9);
47+
}
48+
49+
hMult->Fill(nhits);
50+
if(fabs(fabs(y)-12)<4){
51+
int xid = nmax*(x-minx)/(maxx - minx);
52+
if(xid>=0 && xid<nmax) hMultX[xid]->Fill(nhits);
53+
}
54+
}
55+
}
56+
57+
TGaxis::SetMaxDigits(2);
58+
glx_drawDigi();
59+
glx_canvasAdd(glx_cdigi);
60+
61+
glx_canvasSave(1,0);
62+
63+
TGaxis::SetMaxDigits(4);
64+
glx_canvasAdd("hPoint",500,500);
65+
hPoint->Draw("colz");
66+
67+
glx_canvasAdd("hMult",800,400);
68+
hMult->Draw();
69+
70+
glx_canvasAdd("hMultX");
71+
72+
TGraph *gMult = new TGraph();
73+
for(auto i=0; i<nmax; i++){
74+
double nph = glx_fit(hMultX[i],40,50,30).X();
75+
hMultX[i]->Draw();
76+
77+
double xpos = minx + 0.5*(maxx - minx)/nmax + i*(maxx - minx)/nmax;
78+
glx_waitPrimitive("hMultX");
79+
gMult->SetPoint(i,xpos,nph);
80+
}
81+
82+
glx_canvasAdd("gMultX",800,400);
83+
gMult->GetXaxis()->SetRangeUser(-110,110);
84+
gMult->GetYaxis()->SetRangeUser(0,150);
85+
gMult->GetXaxis()->SetTitle("x [cm]");
86+
gMult->GetYaxis()->SetTitle("detected photons [#]");
87+
gMult->SetMarkerStyle(20);
88+
gMult->SetMarkerSize(0.8);
89+
gMult->Draw("APL");
90+
91+
glx_canvasAdd("hEnergy");
92+
hEnergy->Draw();
93+
94+
glx_canvasSave(0);
95+
96+
}
97+

‎dirc/drawNpho.C‎

Lines changed: 73 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,73 @@
1+
#define glx__sim
2+
#include "TFile.h"
3+
#include "TTree.h"
4+
#include "../../../../sim-recon/master/src/plugins/Analysis/pid_dirc/DrcEvent.h"
5+
#include "glxtools.C"
6+
7+
// plots 1d histogram with N pho / track as a function of the theta angle
8+
// x range of interest is [-100, 100] cm.
9+
void drawNpho(TString infile="data/etaprime2300_all.root"){
10+
11+
const Int_t nbins = 50;
12+
13+
TH1F *hb[nbins];
14+
15+
for(Int_t i=0; i<nbins; i++){
16+
hb[i] = new TH1F(Form("hb_%d", i), Form("hb_%d", i), 100, 0., 200.);
17+
}
18+
19+
TH2F* hbary = new TH2F("hbary", "bar #; y [cm];", 49, -0.5, 48.5, 50, -100., 100.);
20+
TH1F* hNpho = new TH1F("hNpho","; x [cm]", nbins, -100.,100.);
21+
Double_t bin = 200./nbins;
22+
23+
if(!glx_init(infile,1,"data/drawNpho")) return;
24+
25+
DrcHit hit;
26+
for (Int_t ievent=0; ievent<glx_ch->GetEntries(); ievent++){
27+
glx_nextEvent(ievent,1);
28+
cout<<"parent = "<<glx_event->GetParent()<<endl;
29+
if(glx_event->GetParent() > 0) continue;
30+
31+
Double_t x = glx_event->GetPosition().X();
32+
Double_t y = glx_event->GetPosition().Y();
33+
cout<<"x = "<<x<<", y = "<<y<<endl;
34+
if(x < 110. && y < 110.){
35+
Int_t bar = glx_event->GetId();
36+
Int_t N = glx_event->GetHitSize();
37+
Int_t binnum = (Int_t)((x+100.)/bin) + 1;
38+
cout<<ievent<<" , x = "<<x<<", y = "<<y<<", N = "<<N<<", bin = "<<binnum<<", bar - "<<bar<<endl;
39+
if(/*(bar == 13 || bar == 12 || bar == 25 || bar == 26) &&*/ (binnum <= 100 && binnum > 0)) { // take only 4 middle bars
40+
hb[binnum]->Fill(N);
41+
}
42+
hbary->Fill(bar, y);
43+
44+
/* for(Int_t h=0; h<glx_event->GetHitSize(); h++){
45+
hit = glx_event->GetHit(h);
46+
Int_t pmt = hit.GetPmtId();
47+
Int_t pix = hit.GetPixelId();
48+
TVector3 gpos = hit.GetPosition();
49+
Double_t time = hit.GetLeadTime();
50+
if(pmt<102) glx_hdigi[pmt]->Fill(pix%8, 7-pix/8);
51+
}*/
52+
}
53+
}
54+
55+
for(Int_t i=0; i<nbins; i++){
56+
if(hb[i]->GetEntries() > 0.){
57+
hNpho->SetBinContent(i+1, hb[i]->GetMaximum());
58+
}
59+
}
60+
61+
gStyle->SetOptStat(0);
62+
TCanvas* c0 = new TCanvas("c0","c0", 500,500);
63+
hbary->Draw("colz");
64+
65+
TCanvas* c = new TCanvas("c","c",500,500);
66+
hNpho->Draw();
67+
/*
68+
TCanvas* ch = new TCanvas("ch","ch",500,500);
69+
for(Int_t j=0; j<nbins; j++){
70+
hb[j]->Draw();
71+
ch->WaitPrimitive();
72+
}*/
73+
}

0 commit comments

Comments
 (0)