Skip to content
Snippets Groups Projects
Commit 19ad3191 authored by Melanie Cardona's avatar Melanie Cardona Committed by Eric Pooser
Browse files

wrote script analyzing epics tree to calculate position of the beam at the target (#247)

* edited parameters file and wrote script analyzing epics tree to calculate position of the beam at the target

* moved EPICS macros to UTIL
parent dfe19d16
No related branches found
No related tags found
No related merge requests found
; BPM calibration constants
; =========================
guse_bpm_in_recon = 0 ; if 1 use bpm information for reconstruction
guse_bpmc = 1 ; if 1 use all 3 bpm info, if 0 use bpm A and B
gbpm_sample = 500 ; number of events, which are used for average
; beam position ( <5000 ). Optimal value dep. on rate
; the names are based on Paul Gueye's 'Status of the actual Beam Position
; Monitors in the Hall C Beamline', December 1, 1995.
; (mkj) apr-4-03 kappa,alpha and off for 3rd bpm are guesses,
; need to find right values.
;
; gbpm_kappa = 1.85 ,1.85,1.85 ; sensitivity in cm
gbpm_kappa = 2.00 ,1.95,1.85 ; modified mkj 4/9/03
gbpm_alpha_x = 1.8192,0.7330 ,0.7 ; calibration gain: we may get them from the
gbpm_alpha_y = 1.0063,0.8935 ,0.7 ; EPICS events, right now that's according to Paul's note
; the following offsets get added to the calculated positions
; gbpm_x_off = -0.011,-0.004,0. ; in cm: survey according to Paul's note
; gbpm_y_off = +0.052,+0.056,0.
; From Dahlberg survey DT_C853 Apr 7, 2003
gbpm_x_off = 0.017+.412-.094 , 0.061-.353-.055 ,0.037-0.09
gbpm_y_off = 0.037+.381+.003 ,-0.040-.211-.033 ,0.028+.219
;From Dave Gaskell's survey
;gbpm_x_off = -0.284, -0.288, -0.326
;gbpm_y_off = -0.021, -0.020, -0.0045
; average beam positions: only used if 'guse_bpm_in_recon=0'
; gbeam_xoff = +0.18
gbeam_xoff = +0.00
gbeam_xpoff = +0.00
gbeam_yoff = +0.00
gbeam_ypoff = +0.00
; spectrometers would like to see this positions (from optics studies)
; we treat SOS and HMS the same, however, we may change this if necessary
; gspec_xoff = +0.18
gspec_xoff = +0.00
gspec_xpoff = +0.00
gspec_yoff = +0.00
gspec_ypoff = +0.00
; Pedestals for BPM ADCs: from cosmic run #10933, 10/3/96, book XI-137
gbpm_xp_ped = 431.9,331.0,0. ; these pedestals have to be determined
gbpm_xm_ped = 514.7,350.6,0. ; from runs with no beam in the cavities,
gbpm_yp_ped = 406.9,358.8,0. ; e.g. cosmic runs.
gbpm_ym_ped = 499.7,292.7,0. ;
; positions of BPMs relative to target (from Paul's note)
gbpm_zpos = 345.5,166.3,0. ; cm
;positions of BPMs relative to target (from Mark 2017)
gbpm_zpos = 370.82, 224.96, 129.3 ; cm
; Fast Raster calibration constants
; =================================
; flag if 1 fast raster data used with average beam pos and angles in
calculating beam position
; if 0 then only use average beam pos and angles in calculating beam position
; Various fast raster quantities: gUse* are flags
gusefr = 1 ; if 1 correct for FRY in reconstruction
guse_frdefault = 1 ; if 1 do no phase correction (default)
gfr_cal_mom = 6.4 ; = beam momnetum during calibration run
; gfrx_adcpercm = 11328. ; = FR channels per cm deflection on target
; gfry_adcpercm = 10708. ; from run 9981, August 25, 1996, book X-22
;gfr_cal_mom = 2.038 ; = beam momnetum during calibration run
;gfrx_adcpercm = 4364.7 ; = FR channels per cm deflection on target
;gfry_adcpercm = 5471.9 ; from harp scan (H00A), June 19, 2004
gfrxa_adcpercm = 97666.7;
gfrxb_adcpercm = 97333.3;
gfrya_adcpercm = 113667;
gfryb_adcpercm = 114333;
gfrxa_adc_zero_offset = 67800;
gfrxb_adc_zero_offset = 69050;
gfrya_adc_zero_offset = 67700;
gfryb_adc_zero_offset = 67000;
gfrx_adcmax = 1000 ; ADC amplitude in channels.
gfrx_maxsize = 0.1 ; fast raster amplitude in centimeter.
gfr_cal_mom = 6.4 ; = beam momnetum during calibration run
gfry_adcmax = 1000 ; ADC amplitude in channels.
gfry_maxsize = 0.1 ; fast raster amplitude in centimeter.
gfrxa_adcpercm = 97666.7;
gfrxb_adcpercm = 97333.3;
gfrya_adcpercm = 113667;
gfryb_adcpercm = 114333;
; The latest FR phase analysis from spring '96 showed, that there is no
; measurable phase shift. During early running (E91-13, E89-12) the FRY-phase
; was determined to be 5.8 degree.
gfrxa_adc_zero_offset = 67800;
gfrxb_adc_zero_offset = 69050;
gfrya_adc_zero_offset = 67700;
gfryb_adc_zero_offset = 67000;
; positions of FR magnets relative to target
gfrx_dist =1375 ; cm
gfrx_dist = 1375 ; cm
gfry_dist = 1337 ; cm
protorootfile ../ROOTfiles/shms_replay_production_XXXXX_1000000.root
guicolor orange
canvassize 1000 1200
newpage 2 3
title EPICS BPM Tree
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07A.XRAW")
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07A.YRAW")
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07B.XRAW")
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07B.YRAW")
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07C.XRAW")
macro UTIL/BEAMLINE/epicsTree.C("IPM3H07C.YRAW")
newpage 2 1
title BPM Mean Positions for Run
macro UTIL/BEAMLINE/plot_epics.C("XZmean")
macro UTIL/BEAMLINE/plot_epics.C("YZmean")
newpage 4 4
title BPM Positions per Event
macro UTIL/BEAMLINE/plot_epics.C("grXZe0")
macro UTIL/BEAMLINE/plot_epics.C("grYZe0")
macro UTIL/BEAMLINE/plot_epics.C("grXZe1")
macro UTIL/BEAMLINE/plot_epics.C("grYZe1")
macro UTIL/BEAMLINE/plot_epics.C("grXZe2")
macro UTIL/BEAMLINE/plot_epics.C("grYZe2")
macro UTIL/BEAMLINE/plot_epics.C("grXZe3")
macro UTIL/BEAMLINE/plot_epics.C("grYZe3")
macro UTIL/BEAMLINE/plot_epics.C("grXZe4")
macro UTIL/BEAMLINE/plot_epics.C("grYZe4")
macro UTIL/BEAMLINE/plot_epics.C("grXZe5")
macro UTIL/BEAMLINE/plot_epics.C("grYZe5")
macro UTIL/BEAMLINE/plot_epics.C("grXZe6")
macro UTIL/BEAMLINE/plot_epics.C("grYZe6")
newpage 2 2
title Beam Position and Angle at Target
macro UTIL/BEAMLINE/target.C("BeamXt")
macro UTIL/BEAMLINE/target.C("BeamYt")
macro UTIL/BEAMLINE/target.C("BeamXAt")
macro UTIL/BEAMLINE/target.C("BeamYAt")
void epicsTree(TString leafname)
{
TFile *f = new TFile("../ROOTfiles/shms_replay_production_452_1000000.root");
TTree *E = (TTree*)gDirectory->Get("E");
E->Draw(leafname);
}
void plot_epics(TString graphname)
{
TFile *f = new TFile("../ROOTfiles/shms_replay_production_452_1000000.root");
TTree *E = (TTree*)gDirectory->Get("E");
auto nevents = E->GetEntries();
TH1D* hAX;
hAX = new TH1D("hAX","BPMA X",100, -5.0,0.0);
TH1D* hAY;
hAY = new TH1D("hAY","BPMA Y",100, 0.0,6.0);
TH1D* hBX;
hBX = new TH1D("hBX","BPMB X",100, -5.0,0.0);
TH1D* hBY;
hBY = new TH1D("hBY","BPMB Y",100, 0.0,6.0);
TH1D* hCX;
hCX = new TH1D("hCX","BPMC X",100, -5.0,0.0);
TH1D* hCY;
hCY = new TH1D("hCY","BPMC Y",100, 0.0,6.0);
// offsets from Gaskell's survey values
const Double_t AXSOF = -2.84;
const Double_t AYSOF = -0.21;
const Double_t BXSOF = -2.88;
const Double_t BYSOF = -0.20;
const Double_t CXSOF = -3.26;
const Double_t CYSOF = -0.045;
// extra offsets of BPM 07B calculated from line drawn between RAW-SOF BPM A and C points from run 452 with 500000 events
const Double_t BXos = 0.929;
const Double_t BYos = -0.588;
//get mean values of x,y position for each BPM for the run
E->Draw("IPM3H07A.XRAW>>hAX");
Double_t meanAX = hAX->GetMean() - AXSOF;
cout << "Mean A XPOS = " << meanAX << endl;
E->Draw("IPM3H07A.YRAW>>hAY");
Double_t meanAY = hAY->GetMean() - AYSOF;
cout << "Mean A YPOS = " << meanAY << endl;
E->Draw("IPM3H07B.XRAW>>hBX");
Double_t meanBX = hBX->GetMean() - BXSOF - BXos;
cout << "Mean B XPOS = " << meanBX << endl;
E->Draw("IPM3H07B.YRAW>>hBY");
Double_t meanBY = hBY->GetMean() - BYSOF - BYos;
cout << "Mean B YPOS = " << meanBY << endl;
E->Draw("IPM3H07C.XRAW>>hCX");
Double_t meanCX = hCX->GetMean() - CXSOF;
cout << "Mean C XPOS = " << meanCX << endl;
E->Draw("IPM3H07C.YRAW>>hCY");
Double_t meanCY = hCY->GetMean() - CYSOF;
cout << "Mean C YPOS = " << meanCY << endl;
Double_t BPMX[3] = {meanAX, meanBX, meanCX};
Double_t BPMY[3] = {meanAY, meanBY, meanCY};
Double_t BPMZ[3] = {3708.2, 2249.6, 1293};
TF1 *line = new TF1("line", "[0]*x+[1]",-10.0, 10.0);
line->SetParameters(10000,1000);
TGraph *XZmean;
XZmean = new TGraph(3,BPMX,BPMZ);
XZmean->SetName("XZmean");
gDirectory->GetList()->Add(XZmean); //need to add TGraph to current directory
TGraph *YZmean;
YZmean = new TGraph(3,BPMY,BPMZ);
YZmean->SetName("YZmean");
gDirectory->GetList()->Add(YZmean);
TGraph* grmean;
if(graphname == "XZmean"){
grmean = (TGraph*) gDirectory->Get("XZmean");
grmean->GetXaxis()->SetTitle("X Position (mm)");
grmean->GetYaxis()->SetTitle("Z Position (mm)");
grmean->SetTitle("BPM Mean X for Run");
grmean->GetXaxis()->SetTitleOffset(1.25);
grmean->GetYaxis()->SetTitleOffset(1.50);
grmean->Draw("AP*");
grmean->Fit("line");
Double_t slopeX = line->GetParameter(0);
Double_t intX = line->GetParameter(1);
Double_t meanBeamX = (-1)*intX/slopeX;
cout << "The mean beam X position is: " << meanBeamX << endl;
}
if(graphname == "YZmean"){
grmean = (TGraph*) gDirectory->Get("YZmean");
grmean->GetXaxis()->SetTitle("Y Position (mm)");
grmean->GetYaxis()->SetTitle("Z Position (mm)");
grmean->SetTitle("BPM Mean Y for Run");
grmean->GetXaxis()->SetTitleOffset(1.25);
grmean->GetYaxis()->SetTitleOffset(1.50);
grmean->Draw("AP*");
grmean->Fit("line");
Double_t slopeY = line->GetParameter(0);
Double_t intY = line->GetParameter(1);
Double_t meanBeamY = (-1)*intY/slopeY;
cout << "The mean beam Y position is: " << meanBeamY << endl;
}
Double_t AX[nevents];
Double_t AY[nevents];
Double_t BX[nevents];
Double_t BY[nevents];
Double_t CX[nevents];
Double_t CY[nevents];
Double_t *xE[nevents];
Double_t *yE[nevents];
Double_t slopeX[nevents];
Double_t slopeY[nevents];
Double_t intX[nevents];
Double_t intY[nevents];
Double_t BeamX[nevents];
Double_t BeamY[nevents];
//get values of x,y position for each BPM per event
for(Int_t i=0; i<nevents; i++){
E->GetEntry(i);
AX[i] = E->GetLeaf("IPM3H07A.XRAW")->GetValue(0) - AXSOF;
cout << "A XPOS Event " << i << " = " << AX[i] << endl;
AY[i] = E->GetLeaf("IPM3H07A.YRAW")->GetValue(0) - AYSOF;
cout << "A YPOS Event " << i << " = " << AY[i] << endl;
BX[i] = E->GetLeaf("IPM3H07B.XRAW")->GetValue(0) - BXSOF - BXos;
cout << "B XPOS Event " << i << " = " << BX[i] << endl;
BY[i] = E->GetLeaf("IPM3H07B.YRAW")->GetValue(0) - BYSOF - BYos;
cout << "B YPOS Event " << i << " = " << BY[i] << endl;
CX[i] = E->GetLeaf("IPM3H07C.XRAW")->GetValue(0) - CXSOF;
cout << "C XPOS Event " << i << " = " << CX[i] << endl;
CY[i] = E->GetLeaf("IPM3H07C.YRAW")->GetValue(0) - CYSOF;
cout << "C YPOS Event " << i << " = " << CY[i] << endl;
}
Int_t i;
Int_t j;
Int_t n = 3;
for(i=0; i<nevents; i++){
xE[i] = new Double_t[n];
yE[i] = new Double_t[n];
}
for(i=0; i<nevents; i++){
for(j=0; j<1; j++){
xE[i][j] = AX[i];
yE[i][j] = AY[i];
}
for(j=1; j<2; j++){
xE[i][j] = BX[i];
yE[i][j] = BY[i];
}
for(j=2; j<3; j++){
xE[i][j] = CX[i];
yE[i][j] = CY[i];
}
}
TGraph *grX[nevents];
for(i=0; i<nevents; i++){
grX[i] = new TGraph(n,xE[i],BPMZ);
grX[i]->SetName(TString::Format("grXZe%d",i));
gDirectory->GetList()->Add(grX[i]);
}
TGraph *grY[nevents];
for(i=0; i<nevents; i++){
grY[i] = new TGraph(n,yE[i],BPMZ);
grY[i]->SetName(TString::Format("grYZe%d",i));
gDirectory->GetList()->Add(grY[i]);
}
for (i=0; i<nevents; i++){
if(graphname == (TString::Format("grXZe%d",i))){
grmean = (TGraph*) gDirectory->Get(TString::Format("grXZe%d",i));
grmean->GetXaxis()->SetTitle("X Position (mm)");
grmean->GetYaxis()->SetTitle("Z Position (mm)");
grmean->SetTitle(TString::Format("BPM X Pos Event %d",i));
grmean->GetXaxis()->SetTitleOffset(1.25);
grmean->GetYaxis()->SetTitleOffset(1.50);
grmean->Draw("AP*");
grmean->Fit("line");
}
if(graphname == (TString::Format("grYZe%d",i))){
grmean = (TGraph*) gDirectory->Get(TString::Format("grYZe%d",i));
grmean->GetXaxis()->SetTitle("Y Position (mm)");
grmean->GetYaxis()->SetTitle("Z Position (mm)");
grmean->SetTitle(TString::Format("BPM Y Pos Event %d",i));
grmean->GetXaxis()->SetTitleOffset(1.25);
grmean->GetYaxis()->SetTitleOffset(1.50);
grmean->Draw("AP*");
grmean->Fit("line");
}
}
}
void target(TString histoname)
{
TFile *f = new TFile("../ROOTfiles/shms_replay_production_452_1000000.root");
const Int_t n = 3;
Int_t i;
Int_t j;
Double_t z[n] = {3708.2, 2249.6, 1293}; //positions from target in mm
TTree *E = (TTree*)gDirectory->Get("E");
auto nevents = E->GetEntries();
Double_t AX[nevents];
Double_t AY[nevents];
Double_t BX[nevents];
Double_t BY[nevents];
Double_t CX[nevents];
Double_t CY[nevents];
Double_t *xE[nevents];
Double_t *yE[nevents];
Double_t slopeX[nevents];
Double_t slopeY[nevents];
Double_t intX[nevents];
Double_t intY[nevents];
Double_t BeamX[nevents];
Double_t BeamY[nevents];
Double_t theta_xz[nevents]; // angle between x and z at target
Double_t theta_yz[nevents]; // angle between y and z at target
// offsets from Gaskell's survey values
const Double_t AXSOF = -2.84;
const Double_t AYSOF = -0.21;
const Double_t BXSOF = -2.88;
const Double_t BYSOF = -0.20;
const Double_t CXSOF = -3.26;
const Double_t CYSOF = -0.045;
//offsets of BPM 07B calculated from line drawn between RAW-SOF BPM A and C points from run 452 with 500000 events
const Double_t BXos = 0.929;
const Double_t BYos = -0.588;
//get values of x,y position for each BPM per event
for(i=0; i<nevents; i++){
E->GetEntry(i);
AX[i] = E->GetLeaf("IPM3H07A.XRAW")->GetValue(0) - AXSOF;
cout << "A XPOS Event " << i << " = " << AX[i] << endl;
AY[i] = E->GetLeaf("IPM3H07A.YRAW")->GetValue(0) - AYSOF;
cout << "A YPOS Event " << i << " = " << AY[i] << endl;
BX[i] = E->GetLeaf("IPM3H07B.XRAW")->GetValue(0) - BXSOF - BXos;
cout << "B XPOS Event " << i << " = " << BX[i] << endl;
BY[i] = E->GetLeaf("IPM3H07B.YRAW")->GetValue(0) - BYSOF - BYos;
cout << "B YPOS Event " << i << " = " << BY[i] << endl;
CX[i] = E->GetLeaf("IPM3H07C.XRAW")->GetValue(0) - CXSOF;
cout << "C XPOS Event " << i << " = " << CX[i] << endl;
CY[i] = E->GetLeaf("IPM3H07C.YRAW")->GetValue(0) - CYSOF;
cout << "C YPOS Event " << i << " = " << CY[i] << endl;
}
TF1 *line = new TF1("line", "[0]*x+[1]",-10.0, 10.0);
line->SetParameters(10000,10000);
for(i=0; i<nevents; i++){
xE[i] = new Double_t[n];
yE[i] = new Double_t[n];
}
for(i=0; i<nevents; i++){
for(j=0; j<1; j++){
xE[i][j] = AX[i]; //BPM A x positions per event to feed into TGraph array
yE[i][j] = AY[i]; //BPM A y positions per event to feed into TGraph array
}
for(j=1; j<2; j++){
xE[i][j] = BX[i];
yE[i][j] = BY[i];
}
for(j=2; j<3; j++){
xE[i][j] = CX[i];
yE[i][j] = CY[i];
}
}
TGraph *grX[nevents];
for(i=0; i<nevents; i++){
grX[i] = new TGraph(n,xE[i],z);
grX[i]->Fit("line");
slopeX[i] = line->GetParameter(0);
intX[i] = line->GetParameter(1);
BeamX[i] = (-1)*(intX[i]/slopeX[i]); //x position of the beam at the target
theta_xz[i] = TMath::Tan((BeamX[i])/z[0]); //x angle of beam relative to z axis
}
TGraph *grY[nevents];
for(i=0; i<nevents; i++){
grY[i] = new TGraph(n,yE[i],z);
grY[i]->Fit("line");
slopeY[i] = line->GetParameter(0);
intY[i] = line->GetParameter(1);
BeamY[i] = (-1)*(intY[i]/slopeY[i]); //y position of the beam at the target
theta_yz[i] = TMath::Tan((BeamY[i])/z[0]); //y angle of beam relative to z axis
}
TH1D *BeamXt = new TH1D("BeamXt","Beam X Position at Target",100, -12.0, 12.0);
TH1D *BeamYt = new TH1D("BeamYt","Beam Y Position at Target",100, -12.0, 12.0);
TH1D *BeamXAt = new TH1D("BeamXAt","Beam X Angle at Target",100, -12.0, 12.0);
TH1D *BeamYAt = new TH1D("BeamYAt","Beam Y Angle at Target",100, -12.0, 12.0);
TH1D* h1;
if(histoname.Contains("Xt")){
h1 = (TH1D*) gDirectory->Get("BeamXt");
for (i=0; i<nevents; i++) {
h1->Fill(BeamX[i]);
}
h1->GetXaxis()->SetTitle("Beam X Position (mm)");
h1->Draw();
}
if(histoname.Contains("Yt")){
h1 = (TH1D*) gDirectory->Get("BeamYt");
for (i=0; i<nevents; i++) {
h1->Fill(BeamY[i]);
}
h1->GetXaxis()->SetTitle("Beam Y Position (mm)");
h1->Draw();
}
if(histoname.Contains("XAt")){
h1 = (TH1D*) gDirectory->Get("BeamXAt");
for (i=0; i<nevents; i++) {
theta_xz[nevents];
h1->Fill(theta_xz[i]*1000); //convert to mrad
}
h1->GetXaxis()->SetTitle("Beam X Angle (mrad)");
h1->Draw();
}
if(histoname.Contains("YAt")){
h1 = (TH1D*) gDirectory->Get("BeamYAt");
for (i=0; i<nevents; i++) {
theta_xz[nevents];
h1->Fill(theta_yz[i]*1000); //convert to mrad
}
h1->GetXaxis()->SetTitle("Beam Y Angle (mrad)");
h1->Draw();
}
}
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