You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
462 lines
20 KiB
462 lines
20 KiB
#include <string.h>
|
|
#include <stdio.h>
|
|
#include <iostream>
|
|
#include <vector>
|
|
#include <utility>
|
|
|
|
#include <TFile.h>
|
|
#include <TTree.h>
|
|
#include <TSystemDirectory.h>
|
|
#include <string.h>
|
|
#include <TFile.h>
|
|
#include <TTree.h>
|
|
#include <TH1.h>
|
|
#include <TH2.h>
|
|
#include <TNtuple.h>
|
|
#include <iostream>
|
|
#include <fstream>
|
|
#include <TRandom.h>
|
|
#include <TROOT.h>
|
|
#include <algorithm>
|
|
#include <assert.h>
|
|
#include <math.h>
|
|
#include <stdlib.h>
|
|
#include <time.h>
|
|
#include <sys/time.h>
|
|
#include <TCanvas.h>
|
|
#include <TF1.h>
|
|
#include <TProfile.h>
|
|
#include <TVectorD.h>
|
|
|
|
|
|
#include <cstdlib>
|
|
#include <cmath>
|
|
#include <vector>
|
|
|
|
#include <TApplication.h>
|
|
|
|
|
|
|
|
using namespace std;
|
|
|
|
|
|
int main(int argc, char **argv)
|
|
{
|
|
|
|
// Working directories
|
|
const char *dirname = "/work/leverington/beamprofilemonitor/hitdata/HIT_01-12-2016/with_timestamp/";
|
|
|
|
|
|
const char *listfile = NULL;
|
|
int ion = atoi(argv[1]);
|
|
if (ion==1) {listfile = "inlistE255_p";}
|
|
else if (ion==2) {listfile = "inlistE21_p";}
|
|
else if (ion==3) {listfile = "inlistE255_c";}
|
|
else if (ion==4) {listfile = "inlistE21_c";}
|
|
else {cerr<< argv[1] << " not allowed" << endl;
|
|
cerr << "specify (1, 2, 3, or 4) for E255_p, E21_p, E255_c, or E21_c" << endl;
|
|
return 1;
|
|
}
|
|
ifstream in;
|
|
ofstream ic1energyfile,ic2energyfile, beamenergyfile,beamenergyfile_iccorr, noisefile, focusfile,focusHITfile, skewfile, kurtfile,focusFitfile,chi2Fitfile;
|
|
ofstream ic1energyfile_diff, beamenergyfile_diff, focusFitfile_diff, posFitfile_diff, peakFitfile;
|
|
|
|
TFile *f;
|
|
TTree *tree;
|
|
Int_t nfiles = 0;
|
|
Int_t EList = 0;
|
|
Double_t Energy = 0.;
|
|
Double_t F1 = 0.;
|
|
Double_t Intensity = 0.;
|
|
Double_t SP = 0.;
|
|
char Filename[50];
|
|
Int_t id = 0;
|
|
Double_t signalmean = 0.;
|
|
Double_t signalmeanerror = 0.;
|
|
Double_t signalmean_corr = 0.;
|
|
Double_t signalmeanerror_corr = 0.;
|
|
Double_t signalmean_offsetcorr = 0.;
|
|
Double_t signalmeanerror_offsetcorr = 0.;
|
|
|
|
Double_t ic1mean = 0.;
|
|
Double_t ic1meanerror = 0.;
|
|
Double_t ic2mean = 0.;
|
|
Double_t ic2meanerror = 0.;
|
|
|
|
Double_t sp_air = 0.;
|
|
Double_t sp_ps = 0.;
|
|
Double_t sp2_air = 0.;
|
|
Double_t sp2_ps = 0.;
|
|
|
|
Double_t ic_beam_ratio = 0.;
|
|
Double_t ic_beam_ratioerror = 0.;
|
|
|
|
Double_t totalcurrent = 0.;
|
|
Double_t totaltime = 0.;
|
|
Double_t intcorr = 0.;
|
|
Double_t icbpmcorrelation = 0.;
|
|
|
|
Double_t noisemean = 0.;
|
|
Double_t noiseRMS = 0.;
|
|
Double_t focusmean = 0.;
|
|
Double_t focusRMS = 0.;
|
|
Double_t skewmean = 0.;
|
|
Double_t skewRMS = 0.;
|
|
Double_t kurtmean = 0.;
|
|
Double_t kurtRMS = 0.;
|
|
Double_t focusFitmean = 0.;
|
|
Double_t focusFitRMS = 0.;
|
|
Double_t chi2Fitmean = 0.;
|
|
Double_t chi2FitRMS = 0.;
|
|
Double_t peakFitmean = 0.;
|
|
Double_t peakFitRMS = 0.;
|
|
Double_t focusdiffFitmean = 0.;
|
|
Double_t focusdiffFitRMS = 0.;
|
|
Double_t focusdiffFitmeanerror = 0.;
|
|
Double_t focusdiffFitRMSerror = 0.;
|
|
Double_t posdiffFitmean = 0.;
|
|
Double_t posdiffFitRMS = 0.;
|
|
Double_t posdiffFitmeanerror = 0.;
|
|
Double_t posdiffFitRMSerror = 0.;
|
|
Double_t signaldiffmean = 0.;
|
|
Double_t signaldiffRMS = 0.;
|
|
Double_t signaldiffmeanerror = 0.;
|
|
Double_t signaldiffRMSerror = 0.;
|
|
Double_t ic1diffmean = 0.;
|
|
Double_t ic1diffRMS = 0.;
|
|
Double_t ic1diffmeanerror = 0.;
|
|
Double_t ic1diffRMSerror = 0.;
|
|
Double_t focusdiffcorr = 0.;
|
|
|
|
cout << "Trying to open"<< Form(": %s%s",dirname, listfile)<<"... " << endl;
|
|
|
|
in.open(Form("%s%s.prn",dirname, listfile));
|
|
|
|
ic1energyfile.open(Form("%s_ic1.txt",listfile ));
|
|
ic2energyfile.open(Form("%s_ic1ratio.txt",listfile ));
|
|
|
|
beamenergyfile.open(Form("%s_bpm1.txt",listfile ));
|
|
beamenergyfile_iccorr.open(Form("%s_bpm1ratio.txt",listfile ));
|
|
|
|
noisefile.open(Form("%s_noise.txt",listfile ));
|
|
focusfile.open(Form("%s_focus1.txt",listfile ));
|
|
focusHITfile.open(Form("%s_focusHIT.txt",listfile ));
|
|
skewfile.open(Form("%s_skew1.txt",listfile ));
|
|
kurtfile.open(Form("%s_kurt1.txt",listfile ));
|
|
focusFitfile.open(Form("%s_focusfit.txt",listfile ));
|
|
chi2Fitfile.open(Form("%s_chi2fit.txt",listfile ));
|
|
ic1energyfile_diff.open(Form("%s_ic1_diff.txt",listfile ));
|
|
beamenergyfile_diff.open(Form("%s_bpm1_diff.txt",listfile ));
|
|
focusFitfile_diff.open(Form("%s_focusfit_diff.txt",listfile ));
|
|
posFitfile_diff.open(Form("%s_posfit_diff.txt",listfile ));
|
|
peakFitfile.open(Form("%s_peakfit.txt",listfile ));
|
|
|
|
TF1 * f_sp_air = new TF1("f_sp_air","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in Air [MeV cm2 /g]
|
|
f_sp_air->SetParameters(303.746, -0.873697,1.01335); //protons between 50 and 300 MeV
|
|
|
|
TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
|
|
f_sp_ps->SetParameters(361.936, -0.892255, 1.19133); //protons between 50 and 220 MeV
|
|
|
|
TF1 * f_h2sp_air = new TF1("f_h2sp_air","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of helium in air [MeV cm2 /g]
|
|
f_h2sp_air->SetParameters(1236.51, -0.880225, 4.17041); //helium between 50 and 300 MeV
|
|
|
|
TF1 * f_h2sp_ps = new TF1("f_h2sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of helium in polystyrene [MeV cm2 /g]
|
|
f_h2sp_ps->SetParameters(1387.08, -0.882686,4.60833); //heelium between 50 and 300 MeV
|
|
|
|
TF1 * f_c12sp_air = new TF1("f_c12sp_air","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of carbon in air [MeV cm2 /g]
|
|
f_c12sp_air->SetParameters(13291.2, -0.925464,45.3876); //carbon between 80 and 480 MeV
|
|
|
|
TF1 * f_c12sp_ps = new TF1("f_c12sp_ps","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of carbon in polystyrene [MeV cm2 /g]
|
|
f_c12sp_ps->SetParameters(14538.9, -0.922423,49.2859); //carbon between 80 and 480 MeV
|
|
|
|
TF1 * f_o16sp_air = new TF1("f_o16sp_air","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of oxygen in air [MeV cm2 /g]
|
|
f_o16sp_air->SetParameters(24624.9, -0.925425,80.6828); //oxygen between 80 and 480 MeV
|
|
|
|
TF1 * f_o16sp_ps = new TF1("f_o16sp_ps","[0]*pow(x,[1])+[2]", 80,480); //stopping power of oxygen in polystyrene [MeV cm2 /g]
|
|
f_o16sp_ps->SetParameters(26465.6, -0.928309,88.6728); //oxygen between 80 and 480 MeV
|
|
|
|
///// ic1/SP intensity correction factor = 19.8E6+/-0.1E6 particles/s per nA/(Mevcm2/g)
|
|
|
|
if (!in) {
|
|
cout << Form("%s%s",dirname, listfile) << " not found." << endl;
|
|
return 0;}
|
|
cout << Form("%s%s",dirname, listfile) << " opened." << endl;
|
|
|
|
in >> EList >> Energy >> F1 >> Intensity >> SP >> Filename >> id;//prime the while loop
|
|
while (!in.eof()) {
|
|
TFile * file = TFile::Open(Form("%s%s",dirname,Filename));
|
|
TFile * outfile = TFile::Open(Form("plots/%s_plots",Filename),"RECREATE");
|
|
cout << Filename << endl;
|
|
// TH1D* signalDist = 0;
|
|
// file->GetObject("signalDist",signalDist);
|
|
// signalmean = signalDist->GetMean();
|
|
// signalmeanerror = signalDist->GetMeanError();
|
|
// cout << signalmean << " " << signalmeanerror << endl;
|
|
// delete signalDist;
|
|
|
|
// TCanvas *c1 = new TCanvas("c1");
|
|
TTree *data;
|
|
file->GetObject("newdata", data);
|
|
TH2D *h_beamSignal_channel = (TH2D*)file->Get("th2d_beamSignal_channel");
|
|
TH2D *h_bkg_channel = (TH2D*)file->Get("th2d_bkg_channel");
|
|
|
|
noisemean = h_bkg_channel->GetMean(2);
|
|
noiseRMS = h_bkg_channel->GetRMS(2);
|
|
TF1 * f1 = new TF1("f1", "gaus(0)",-5,5);
|
|
|
|
|
|
TVectorD *v = (TVectorD*)file->Get("icinfo");
|
|
// v->Print();
|
|
totalcurrent = (*v)[0];
|
|
totaltime = (*v)[1];
|
|
// data->Print();
|
|
TH1D* h_beamSignal_1 = new TH1D("h_beamSignal_1","h_beamSignal_1",510,-500,50000);
|
|
TH1D* h_ic1_1 = new TH1D("h_ic1_1","h_ic1_1",510,-50,500);
|
|
TH1D* h_ic2_1 = new TH1D("h_ic2_1","h_ic2_1",510,-50,500);
|
|
TH1D* h_beamSignal_diff = new TH1D("h_beamSignal_diff","h_beamSignal_diff",200,-1.5,1.5);
|
|
TH1D* h_ic1_diff = new TH1D("h_ic1_diff","h_ic1_diff",200,-0.5,0.5);
|
|
|
|
|
|
|
|
TH1D* h_ic_beam_ratio = new TH1D("h_ic_beam_ratio","h_ic_beam_ratio",500,0,500.);
|
|
TH2D* h_beamSignal_v_ic = new TH2D("h_beamSignal_v_ic","h_beamSignal_v_ic",510,-50,500,510,-500,50000);
|
|
TH2D* h_beamSignal_v_beamSignal2 = new TH2D("h_beamSignal_v_beamSignal2","h_beamSignal_v_beamSignal2",510,-500,50000,510,-500,50000);
|
|
|
|
TH1D* h_beamFocusX_1 = new TH1D("h_beamFocusX_1","h_beamFocusX_1",200,0,50.);
|
|
TH1D* h_beamSkewX_1 = new TH1D("h_beamSkewX_1","h_beamSkewX_1",200,0,50.);
|
|
TH1D* h_beamKurtX_1 = new TH1D("h_beamKurtX_1","h_beamKurtX_1",200,0,50.);
|
|
TH1D* h_beamFocusX_fit = new TH1D("h_beamFocusX_fit","h_beamFocusX_fit",200,0,50.);
|
|
TH1D* h_beamChi2_fit = new TH1D("h_beamChi2_fit","h_beamChi2_fit",200,0,50.);
|
|
TH1D* h_beamChi2_fit2 = new TH1D("h_beamChi2_fit2","h_beamChi2_fit2",200,0,50.);
|
|
TH1D* h_beamFocusX_corr = new TH1D("h_beamFocusX_corr","h_beamFocusX_corr",200,0.0,2.0);
|
|
|
|
// TH1D* h_beamFocusX_diff = new TH1D("h_beamFocusX_diff","h_beamFocusX_diff",200,-.4,0.4);
|
|
TH1D* h_beamFocusX_diff = new TH1D("h_beamFocusX_diff","h_beamFocusX_diff",2000,-20,20);
|
|
TH1D* h_beamPosX_diff = new TH1D("h_beamPosX_diff","h_beamPosX_diff",400,-5.,5.);
|
|
|
|
Double_t peak = 0.;
|
|
TH1D* h_beamPeakX_fit = new TH1D("h_beamPeakX_fit","h_beamPeakX_fit",1000,0,8000.);
|
|
|
|
|
|
|
|
|
|
|
|
if (id>=0&&id<=26){
|
|
sp_air = f_sp_air->Eval(Energy)*(Intensity/1.0E6);
|
|
sp_ps = f_sp_ps->Eval(Energy)*(Intensity/1.0E6);
|
|
sp2_air = f_sp_air->Eval(Energy);
|
|
sp2_ps = f_sp_ps->Eval(Energy);
|
|
}
|
|
else if (id>=27&&id<=52){
|
|
sp_air = f_h2sp_air->Eval(Energy)*(Intensity/1.0E6);
|
|
sp_ps = f_h2sp_ps->Eval(Energy)*(Intensity/1.0E6);
|
|
sp2_air = f_h2sp_air->Eval(Energy);
|
|
sp2_ps = f_h2sp_ps->Eval(Energy);
|
|
}
|
|
else if (id>=53&&id<=78){
|
|
sp_air = f_c12sp_air->Eval(Energy)*(Intensity/1.0E6);
|
|
sp_ps = f_c12sp_ps->Eval(Energy)*(Intensity/1.0E6);
|
|
sp2_air = f_c12sp_air->Eval(Energy);
|
|
sp2_ps = f_c12sp_ps->Eval(Energy);
|
|
}
|
|
else if (id>=79&&id<=99){
|
|
sp_air = f_o16sp_air->Eval(Energy)*(Intensity/1.0E6);
|
|
sp_ps = f_o16sp_ps->Eval(Energy)*(Intensity/1.0E6);
|
|
sp2_air = f_o16sp_air->Eval(Energy);
|
|
sp2_ps = f_o16sp_ps->Eval(Energy);
|
|
}
|
|
else {
|
|
sp_air = -1.;
|
|
sp_ps = -1.;
|
|
sp2_air = -1.;
|
|
sp2_ps = -1.;
|
|
}
|
|
|
|
|
|
// cout << Form("beamPeakX_fit*sqrt(2)*%F/2.3548",F1) << endl;
|
|
// data->Project("h_beamSignal_1","beamSignal_1","beamon==1&&beamSignal_1>50");
|
|
// data->Project("h_beamSignal_1", Form("beamPeakX_fit*sqrt(2)*%F/2.3548",F1) ,"beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
data->Project("h_beamSignal_1", "beamPeakX_fit*sqrt(2)*beamFocusX_fit/2.3548" ,"beamon==1&&ic1_1>0.001&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
|
|
//second layer shows a factor 1.08 more signal
|
|
data->Project("h_beamSignal_diff", "((beamPeakX_fit*sqrt(2)*beamFocusX_fit/2.3548) - (beamPeakX_fit2*sqrt(2)*beamFocusX_fit2/2.3548)/1.08)/(beamPeakX_fit*sqrt(2)*beamFocusX_fit/2.3548)" ,"beamon==1&&ic1_1>0.001&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
|
|
data->Project("h_ic1_diff", " 2*(ic1_1 -ic2_1)/(ic1_1 +ic2_1)" ,"beamon==1&&ic1_1>0.001&&ic2_1>0.001&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
|
|
data->Project("h_beamFocusX_1","beamFocusX_1","beamon==1&&ic1_1>0.05");
|
|
|
|
data->Project("h_beamFocusX_corr","(beamFocusX_fit/beamFocusX_fit2)","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
focusdiffcorr = h_beamFocusX_corr->GetMean();
|
|
data->Project("h_beamSkewX_1","beamSkewX_1","beamon==1&&ic1_1>0.005");
|
|
data->Project("h_beamKurtX_1","beamKurtX_1","beamon==1&&ic1_1>0.005");
|
|
data->Project("h_beamFocusX_fit","beamFocusX_fit","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
data->Project("h_beamChi2_fit","beamChi2_fit","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
// data->Project("h_beamFocusX_diff",Form("(beamFocusX_fit/%f-beamFocusX_fit2)/((beamFocusX_fit/%f+beamFocusX_fit2)/2)",focusdiffcorr,focusdiffcorr),"beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
data->Project("h_beamFocusX_diff",Form("beamFocusX_fit/%f-beamFocusX_fit2",focusdiffcorr),"beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
data->Project("h_beamPosX_diff","beamPosX_fit2+beamPosX_fit-128","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamChi2_fit2<200&&beamPeakX_fit>20&&beamPeakX_fit2>20");
|
|
data->Project("h_beamPeakX_fit","beamPeakX_fit","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
peakFitmean = h_beamPeakX_fit->GetMean();
|
|
peakFitRMS = h_beamPeakX_fit->GetRMS();
|
|
|
|
|
|
focusmean = h_beamFocusX_1->GetMean();
|
|
focusRMS = h_beamFocusX_1->GetRMS();
|
|
peak = h_beamFocusX_diff->GetEntries()*sqrt(2)*(h_beamFocusX_diff->GetRMS());
|
|
f1->SetParameters(200,h_beamFocusX_diff->GetMean(),h_beamFocusX_diff->GetRMS());
|
|
f1->SetParLimits(0,0,100000);
|
|
f1->SetParLimits(1,-10,10);
|
|
f1->SetParLimits(2,0,5);
|
|
|
|
|
|
|
|
h_beamFocusX_diff->Fit(f1,"");
|
|
// focusdiffFitmean = h_beamFocusX_diff->GetMean();
|
|
// focusdiffFitmeanerror = h_beamFocusX_diff->GetMeanError();
|
|
focusdiffFitmean =f1->GetParameter(1);
|
|
focusdiffFitmeanerror = f1->GetParError(1);
|
|
focusdiffFitRMS =f1->GetParameter(2);
|
|
focusdiffFitRMSerror =f1->GetParError(2);
|
|
// focusdiffFitRMS = h_beamFocusX_diff->GetRMS();
|
|
// focusdiffFitRMSerror = h_beamFocusX_diff->GetRMSError();
|
|
peak = h_beamPosX_diff->GetEntries()*sqrt(2)*(h_beamPosX_diff->GetRMS());
|
|
|
|
f1->SetParameters(peak, h_beamPosX_diff->GetMean(),h_beamPosX_diff->GetRMS());
|
|
h_beamPosX_diff->Fit(f1,"Q");
|
|
// posdiffFitmean = h_beamPosX_diff->GetMean();
|
|
// posdiffFitmeanerror = h_beamPosX_diff->GetMeanError();
|
|
// posdiffFitRMS = h_beamPosX_diff->GetRMS();
|
|
// posdiffFitRMSerror = h_beamPosX_diff->GetRMSError();
|
|
posdiffFitmean =f1->GetParameter(1);
|
|
posdiffFitmeanerror = f1->GetParError(1);
|
|
posdiffFitRMS =f1->GetParameter(2);
|
|
posdiffFitRMSerror =f1->GetParError(2);
|
|
|
|
skewmean = h_beamSkewX_1->GetMean();
|
|
skewRMS = h_beamSkewX_1->GetRMS();
|
|
kurtmean = h_beamKurtX_1->GetMean();
|
|
kurtRMS = h_beamKurtX_1->GetRMS();
|
|
focusFitmean = h_beamFocusX_fit->GetMean();
|
|
focusFitRMS = h_beamFocusX_fit->GetRMS();
|
|
chi2Fitmean = h_beamChi2_fit->GetMean();
|
|
chi2FitRMS = h_beamChi2_fit->GetRMS();
|
|
|
|
|
|
data->Project("h_ic1_1","ic1_1","ic1_1>0.01");
|
|
data->Project("h_ic2_1","ic2_1","ic2_1>0.01");
|
|
data->Project("h_ic1_diff","(ic1_1-ic2_1)/ic1_1","ic1_1>0.01&&ic2_1>0.01");
|
|
|
|
|
|
|
|
|
|
|
|
|
|
data->Project("h_beamSignal_v_ic","(beamPeakX_fit*sqrt(2)*beamFocusX_fit/2.3548):ic1_1","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
icbpmcorrelation = h_beamSignal_v_ic->GetCorrelationFactor();
|
|
|
|
data->Project("h_ic_beam_ratio","(beamPeakX_fit*sqrt(2)*beamFocusX_fit/2.3548)/ic1_1","beamon==1&&ic1_1>0.005&&beamChi2_fit<200&&beamPeakX_fit>20");
|
|
|
|
signalmean = h_beamSignal_1->GetMean();
|
|
signalmeanerror = h_beamSignal_1->GetMeanError();
|
|
|
|
|
|
peak = h_beamSignal_diff->GetEntries()*sqrt(2)*(h_beamSignal_diff->GetRMS());
|
|
f1->SetParameters(peak, h_beamSignal_diff->GetMean(), h_beamSignal_diff->GetRMS());
|
|
h_beamSignal_diff->Fit(f1,"Q");
|
|
// signaldiffmean = h_beamSignal_diff->GetMean();
|
|
// signaldiffmeanerror = h_beamSignal_diff->GetMeanError();
|
|
// signaldiffRMS = h_beamSignal_diff->GetRMS();
|
|
// signaldiffRMSerror = h_beamSignal_diff->GetRMSError();
|
|
signaldiffmean =f1->GetParameter(1);
|
|
signaldiffmeanerror = f1->GetParError(1);
|
|
signaldiffRMS =f1->GetParameter(2);
|
|
signaldiffRMSerror =f1->GetParError(2);
|
|
|
|
|
|
ic1mean = h_ic1_1->GetMean();
|
|
ic1meanerror = h_ic1_1->GetMeanError();
|
|
|
|
// peak = h_ic1_diff->GetEntries()*sqrt(2)*(h_ic1_diff->GetRMS());
|
|
//f1->SetParameters(peak, h_ic1_diff->GetMean(), h_ic1_diff->GetRMS());
|
|
// h_ic1_diff->Fit(f1,"");
|
|
ic1diffmean = h_ic1_diff->GetMean();
|
|
ic1diffmeanerror = h_ic1_diff->GetMeanError();
|
|
ic1diffRMS = h_ic1_diff->GetRMS();
|
|
ic1diffRMSerror = h_ic1_diff->GetRMSError();
|
|
// ic1diffmean =f1->GetParameter(1);
|
|
// ic1diffmeanerror = f1->GetParError(1);
|
|
// ic1diffRMS =f1->GetParameter(2);
|
|
// ic1diffRMSerror =f1->GetParError(2);
|
|
|
|
|
|
ic2mean = h_ic2_1->GetMean();
|
|
ic2meanerror = h_ic2_1->GetMeanError();
|
|
|
|
ic_beam_ratio=h_ic_beam_ratio->GetMean();
|
|
ic_beam_ratioerror=h_ic_beam_ratio->GetRMS();
|
|
// intcorr = (Intensity /19.8E6) / (totalcurrent/sp2_air/totaltime);
|
|
intcorr = (Intensity / 19.8E6) / (ic1mean/sp2_air)/ 3.11;
|
|
|
|
|
|
|
|
cout << id << " " <<sp2_air << " " << ic1mean/sp2_air << " " << intcorr << endl;
|
|
// cout << id <<" " << sp2_ps << " " <<signalmean/(Intensity/1.0E6)/(ic1mean/sp_air) << " " << signalmeanerror/(Intensity/1.0E6)/(ic1mean/sp_air) << endl;
|
|
|
|
ic1energyfile << Intensity/1.0E6*intcorr << " " << ic1mean *intcorr << " " << ic1meanerror*intcorr<< " " << Filename << endl;
|
|
// ic2energyfile << Intensity/1.0E6*intcorr << " " << ic2mean / (Intensity/1.0E6) / (0.659869 + 1.50842 * sp2_air) << " " << ic2meanerror/sp_air << " " << Filename << endl;
|
|
|
|
beamenergyfile << Intensity/1.0E6*intcorr << " " << signalmean*intcorr << " " << signalmeanerror*intcorr << endl;
|
|
beamenergyfile_iccorr << Intensity/1.0E6*intcorr << " " << (ic_beam_ratio*sp2_air/sp2_ps)/84.5 << " " << (ic_beam_ratioerror*sp2_air/sp2_ps)/84.5 << endl;
|
|
|
|
ic2energyfile << Intensity/1.0E6*intcorr << " " << icbpmcorrelation << endl;
|
|
noisefile << Intensity/1.0E6*intcorr << " " << noisemean << " " << noiseRMS << endl;
|
|
focusfile << Intensity/1.0E6*intcorr << " " << focusmean*0.8 << " " << focusRMS*0.8 << endl;
|
|
focusHITfile << Intensity/1.0E6*intcorr << " " << F1 << endl;
|
|
skewfile << Intensity/1.0E6*intcorr << " " << skewmean << " " << skewRMS << endl;
|
|
kurtfile << Intensity/1.0E6*intcorr << " " << kurtmean << " " << kurtRMS << endl;
|
|
focusFitfile << Intensity/1.0E6*intcorr << " " << focusFitmean*0.8 << " " << focusFitRMS*0.8 << endl;
|
|
chi2Fitfile << Intensity/1.0E6*intcorr << " " << chi2Fitmean << " " << chi2FitRMS << endl;
|
|
|
|
ic1energyfile_diff << Intensity/1.0E6*intcorr << " " << ic1diffRMS/sqrt(2) << " " << ic1diffRMSerror/sqrt(2) << " " << Filename << endl;
|
|
beamenergyfile_diff << Intensity/1.0E6*intcorr << " " << signaldiffRMS/sqrt(2) << " " << signaldiffRMSerror/sqrt(2) << endl;
|
|
// focusFitfile_diff << Intensity/1.0E6*intcorr << " " << focusdiffFitRMS*0.8/sqrt(2) << " " << focusdiffFitRMSerror*0.8/sqrt(2) << endl;
|
|
focusFitfile_diff << Intensity/1.0E6*intcorr << " " << focusdiffFitRMS/sqrt(2) << " " << focusdiffFitRMSerror/sqrt(2) << endl;
|
|
|
|
posFitfile_diff << Intensity/1.0E6*intcorr << " " << posdiffFitRMS*0.8/sqrt(2) << " " << posdiffFitRMSerror*0.8/sqrt(2) << endl;
|
|
peakFitfile << Intensity/1.0E6*intcorr << " " << (peakFitmean/2.)/6.7 << " " << (peakFitRMS/2.)/6.7 << endl;
|
|
|
|
|
|
|
|
|
|
// h_beamSignal_1->Draw("");
|
|
// c1->SaveAs(Form("figs/h_beamsignal_%i.pdf",nfiles));
|
|
// c1->Destructor();
|
|
|
|
// f->Close();
|
|
outfile->Write();
|
|
outfile->Close();
|
|
file->Close();
|
|
in >> EList >> Energy >> F1 >> Intensity >> SP >> Filename >> id;
|
|
// cout << EList << Energy << F1 << Intensity << SP << Filename << id << endl;
|
|
nfiles++;
|
|
|
|
}
|
|
ic1energyfile.close();
|
|
ic2energyfile.close();
|
|
|
|
beamenergyfile.close();
|
|
beamenergyfile_iccorr.close();
|
|
|
|
noisefile.close();
|
|
focusfile.close();
|
|
focusHITfile.close();
|
|
skewfile.close();
|
|
kurtfile.close();
|
|
in.close();
|
|
focusFitfile.close();
|
|
chi2Fitfile.close();
|
|
ic1energyfile_diff.close();
|
|
focusFitfile_diff.close();
|
|
beamenergyfile_diff.close();
|
|
posFitfile_diff.close();
|
|
return 0;
|
|
}
|