leverington/hit_analyse_v2
1
0
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity
master
hit_analyse_v2/Scripts_20161126/beammon.h_old
Blake Leverington af0b3f1ed4 all files first push
2021-01-25 12:39:07 +01:00

2171 lines
60 KiB
Plaintext
Raw Permalink Blame History

#include <TFile.h>
#include <TTree.h>
#include <stdio.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TProfile.h>
#include <TCanvas.h>
#include <TF1.h>
#include <TMath.h>
#include <vector>
#include <utility>
#include <TString.h>
#include <TDirectory.h>
#include <TString.h>
#include <TVectorD.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <cmath>
#include <TPaveStats.h>
#include <TText.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TObjArray.h>
#include <TLegend.h>
using namespace std;
double globalUpperBound = 63;
double globalLowerBound = 1;
double measurementsPerSecond = 3000;
int Fit_Fail_Counts = 0;
int signalScale = 1;
int alignint = 1;
int alignintIc = 1;
double discardUpperThreshold = 0.01;
double discardLowerThreshold = -0.015;
Double_t fitfunc(Double_t *x, Double_t *par)
{
Double_t arg1 = 0;
if (par[2] != 0)
arg1 = (x[0] - par[1]) / par[2];
Double_t fitval1 = par[0] * TMath::Exp(-0.5 * arg1 * arg1);
Double_t arg2 = 0;
if (par[4] != 0)
arg2 = (x[0] - par[1]) / par[4];
Double_t fitval2 = par[3] * TMath::Exp(-0.5 * arg2 * arg2);
Double_t fitsum = par[5] + fitval1 + fitval2;
return fitsum;
}
class Beammon
{
private:
TTree *data;
TTree *newdata;
TBranch *b_time;
TBranch *b_ic1;
TBranch *b_ic2;
TBranch *b_mw1_focusx;
TBranch *b_mw1_focusy;
TBranch *b_mw2_focusx;
TBranch *b_mw2_focusy;
TBranch *b_mw1_posx;
TBranch *b_mw1_posy;
TBranch *b_mw2_posx;
TBranch *b_mw2_posy;
TBranch *b_channels[64];
TBranch *bb_ic1;
TBranch *bb_ic2;
TBranch *bb_mw1_focusx;
TBranch *bb_mw2_focusx;
TBranch *bb_mw1_posx;
TBranch *bb_mw2_posx;
TBranch *b_beamPosX_1;
TBranch *b_beamFocusX_1;
TBranch *b_beamSignal_1;
TBranch *b_beamPosX_2;
TBranch *b_beamFocusX_2;
TBranch *b_beamSignal_2;
TBranch *b_beamOn;
double beamPosX_1;
double beamFocusX_1;
double beamSignal_1;
double beamPosX_2;
double beamFocusX_2;
double beamSignal_2;
double beamOn;
double sigma_pos = 0.1;
char *file;
TFile *fileIn;
TFile *fileOut;
TObjArray objectArray;
bool isSipm;
bool firstMat;
double channelWidth;
double shift;
long nevents;
int nBaselineEvents;
int nPrelimEvents;
int highestChannelBin;
int highestTimeBin;
double beamOnAdcVal;
double referenceIntegral;
double beamOnLevel;
double prelimMean;
double prelimSigma;
double highVal;
double readOutFreq;
int intTime;
double chi2ndf;
double fwhm;
double mean_ave;
double sigma_ave;
double signalToParticles = 1; // multiply signal by this amount, so that mean of singal (photodiode) and icChamber distribution are the same
double signalDistScale = 1; // icChamber signal distribution has more entries than signal (photodiode) one, so increase bin content in signal dist. by 1*this factor, to achieve the same integral
double time;
int startTime;
int stopTime;
int timesteps;
double intSignal;
int beamOnTime;
double baselineIntegral;
double ic1;
double ic2;
double mw1_focusx;
double mw1_focusy;
double mw2_focusx;
double mw2_focusy;
double mw2_posx;
double mw2_posy;
double mw1_posx;
double mw1_posy;
double icIntegral;
double completeIntegral;
double ic1_1;
double ic2_1;
double mw1_focusx_1;
double mw2_focusx_1;
double mw2_posx_1;
double mw1_posx_1;
vector<pair<int, int> > spillTimes;
double channels[64];
double baseline[64];
double channelErrors[64];
double current;
TGraph *gChannelErrors;
TH1D *meanHistOn;
TH1D *widthHistOn;
TH2D *timePosHist;
TH1D *prelimBeam;
TH2D *timeCenters;
TH2D *timeCenters_MWPC;
TH2D *timeWidths;
TH2D *timeWidths_MWPC;
TH1D *signalDist;
TH1D *signalNoiseDist;
TH1D *signalNoiseChannel;
TH2D *simpleView;
TH1D *timeIntSignal;
TH1D *timeIntSignal_MWPC;
TH1D *beamPosHist;
TH1D *beamWidthHist;
TH1D *icSignalDist;
TGraphErrors *beamPosTimeGraph;
TGraphErrors *signalTimeGraph;
TGraphErrors *sigIcScatterGraph;
TGraphErrors *mwPosScatterGraph;
TH2D *sigIcScatterPlot[20];
TH2D *sigIcScatterPlot_ic2;
TH2D *mwPosScatterPlot[20];
TH1D *ic1Hist;
TH1D *ic2Hist;
TH1D *signalComp;
TH1D *positionComp;
TH1D *widthComp;
TH1D *mw1_focusxHist;
TH1D *mw1_focusyHist;
TH1D *mw2_focusxHist;
TH1D *mw2_focusyHist;
TH1D *mw1_posxHist;
TH1D *mw1_posyHist;
TH1D *mw2_posxHist;
TH1D *mw2_posyHist;
TH1D *channelSig[64];
TH1D *chisquareDist;
TH1D *mwPosX;
TH1D *mwFocusX;
TH1D *beamPosX;
TH1D *beamFocusX;
TH1D *sigHist;
TH1D *sigHist_weighted;
TH1D *fibreHist;
TH1D *fibreHist_weighted;
TH1D *mwpcHist;
TH2D *Center_Signal;
TH2D *Diff_timeCenters;
TH2D *Diff_Signal;
TH2D *Diff_ic1;
TH1D *Diff_hist;
TH2D *deltaSig_time;
TH2D *deltaSig_Sig;
TH1D *SlicesAvg;
TH2D *SignalErrHist;
void calcBaseline(bool useFirstEvents);
void calcBaselineIntegral();
void prelimProfile();
TH1D *AverageHist(int dir, int timeStart = -1, int timeStop = -1);
TH1D *TimeAverage(int timeStart, int timeStop);
TH1D *TimeIntSignal(int timeStart = -1, int timeStop = -1);
bool isAboveThreshold(TH1D *slice);
bool isAboveThreshold(double val);
bool isAboveThreshold(TGraphErrors *slice);
TH1D *calcTimeCenters(int start, int stop);
double calcBeamTime();
// Tom Add
public:
Beammon();
~Beammon();
void Initialize(const char *dirname, const char *filename, int baselineEvents = 1000, int prelimEvents = 10000, double beamLevel = 0.1, bool firstFibermat = true, double readOutFrequency = 1000., int integrationTime = 100);
TH1D *GetPrelimProfile();
void ProcessData();
TH1D *GetChannelAvg(bool beamOn = false);
TH1D *GetChannelSigDist(int channel);
TH1D *GetTimeAvg();
TH1D *GetTimeAvg(int nSpill);
TH1D *GetTimeAvg(int nSpill, int nMeasurements);
TH2D *GetTimeAvgScatter();
TH2D *GetTimeAvgScatter(int nSpill, int nMeasurements);
TH1D *GetChannelSlice(int time = -1, int offset = 0);
TH1D *GetChannelSlice2(int time = -1, int offset = 0);
TGraphErrors *GetChannelGraph(int time = -1, int offset = 0);
TH1D *GetChannelGraphHist(int time, int offset);
TH1D *GetTimeSlice(int channel = -1, int startTime = -1, int stopTime = -1, bool beamOn = false);
TCanvas *GetTimeSliceCanvas(int channel = -1);
TCanvas *GetSignalNoiseChannel(int channel = -1);
TH1D *GetSignalOverTime();
void GetTimeCenters();
void GetTimeWidths();
void PrintSpillTimes();
pair<int, int> GetSpillTimes(int nSpill);
TH2D *GetSpillCenters(int nSpill);
TH1D *GetSpillSlice(int nSpill);
TH1D *GetSignalDistribution();
TCanvas *GetSignalDistributionCanvas();
TH2D *GetSimpleView();
TH1D *GetBaselineHist();
TH1D *icSignalComp();
double GetIntSignal();
int GetBeamOnTime();
double GetSignalStd();
TCanvas *signalOverTime2d();
TCanvas *FitSnapshot(int offset);
TCanvas *FitAverage();
TCanvas *PrintBaseline();
TCanvas *PrintChannelNoise(int channel);
TCanvas *PrintProfile(int offset);
TCanvas *IcFitSlicesY();
TCanvas *PrintSigIcGraph();
TH1D *GetMwPosComp();
void Save();
void Close();
TCanvas *GetSig_e_ic_Hist();
TCanvas *GetSig_e_ic_Hist_Weighted();
int Fit_Fail_Number = 200;
int event_count = 0;
};
Beammon::Beammon() {}
Beammon::~Beammon() {}
void Beammon::Initialize(const char *dirname, const char *filename, int baselineEvents, int prelimEvents, double beamLevel, bool firstFibermat, double readOutFrequency, int integrationTime)
{
signalToParticles = 7937 / 1.443e4; // FOR RUN 26 only
signalDistScale = 1.; //44585./26225;// FOR RUN 26 ONLY
intTime = integrationTime;
readOutFreq = readOutFrequency;
firstMat = firstFibermat;
highestChannelBin = 20;
highestTimeBin = 0;
icIntegral = 0;
completeIntegral = 0;
TString fn = TString(filename);
TString dn = TString(dirname);
if (fn.Contains("sipm"))
{
channelWidth = 0.25;
isSipm = true;
}
else
{
channelWidth = 0.8;
isSipm = false;
}
TString inName = dn + "pin/" + fn;
fileIn = new TFile(inName, "UPDATE");
fileIn->GetObject("t", data);
TString outname = filename;
if (firstMat)
outname.Insert(0, "SAVE_");
else
outname.Insert(0, "SAVE_mat2_");
outname.Insert(0, dirname);
data->SetBranchStatus("ch*", 0);
fileOut = new TFile(outname.Data(), "RECREATE");
newdata = new TTree("newdata", "newdata");
data->SetBranchStatus("ch*", 1);
b_time = data->GetBranch("time");
b_time->SetAddress(&time);
for (int n = 0; n < 64; n++)
{
baseline[n] = 0.;
char bname[20];
if (firstMat)
{
sprintf(bname, "ch%02d", n);
b_channels[n] = data->GetBranch(bname);
b_channels[n]->SetAddress(&channels[n]);
}
else
{
sprintf(bname, "ch%d", n + 64);
printf("%s\n", bname);
b_channels[n] = data->GetBranch(bname);
b_channels[n]->SetAddress(&channels[63 - n]);
}
}
b_ic1 = data->GetBranch("ic1");
b_ic1->SetAddress(&ic1);
b_ic2 = data->GetBranch("ic2");
b_ic2->SetAddress(&ic2);
b_mw1_focusx = data->GetBranch("mw1_focusx");
b_mw1_focusx->SetAddress(&mw1_focusx);
b_mw1_focusy = data->GetBranch("mw1_focusy");
b_mw1_focusy->SetAddress(&mw1_focusy);
b_mw2_focusx = data->GetBranch("mw2_focusx");
b_mw2_focusx->SetAddress(&mw2_focusx);
b_mw2_focusy = data->GetBranch("mw2_focusy");
b_mw2_focusy->SetAddress(&mw2_focusy);
b_mw1_posx = data->GetBranch("mw1_posx");
b_mw1_posx->SetAddress(&mw1_posx);
b_mw1_posy = data->GetBranch("mw1_posy");
b_mw1_posy->SetAddress(&mw1_posy);
b_mw2_posx = data->GetBranch("mw2_posx");
b_mw2_posx->SetAddress(&mw2_posx);
b_mw2_posy = data->GetBranch("mw2_posy");
b_mw2_posy->SetAddress(&mw2_posy);
nevents = data->GetEntries();
startTime = 0;
stopTime = nevents - 1;
timePosHist = new TH2D("timePosHist", "beam position", 64, 0.5, 64.5, nevents, 0, nevents - 1); //display in seconds (y-axis), but with bins of 100ms
SignalErrHist = new TH2D("SignalErrHist", "SignalErrHist", 64, 0.5, 64.5, 200, 0, -1);
chisquareDist = new TH1D("chiqsuareDist", "chiqsuareDist", 100, -1, -1);
prelimBeam = new TH1D("prelimBeam", "preliminary beam profile", 64, 0.5, 64.5);
beamPosHist = new TH1D("beamPosHist", "mean positions", readOutFreq, 0., 64. * channelWidth);
beamWidthHist = new TH1D("beamWidthHist", "beam width", 500, 0., 50.);
beamPosTimeGraph = new TGraphErrors(nevents);
beamPosTimeGraph->SetName("beamPosTimeGraph");
signalTimeGraph = new TGraphErrors(nevents);
signalTimeGraph->SetName("signalTimeGraph");
sigIcScatterGraph = new TGraphErrors();
mwPosScatterGraph = new TGraphErrors();
//sigIcScatterPlot = new TH2D("sigIcScatterPlot", "", 100, 0, -1, 100, 0, -1);
signalComp = new TH1D("signalComp", "signalComp", 200., -1., 1.);
positionComp = new TH1D("positionComp", "positionComp", 200., -3., 3.);
widthComp = new TH1D("widthComp", "widthComp", 200., 0., 3.5);
ic1Hist = new TH1D("ic1Hist", "ic1Hist", nevents, 0., nevents - 1);
ic2Hist = new TH1D("ic2Hist", "ic2Hist", nevents, 0., nevents - 1);
sigHist = new TH1D("sigHist", "sigHist", nevents, 0., nevents - 1);
sigHist_weighted = new TH1D("sigHist_weighted", "sigHist_weighted", nevents, 0., nevents - 1);
fibreHist = new TH1D("fibreHist", "fibreHist", nevents, 0., nevents - 1);
fibreHist_weighted = new TH1D("fibreHist_weighted", "fibreHist_weighted", nevents, 0., nevents - 1);
mwpcHist = new TH1D("mwpcHist", "mwpcHist", nevents, 0., nevents - 1);
mw1_focusxHist = new TH1D("mw1_focusxHist", "mw1_focusxHist", nevents, 0, nevents - 1);
mw1_focusyHist = new TH1D("mw1_focusyHist", "mw1_focusyHist", nevents, 0, nevents - 1);
mw2_focusxHist = new TH1D("mw2_focusxHist", "mw2_focusxHist", nevents, 0, nevents - 1);
mw2_focusyHist = new TH1D("mw2_focusyHist", "mw2_focusyHist", nevents, 0, nevents - 1);
mw1_posxHist = new TH1D("mw1_posxHist", "mw1_posxHist", nevents, 0, nevents - 1);
mw1_posyHist = new TH1D("mw1_posyHist", "mw1_posyHist", nevents, 0, nevents - 1);
mw2_posxHist = new TH1D("mw2_posxHist", "mw2_posxHist", nevents, 0, nevents - 1);
mw2_posyHist = new TH1D("mw2_posyHist", "mw2_posyHist", nevents, 0, nevents - 1);
mwPosX = new TH1D("mwPosX", "", 200, 0, 20);
mwFocusX = new TH1D("mwFocusX", "", 200, 2, 3.5);
beamPosX = new TH1D("beamPosX", "", 200, 20, 30);
beamFocusX = new TH1D("beamFocusX", "", 200, 0, 20);
nBaselineEvents = baselineEvents;
nPrelimEvents = prelimEvents;
beamOnLevel = beamLevel;
beamOnAdcVal = 0;
beamOnTime = 0;
if (!isSipm)
{
// if (fn.Contains("int1"))
calcBaseline(false);
//else
// calcBaseline(true);
}
prelimProfile();
calcBaselineIntegral();
ProcessData();
}
void Beammon::Save()
{
fileIn->Close();
fileOut->cd();
for (int i = 0; i < 5; i++){
PrintProfile(i)->Write();
PrintProfile(i)->Delete();
}
FitAverage()->Write();
FitAverage()->Delete();;
PrintBaseline()->Write();
PrintBaseline()->Delete();//added
//for(int i=0;i<64;i++)
// PrintChannelNoise(i)->Write();
GetSignalOverTime()->Write();
GetSignalOverTime()->Delete();
for (int i = 0; i < 5; i++){
GetChannelSlice2(-1, i)->Write();
GetChannelSlice2(-1, i)->Delete();
}
//fileIn->Close();
GetTimeSliceCanvas(-1)->Write();
timePosHist->ProjectionX()->Write();
GetTimeCenters();
GetSignalDistributionCanvas()->Write();
GetSignalDistribution()->Write();
GetSignalNoiseChannel()->Write();
chisquareDist->Write();
sigIcScatterGraph->Write();
mwPosScatterGraph->Write();
for (int i = 6; i < 15; i++)
{
sigIcScatterPlot[i]->Write();
mwPosScatterPlot[i]->Write();
}
mwPosScatterPlot[11]->Write();
sigIcScatterPlot[11]->Write();
mw1_focusxHist->Write();
mw1_posxHist->Write();
mwPosX->Write();
mwFocusX->Write();
beamPosX->Write();
beamFocusX->Write();
TVectorD v(6);
v[0] = GetIntSignal();
v[1] = GetBeamOnTime();
v[2] = signalDist->GetStdDev();
v[3] = icIntegral;
v[4] = completeIntegral;
v[5] = fwhm;
v.Write("beaminfo");
beamPosTimeGraph->Write();
signalTimeGraph->Write();
gChannelErrors->Write();
ic1Hist->Write();
ic2Hist->Write();
sigHist->Write();
fibreHist->Write();
mwpcHist->Write();
fibreHist_weighted->Write();
sigHist_weighted->Write();
GetSig_e_ic_Hist()->Write();
GetSig_e_ic_Hist_Weighted()->Write();
signalComp->Write();
GetMwPosComp()->Write();
widthComp->Write();
icSignalComp()->Write();
signalOverTime2d();
IcFitSlicesY();
PrintSigIcGraph()->Write();
//GetTimeCenters();
newdata->AutoSave();
Close();
// v.Delete();
}
void Beammon::calcBaseline(bool useFirstEvents)
{
// use the first/last nBaselineEvents events to calculate the average/baseline
if (!useFirstEvents)
{ // last
for (int i = nevents - nBaselineEvents; i < nevents; i++)
{
data->GetEvent(i);
for (int ch = 0; ch < 64; ch++)
{
baseline[ch] += channels[ch] / nBaselineEvents;
}
}
}
else
{
for (int i = 1; i <= nBaselineEvents; i++)
{
data->GetEvent(i);
for (int ch = 0; ch < 64; ch++)
{
baseline[ch] += channels[ch] / nBaselineEvents;
}
}
}
}
TCanvas *Beammon::PrintBaseline()
{
TCanvas *c1 = new TCanvas();
c1->cd();
TGraphErrors *g = new TGraphErrors(64);
for (int k = 0; k < 64; k++)
{
double stdDev = channelSig[k]->GetStdDev();
channelErrors[k] = stdDev;
g->SetPoint(k, k, baseline[k]);
g->SetPointError(k, 0., channelErrors[k]);
}
g->SetTitle("");
g->SetName("baseline");
g->GetXaxis()->SetTitle("channel");
TString label;
label.Form("adc counts / %d #mus", intTime);
g->GetYaxis()->SetTitle(label.Data());
g->Draw("AP");
c1->SetName("baseline");
//delete g;//added
return c1;
}
void Beammon::calcBaselineIntegral()
{
/*baselineIntegral = 0;
for(int i=1; i <= nBaselineEvents;i++){
data->GetEvent(i);
double currentIntegral = 0;
for(int ch=0;ch<64;ch++){
currentIntegral += channels[ch] - baseline[ch];
}
if(currentIntegral > baselineIntegral)
baselineIntegral = currentIntegral;
}
printf("baseline int: %f \n",baselineIntegral);*/
double highVal = 0;
for (int i = 1; i <= nBaselineEvents; i++)
{
data->GetEvent(i);
double val = channels[highestChannelBin] - baseline[highestChannelBin];
if (val > highVal)
highVal = val;
}
beamOnAdcVal = beamOnLevel * highVal;
int maxNoise = floor(highVal);
for (int j = 0; j < 64; j++)
{
TString name;
name.Form("channelSig%d", j);
channelSig[j] = new TH1D(name.Data(), name.Data(), maxNoise * 2, -maxNoise, maxNoise);
}
for (int i = 1; i <= nBaselineEvents; i++)
{
data->GetEvent(i);
for (int ch = 0; ch < 64; ch++)
channelSig[ch]->Fill(channels[ch] - baseline[ch]);
}
gChannelErrors = new TGraph(64);
for (int k = 0; k < 64; k++)
{
double stdDev = channelSig[k]->GetStdDev();
channelErrors[k] = stdDev;
gChannelErrors->SetPoint(k + 1, k + 1, stdDev);
}
}
void Beammon::prelimProfile()
{
highVal = 0;
int histChannel = 0;
double highVal_tmp[64];
for (int i = 1; i < nevents - 1000; i += 1)
{ //nevents/nPrelimEvents){
data->GetEvent(i);
for (int ch = 0; ch < 64; ch++)
{
double val = channels[ch] - baseline[ch];
/* if (val > highVal_tmp[ch])
{
highVal_tmp[ch] = val;
}*/
if (isSipm)
histChannel = 64 - ch;
else
histChannel = ch + 1;
prelimBeam->SetBinContent(histChannel, prelimBeam->GetBinContent(histChannel) + val / nPrelimEvents); //nPrelimEvents);
}
}
TF1 *fit = new TF1("fitfunc", "gaus", 1, 64);
fit->SetParameter(1, prelimBeam->GetMean());
// printf("Mean: %f \n", prelimBeam->GetMean());
// printf("Sigma: %f \n", prelimBeam->GetStdDev());
// cout << prelimBeam->GetStdDev() << endl;
fit->SetParameter(2, prelimBeam->GetStdDev());
prelimBeam->Fit("fitfunc", "Q");
prelimMean = fit->GetParameter(1);
prelimSigma = fit->GetParameter(2);
highestChannelBin = floor(prelimMean);
if (highestChannelBin < 0)
highestChannelBin = 0;
/*
highVal = 0;
for (int i = 0; i < 64; i++)
{
if ( highVal_tmp[i] > highVal)
{
highVal = highVal_tmp[i];
}
}*/
// beamOnAdcVal = beamOnLevel * highVal;
delete fit;
printf("\ncenter channel: %d \n", highestChannelBin);
}
void Beammon::ProcessData()
{
bool beamOn = false;
bool previousBeamStatus = false;
int lastSpillOn = 0;
int lastSpillOff = 0;
int minSpillTime = 3000; // in ms
int histChannel;
int highestSigVal = 0;
for (int i = 1; i < nevents; i++)
{
if (i % 10000 == 00)
printf("Event: %d / %ld \n", i, nevents);
// cout << "Event: " << i << endl;
data->GetEvent(i);
double maxSignal = channels[highestChannelBin] - baseline[highestChannelBin];
/*
for(int j=0;j<64;j++){
double val = channels[j] - baseline[j];
if (val > maxSignal)
maxSignal = val;
}*/
if (isAboveThreshold(maxSignal) && (i > nBaselineEvents))
{
beamOn = true;
if ((i - lastSpillOff) > minSpillTime && !previousBeamStatus)
{
if (lastSpillOff != 0)
spillTimes.push_back(make_pair(lastSpillOn, lastSpillOff));
lastSpillOn = i;
}
}
else
{
beamOn = false;
if (previousBeamStatus)
lastSpillOff = i;
}
previousBeamStatus = beamOn;
int sigVal = 0;
for (int j = 0; j < 64; j++)
{
int chan;
if (isSipm)
{
if (j % 2 == 0)
histChannel = 64 - (j + 1);
else
histChannel = 64 - (j - 1);
}
else
histChannel = j + 1;
double chVal = channels[j] - baseline[j];
sigVal += chVal;
timePosHist->SetBinContent(histChannel, i, chVal);
completeIntegral += chVal;
}
if (sigVal > highestSigVal)
{
highestSigVal = sigVal;
highestTimeBin = i;
}
// cout << ic1 << " " << mw1_focusx << " " << mw1_focusy << endl;
// IONISATION CHAMBERS
double chargeToParts = 1. / 3.2886706587037837e-15 * 1e-9 * 312e-6;
if (!(ic1 != ic1))
ic1Hist->SetBinContent(i, ic1 * chargeToParts);
if (!(ic2 != ic2))
ic2Hist->SetBinContent(i, ic2 * chargeToParts);
// MW CHAMBERS
if (!(mw1_focusx != mw1_focusx))
mw1_focusxHist->SetBinContent(i, mw1_focusx);
mw1_focusyHist->SetBinContent(i, mw1_focusy);
mw2_focusxHist->SetBinContent(i, mw2_focusx);
mw2_focusyHist->SetBinContent(i, mw2_focusy);
if (!(mw1_posx != mw1_posx))
mw1_posxHist->SetBinContent(i, mw1_posx);
mw1_posyHist->SetBinContent(i, mw1_posy);
mw2_posxHist->SetBinContent(i, mw2_posx);
mw2_posyHist->SetBinContent(i, mw2_posy);
if (mw1_posx < 100)
{
mwPosX->Fill(mw1_posx);
mwFocusX->Fill(mw1_focusx);
}
}
spillTimes.push_back(make_pair(lastSpillOn, lastSpillOff));
TH1D *tmp = GetChannelSlice(-1);
tmp->Fit("gaus");
TF1 *fit = tmp->GetFunction("gaus");
referenceIntegral = fit->GetParameter(0); // this is the constant in front of the gaus!! The real integral is constant*sigma*sqrt(2*pi), but because sigma is the same in one measurement it doesn't matter!
signalDist = new TH1D("signalDist", "signal distribution", 200, 0., referenceIntegral * 2.51 * fit->GetParameter(2));
signalNoiseDist = new TH1D("signalNoiseDist", "", 200, 0., referenceIntegral * 2.51 * fit->GetParameter(2));
signalNoiseChannel = new TH1D("signalNoiseChannel", "", 200, 0., tmp->GetMaximum());
calcBeamTime();
printf("HighestBin: %d \n", highestTimeBin);
}
TH1D *Beammon::GetPrelimProfile()
{
if (!prelimBeam)
{
printf("WARNING: preliminary profile not yet generated");
}
return prelimBeam;
}
TH1D *Beammon::TimeAverage(int timeStart, int timeStop)
{
if (timeStart == -1)
timeStart = startTime;
if (timeStop == -1)
timeStop = stopTime;
int count = 0;
TH1D *res = new TH1D("timeAverage", "time average", 64, 0.5, 64.5);
TH1D *tmp;
printf("\n\n start: %d, stop: %d \n\n", timeStart, timeStop);
for (int i = timeStart; i <= timeStop; i++)
{
tmp = GetChannelSlice(i, 0);
if (isAboveThreshold(tmp))
{
res->Add(tmp, 1.);
count++;
}
delete tmp;
}
if (count)
{
for (int n = 1; n <= 64; n++)
{
res->SetBinContent(n, res->GetBinContent(n) / count);
}
return res;
}
else
return 0;
}
TH1D *Beammon::TimeIntSignal(int timeStart, int timeStop)
{
if (timeStart == -1)
timeStart = startTime;
if (timeStop == -1)
timeStop = stopTime;
TH1D *res = new TH1D("timeAverage", "time average", 64, 0.5, 64.5);
TH1D *tmp;
for (int i = timeStart; i < timeStop; i++)
{
tmp = GetChannelSlice(i, 0);
if (isAboveThreshold(tmp))
{
res->Add(tmp, 1.);
beamOnTime++;
}
delete tmp;
}
return res;
}
TH1D *Beammon::AverageHist(int dir, int timeStart, int timeStop)
{ // dir = 0 => x direction, 1 => y
TH2D *hist = timePosHist;
int nbins1, nbins2;
double binStart, binStop;
int firstBin, secondBin;
if (!dir)
{
// Timeaverage
if ((timeStart != -1) && (timeStop != -1))
{
firstBin = timeStart; //hist->GetYaxis()->FindBin((timeStart-startTime)/readOutFreq);
secondBin = timeStop; //hist->GetYaxis()->FindBin((timeStop-startTime)/readOutFreq);
nbins1 = hist->GetNbinsX();
nbins2 = secondBin - firstBin;
binStart = hist->GetXaxis()->GetBinCenter(1);
binStop = hist->GetXaxis()->GetBinCenter(nbins1);
}
else
{
nbins1 = hist->GetNbinsX();
nbins2 = hist->GetNbinsY();
binStart = hist->GetXaxis()->GetBinCenter(1);
binStop = hist->GetXaxis()->GetBinCenter(nbins1);
firstBin = 1;
secondBin = nbins2;
}
}
else
{
nbins1 = hist->GetNbinsY();
nbins2 = hist->GetNbinsX();
binStart = hist->GetYaxis()->GetBinCenter(1);
binStop = hist->GetYaxis()->GetBinCenter(nbins1);
firstBin = 1;
secondBin = hist->GetNbinsX();
}
TH1D *retHist = new TH1D("average", "average", nbins1, binStart, binStop);
for (int i = 1; i <= nbins1; i++)
{
int count = 0;
double avg = 0;
for (int j = firstBin; j <= secondBin; j++)
{
if (!dir)
{
/*bool beam = beamOnVector[j];
if(beam==false)
continue;
//else
//printf("TRUE");*/
avg += hist->GetBinContent(i, j);
}
else
avg += hist->GetBinContent(j, i);
count++;
}
if (count > 0)
avg /= count;
retHist->SetBinContent(i, avg);
}
return retHist;
}
TH1D *Beammon::GetTimeAvg()
{
int n = spillTimes.size();
TH1D *res = new TH1D("timeAverage", "time average", 64, 0.5, 64.5);
for (int i = 0; i < n; i++)
{
TH1D *tmp = TimeAverage(spillTimes[i].first, spillTimes[i].second);
if (!tmp)
printf("\n\n timeAvgHist is Null! \n\n");
if (tmp)
res->Add(tmp, 1. / n);
}
res->SetTitle("average beam profile");
return res;
}
TH1D *Beammon::GetTimeAvg(int nSpill)
{
if (nSpill >= spillTimes.size())
return 0;
TH1D *res = TimeAverage(spillTimes[nSpill].first, spillTimes[nSpill].second);
res->SetTitle("average beam profile");
return res;
}
TH1D *Beammon::GetTimeAvg(int nSpill, int nMeasurements)
{
if (nSpill >= spillTimes.size())
return 0;
TH1D *res = TimeAverage(spillTimes[nSpill].second - nMeasurements, spillTimes[nSpill].second);
TString title = TString::Format("average beam profile, %d ms, spill %d", nMeasurements, nSpill);
TString name = TString::Format("avg%dmsspill%d", nMeasurements, nSpill);
res->SetTitle(title.Data());
res->SetName(name.Data());
return res;
}
TH2D *Beammon::GetTimeAvgScatter(int nSpill, int nMeasurements)
{
int tStart = spillTimes[nSpill].first;
if (nMeasurements == -1)
tStart = spillTimes[nSpill].second - nMeasurements;
TH2D *res = new TH2D("tavgScatter", "tavgScatter", 64, 0.5, 64.5, 400, 0., highVal);
for (int i = spillTimes[nSpill].first; i < spillTimes[nSpill].second; i++)
{
TH1D *tmp = GetChannelSlice(i, 0);
if (!isAboveThreshold(tmp))
continue;
for (int bin = 1; bin <= 64; bin++)
{
res->Fill(bin, tmp->GetBinContent(bin));
}
delete tmp;
}
TString title = TString::Format("beam profile, %d ms, spill %d", nMeasurements, nSpill);
res->SetTitle(title.Data());
return res;
}
TH2D *Beammon::GetTimeAvgScatter()
{
TH2D *res = new TH2D("tavgScatter", "tavgScatter", 64, 0.5, 64.5, 400, 0., highVal);
int n = spillTimes.size();
for (int i = 0; i < n; i++)
{
res->Add(GetTimeAvgScatter(i, -1), 1. / n);
}
return res;
}
TH1D *Beammon::GetChannelSigDist(int channel)
{
TString sName = TString("signal distribution channel ") + TString(channel);
double maxAmp = 0;
for (int n = 0; n < spillTimes.size(); n++)
{
TH1D *tmp = GetSpillSlice(n);
if (tmp->GetMaximum() > maxAmp)
maxAmp = tmp->GetMaximum();
delete tmp;
}
TH1D *res = new TH1D(sName.Data(), sName.Data(), 200, 0., maxAmp);
for (int n = 0; n < spillTimes.size(); n++)
{
TH1D *tmp = GetSpillSlice(n);
for (int bin = 1; bin <= tmp->GetNbinsX(); bin++)
{
res->Fill(tmp->GetBinContent(bin));
//printf("bin: %d, amp: %f \n", bin, tmp->GetBinContent(bin));
}
delete tmp;
}
return res;
}
TH1D *Beammon::GetChannelAvg(bool beamOn)
{
return AverageHist(1, -1, -1);
}
TH1D *Beammon::GetTimeSlice(int channel, int timeStart, int timeStop, bool beamOn)
{
if (timeStart == -1)
timeStart = startTime;
if (timeStop == -1)
timeStop = stopTime;
int nbins = timeStop - timeStart;
TH1D *tmp2 = new TH1D("ts", "timeslice", nbins, (timeStart - startTime) / readOutFreq, (timeStop - startTime) / readOutFreq);
if (channel == -1)
channel = highestChannelBin;
TH1D *tmp = timePosHist->ProjectionY("", channel, channel);
int offset = tmp->GetXaxis()->FindBin(timeStart);
for (int bin = 1; bin <= nbins; bin++)
{
tmp2->SetBinContent(bin, tmp->GetBinContent(bin + offset));
}
TString title = TString::Format("time projection of channel %d", channel);
tmp2->SetTitle(title.Data());
return tmp2;
}
TCanvas *Beammon::GetTimeSliceCanvas(int channel)
{
TCanvas *c1 = new TCanvas();
c1->SetName("timeslice");
int nbins = stopTime - startTime;
TH1D *tmp2 = new TH1D("timeslice", "", nbins, 0., (stopTime - startTime) / readOutFreq);
if (channel == -1)
channel = highestChannelBin;
TH1D *tmp = timePosHist->ProjectionY("", channel, channel);
for (int bin = 1; bin <= nbins; bin++)
{
tmp2->SetBinContent(bin, tmp->GetBinContent(bin));
}
tmp2->GetXaxis()->SetTitle("times [s]");
TString label;
label.Form("adc counts / %d #mus", intTime);
tmp2->GetYaxis()->SetTitle(label.Data());
tmp2->Draw("P");
return c1;
}
TCanvas *Beammon::GetSignalNoiseChannel(int channel)
{
TCanvas *c1 = new TCanvas();
c1->SetName("signalNoiseChannel");
int nbins = stopTime - startTime;
if (channel == -1)
channel = highestChannelBin;
TH1D *tmp = timePosHist->ProjectionY("", channel, channel);
for (int bin = 1; bin <= nbins; bin++)
{
signalNoiseChannel->Fill(tmp->GetBinContent(bin));
}
signalNoiseChannel->GetYaxis()->SetTitle("entries");
TString label;
label.Form("adc counts / %d #mus", intTime);
signalNoiseChannel->GetXaxis()->SetTitle(label.Data());
signalNoiseChannel->Draw();
return c1;
}
TH1D *Beammon::GetChannelSlice(int time, int offset)
{
if (time == -1)
time = highestTimeBin + offset;
//printf("time: %d\n", time);
TH1D *res = timePosHist->ProjectionX("", time, time); //timePosHist->ProjectionX("",timePosHist->GetYaxis()->FindBin(time+offset),timePosHist->GetYaxis()->FindBin(time+offset));
//res->SetError(channelErrors);
res->SetTitle("single measurement (1ms)");
res->SetName("Channel_" + (TString)time);
res->Fit("gaus", "Q0");
double mean = res->GetFunction("gaus")->GetParameter(1);
double sigma = res->GetFunction("gaus")->GetParameter(2);
return res;
}
TH1D *Beammon::GetChannelSlice2(int time, int offset)
{
if (time == -1)
time = highestTimeBin + offset;
//printf("time: %d\n", time);
TH1D *res = timePosHist->ProjectionX("", time, time); //timePosHist->ProjectionX("",timePosHist->GetYaxis()->FindBin(time+offset),timePosHist->GetYaxis()->FindBin(time+offset));
//res->SetError(channelErrors);
res->SetTitle("single measurement (1ms)");
res->SetName("Channel_" + (TString)time);
res->Fit("gaus", "Q0");
double mean = res->GetFunction("gaus")->GetParameter(1);
mean_ave += mean;
double sigma = res->GetFunction("gaus")->GetParameter(2);
sigma_ave += sigma;
return res;
}
TGraphErrors *Beammon::GetChannelGraph(int time, int offset)
{
if (time == -1)
time = highestTimeBin + offset;
double xErr = 0.8 / sqrt(12);
TH1D *res = timePosHist->ProjectionX("", time, time);
TGraphErrors *graph = new TGraphErrors(64);
for (int ch = 1; ch <= 64; ch++)
{
graph->SetPoint(ch, ch, res->GetBinContent(ch));
// graph->SetPointError(ch,xErr,channelErrors[ch-1]);//sqrt(channelErrors[ch-1]*channelErrors[ch-1]+res->GetBinContent(ch)));
graph->SetPointError(ch, xErr, channelErrors[ch - 1]);
}
delete res;
return graph;
}
TH1D *Beammon::GetChannelGraphHist(int time, int offset)
{
if (time == -1)
time = highestTimeBin + offset;
double xErr = 0.8 / sqrt(12);
TH1D *res = timePosHist->ProjectionX("", time, time);
return res;
}
TH1D *Beammon::GetSignalOverTime()
{
TH1D *tmp = timePosHist->ProjectionY("", 1, 64);
tmp->SetTitle("signal over time");
tmp->SetName("signal over time");
return tmp;
}
void Beammon::PrintSpillTimes()
{
printf("%lu spills recorded: \n", spillTimes.size());
// cout << spillTimes.size() << "spills recorded: " << endl;
for (int n = 0; n < spillTimes.size(); n++)
printf("%02d: %f - %f \n", n, (spillTimes[n].first - startTime) / readOutFreq, (spillTimes[n].second - startTime) / readOutFreq);
}
TH1D *Beammon::calcTimeCenters(int start, int stop)
{
int tStart = start;
int tStop = stop;
if (start == -1)
tStart = startTime;
if (stop == -1)
tStop = stopTime;
int nbins = (tStop - tStart);
double mean = 0;
double sigma = 0;
double constant;
TH1D *res = new TH1D("timeIntSignal_tmp", "time integrated signal", 64, 0.5, 64.5);
TCanvas *tmpcanvas = new TCanvas();
int offset = tStart;
TF1 *fit = new TF1("gaus_L", "gaus+[3]", 1, 64);
TF1 *fit2 = new TF1("gaus_2", "gaus", 1, 64);
fit2->SetParameter(0, 2000);
fit2->SetParameter(1, mean_ave);
fit2->SetParameter(2, sigma_ave);
for (int bin = tStart; bin <= tStop; bin++)
{
TGraphErrors *tmp = GetChannelGraph(bin, 0);
TH1D *tmp2 = GetChannelGraphHist(bin, 0);
tmp->Fit("gaus_2", "Q0", "", mean_ave - 3 * sigma_ave, mean_ave + 3 * sigma_ave);
TF1 *func = tmp->GetFunction("gaus_2");
if (tmp2->GetMean() * 0.8 > 0)
fibreHist_weighted->SetBinContent(bin + 1, tmp2->GetMean() * 0.8);
else
fibreHist_weighted->SetBinContent(bin + 1, 0);
if (tmp2->GetEntries() > 0)
sigHist_weighted->SetBinContent(bin - 1, tmp2->GetEntries() * signalToParticles);
else
sigHist_weighted->SetBinContent(bin - 1, 0);
if (func)
{
if (!isAboveThreshold(tmp))
continue;
constant = func->GetParameter(0);
mean = func->GetParameter(1) * channelWidth;
sigma = func->GetParameter(2) * channelWidth;
if (mean <= 0 || mean > 64)
continue;
/*if (constant > 1e5 || mean > (mean_ave * channelWidth + 0.2) || mean < (mean_ave * channelWidth - 0.2))
{
fit->SetParameter(0, 1000);
fit->SetParameter(1, mean_ave);
fit->SetParameter(2, sigma_ave);
fit->SetParameter(3, 0.1);
fit->SetParLimits(0, 100, 2e5);
fit->SetParLimits(1, 1, 64);
fit->SetParLimits(2, 1, 10);
fit->SetParLimits(3, -50, 50);
tmp->Fit("gaus", "Q0", "", mean_ave - 3 * sigma_ave, mean_ave + 3 * sigma_ave);
func = tmp->GetFunction("gaus");
constant = func->GetParameter(0);
mean = func->GetParameter(1) * channelWidth;
sigma = func->GetParameter(2) * channelWidth;
}*/
if (abs(constant) > 1e5)
{
Fit_Fail_Number = Fit_Fail_Number - 1;
if (Fit_Fail_Number > 0)
{
tmpcanvas->cd();
tmp->Draw("AP");
func->Draw("SAME");
tmpcanvas->Modified();
tmpcanvas->Write();
}
continue;
}
//Fill the error hidto
double erry, tmpx;
bool discard = false;
int nPoint = 0;
for (int chn = 1; chn < 65; chn++)
{
tmp->GetPoint(nPoint, tmpx, erry);
SignalErrHist->Fill((double)tmpx, (erry - func->Eval(tmpx)) / (constant * sigma));
if ((erry - func->Eval(tmpx)) / (constant * sigma) > discardUpperThreshold || (erry - func->Eval(tmpx)) / (constant * sigma) < discardLowerThreshold)
{
tmp->RemovePoint(nPoint);
nPoint--;
discard = true;
}
nPoint++;
}
double signal = sigma * 2.51 * constant * signalToParticles;
//Fill the new tree
event_count++;
//data->GetEntry(bin);
beamPosX_1 = mean;
beamFocusX_1 = 2.35 * sigma;
beamSignal_1 = signal;
ic1_1 = ic1Hist->GetBinContent(bin - alignintIc);
ic2_1 = ic2Hist->GetBinContent(bin - alignintIc);
mw1_focusx_1 = mw1_focusxHist->GetBinContent(bin - alignint);
mw1_posx_1 = mw1_posxHist->GetBinContent(bin - alignint);
mw2_focusx_1 = mw2_focusxHist->GetBinContent(bin - alignint);
mw2_posx_1 = mw2_posxHist->GetBinContent(bin - alignint);
beamOn = 1;
newdata->Fill();
if (event_count % 5000 == 0)
{
printf("Signal events: %d ", event_count);
printf("Event number: %d\n", bin);
}
if (signal < 0 || signal != signal)
continue;
//Fill the Histograms
beamPosX->Fill(mean);
fibreHist->Fill(bin, mean);
mwpcHist->Fill(bin, mw1_posxHist->GetBinContent(bin - alignint));
beamFocusX->Fill(2.35 * sigma);
beamPosTimeGraph->SetPoint(bin, bin, mean);
beamPosTimeGraph->SetPointError(bin, 0., func->GetParError(1));
timeCenters->Fill((bin) / readOutFreq, mean);
timeWidths->Fill((bin) / readOutFreq, sigma * 2.35);
timeCenters_MWPC->Fill((bin) / readOutFreq, mw1_posxHist->GetBinContent(bin - alignint));
timeWidths_MWPC->Fill((bin) / readOutFreq, mw1_focusxHist->GetBinContent(bin - alignint));
Diff_timeCenters->Fill((bin) / readOutFreq, mean - mw1_posxHist->GetBinContent(bin - alignint) - shift);
if (mw1_posxHist->GetBinContent(bin - alignint) < 20)
{
mwPosScatterGraph->SetPoint(mwPosScatterGraph->GetN(), mw1_posxHist->GetBinContent(bin - alignint), mean);
//positionComp->Fill(mean - mw1_posxHist->GetBinContent(bin-4) +111.5415 - 32.02);
for (int i = 0; i < 20; i++)
mwPosScatterPlot[i]->Fill(mw1_posxHist->GetBinContent(bin - i + 10), mean);
}
double sigError = 2.51 * sqrt(sigma * sigma * func->GetParError(0) * func->GetParError(0) + constant * constant * func->GetParError(2) * func->GetParError(2)) * signalToParticles;
if (signal > 1e5 || sigError > signal)
{
//signal = 0;
sigError = 0;
}
if (signal)
{
signalDist->Fill(signal);
if (ic1Hist->GetBinContent(bin - alignintIc))
{
sigHist->Fill(bin - alignintIc - 1, signal);
signalComp->Fill((signal - ic1Hist->GetBinContent(bin - alignintIc)) / signal);
double x, y, totalSignal = 0, signalError = 0;
for (int k = 0; k < 63; k++)
{
tmp->GetPoint(k, x, y);
signalError += tmp->GetErrorY(k) * tmp->GetErrorY(k);
totalSignal += y * signalToParticles;
}
signalError = sqrt(signalError) * signalToParticles;
sigIcScatterGraph->SetPoint(sigIcScatterGraph->GetN(), ic1Hist->GetBinContent(bin - alignintIc), totalSignal);
sigIcScatterGraph->SetPointError(sigIcScatterGraph->GetN(), 0, signalError);
for (int i = 0; i < 20; i++)
sigIcScatterPlot[i]->Fill(ic1Hist->GetBinContent(bin - i + 10), signal);
Center_Signal->Fill(signal, mean);
// sigIcScatterPlot[14]->Fill(ic1Hist->GetBinContent(bin - 14 + 10), signal);
sigIcScatterPlot_ic2->Fill(ic2Hist->GetBinContent(bin - 14 + 10), signal);
}
}
else
{
//ic1Hist->SetBinContent(bin, 0);
}
widthComp->Fill(sigma * 2.51 - mw1_focusxHist->GetBinContent(bin - alignint));
//Center_Signal->Fill(totalSignal, mean);
Diff_hist->Fill(mean - mw1_posxHist->GetBinContent(bin - alignint) - shift);
Diff_Signal->Fill(signal, mean - mw1_posxHist->GetBinContent(bin - alignint) - shift);
Diff_ic1->Fill(ic1Hist->GetBinContent(bin - alignintIc), mean - mw1_posxHist->GetBinContent(bin - alignint) - shift);
icIntegral += ic1Hist->GetBinContent(bin);
signalTimeGraph->SetPoint(bin, bin, signal);
signalTimeGraph->SetPointError(bin, 0., sigError);
beamOnTime++;
double c2 = func->GetChisquare() / func->GetNDF();
if (c2 < 100)
chisquareDist->Fill(c2);
}
delete func;
delete tmp;
delete tmp2;
}
delete tmpcanvas;
return res;
}
void Beammon::GetTimeCenters()
{
TH1D *tmp_TimeCenter;
double mean_MWPC = mwPosX->GetMean();
shift = mean_ave * 0.8 / 5 - mwPosX->GetMean();
double max_ic1 = ic1Hist->GetMaximum();
beamOnTime = 0;
int nbins = stopTime - startTime;
mean_ave = prelimBeam->GetMean();
sigma_ave = prelimBeam->GetStdDev();
printf("Mean: %f \n", mean_ave);
printf("Sigma: %f \n", sigma_ave);
timeCenters = new TH2D("timeCenters", "", nbins / 100, 0., (stopTime - startTime) / readOutFreq, 200, (0.8 * mean_ave - 6. * sigma_pos), (0.8 * mean_ave + 6 * sigma_pos));
timeWidths = new TH2D("timeWidths", "", nbins / 100, 0., (stopTime - startTime) / readOutFreq, 200, 0.4 * sigma_ave * 2.35, 0.95 * sigma_ave * 2.35);
timeIntSignal = new TH1D("timeIntSignal", "time integrated signal", 64, 0.5, 64.5);
timeCenters_MWPC = new TH2D("timeCenters_MWPC", "", nbins / 100, 0., (stopTime - startTime) / readOutFreq, 200, mean_MWPC - 5 * sigma_pos, mean_MWPC + 5 * sigma_pos);
timeWidths_MWPC = new TH2D("timeWidths_MWPC", "", nbins / 100, 0., (stopTime - startTime) / readOutFreq, 200, 1.5, mw1_focusxHist->GetMaximum() * 1.05);
timeIntSignal_MWPC = new TH1D("timeIntSignal_MWPC", "time integrated signal", 64, 0.5, 64.5);
for (int j = 0; j < 20; j++)
{
TString name;
name.Form("sigIcScatterPlot%d", j);
sigIcScatterPlot[j] = new TH2D(name.Data(), name.Data(), 200, 0, -1, 200, 0, -1);
TString name2;
name2.Form("mwPosScatterPlot%d", j);
mwPosScatterPlot[j] = new TH2D(name2.Data(), name2.Data(), 200, mean_MWPC - 0.1 * sigma_ave, mean_MWPC + 0.1 * sigma_ave, 200, 0.8 * (mean_ave - 0.15 * sigma_ave), 0.8 * (mean_ave + 0.15 * sigma_ave));
}
sigIcScatterPlot_ic2 = new TH2D("sigIcScatterPlot_ic2", "sigIcScatterPlot_ic2", 200, 0, -1, 200, 0, -1);
Center_Signal = new TH2D("Center_Signal", "", 200, 0., max_ic1 - 5000, 200, (0.8 * mean_ave - 5 * sigma_pos), (0.8 * mean_ave + 5 * sigma_pos));
Diff_timeCenters = new TH2D("Diff_time", "", nbins / 100, 0., (stopTime - startTime) / readOutFreq, 200, -6 * sigma_pos, 6 * sigma_pos);
Diff_hist = new TH1D("Diff_hist", "l", 200, -2, 2.);
Diff_Signal = new TH2D("Diff_Signal", "", 200, 0., max_ic1 - 5000, 200, -6 * sigma_pos, 6 * sigma_pos);
Diff_ic1 = new TH2D("Diff_ic1", "", 200, 0., max_ic1 - 5000, 200, -6 * sigma_pos, 6 * sigma_pos);
b_beamPosX_1 = newdata->Branch("beamPosX_1", &beamPosX_1);
b_beamFocusX_1 = newdata->Branch("beamFocusX_1", &beamFocusX_1);
b_beamSignal_1 = newdata->Branch("beamSignal_1", &beamSignal_1);
bb_ic1 = newdata->Branch("ic1_1", &ic1_1);
bb_ic2 = newdata->Branch("ic2_1", &ic2_1);
bb_mw1_focusx = newdata->Branch("mw1_focusx_1", &mw1_focusx_1);
bb_mw1_posx = newdata->Branch("mw1_posx_1", &mw1_posx_1);
bb_mw2_focusx = newdata->Branch("mw2_focusx_1", &mw2_focusx_1);
bb_mw2_posx = newdata->Branch("mw2_posx_1", &mw2_posx_1);
b_beamOn = newdata->Branch("beamOn", &beamOn);
newdata->SetBranchStatus("*", 1);
for (int i = 0; i < spillTimes.size(); i++)
{
tmp_TimeCenter = calcTimeCenters(spillTimes[i].first + 200, spillTimes[i].second);
timeIntSignal->Add(tmp_TimeCenter, 1.);
delete tmp_TimeCenter;
beamPosX_1 = 0;
beamFocusX_1 = 0;
beamSignal_1 = 0;
ic1_1 = 0;
mw1_focusx_1 = 0;
mw1_posx_1 = 0;
beamOn = 0;
newdata->Fill();
}
SignalErrHist->Write();
printf("Fit_Fail_Counts=%d\n", Fit_Fail_Counts);
//Draw timeCenter
timeCenters->GetXaxis()->SetTitle("time [s]");
TString label;
label.Form("mean position[mm]");
timeCenters->GetYaxis()->SetTitle(label.Data());
timeCenters->Write();
//Draw timeCenter_MWPC
timeCenters_MWPC->GetXaxis()->SetTitle("time [s]");
TString label2;
label2.Form("mean position[mm]");
timeCenters_MWPC->GetYaxis()->SetTitle(label.Data());
timeCenters_MWPC->Write();
//Draw timeWidths
timeWidths->GetXaxis()->SetTitle("time [s]");
TString label12;
label12.Form("width fiber[mm]");
timeWidths->GetYaxis()->SetTitle(label12.Data());
timeWidths->Write();
//Draw timeWidths_MWPC
timeWidths_MWPC->GetXaxis()->SetTitle("time [s]");
TString label13;
label13.Form("width MWPC[mm]");
timeWidths_MWPC->GetYaxis()->SetTitle(label12.Data());
timeWidths_MWPC->Write();
Center_Signal->GetXaxis()->SetTitle("Signal [particles]");
TString label123;
label123.Form("Position_Fiber[mm]");
Center_Signal->GetYaxis()->SetTitle(label123.Data());
Center_Signal->Write();
Center_Signal->ProfileX()->Write();
//Draw Diff_timeCenters
Diff_timeCenters->GetXaxis()->SetTitle("time [s]");
TString label3;
label3.Form("mean position diff[mm]");
Diff_timeCenters->GetYaxis()->SetTitle(label.Data());
Diff_timeCenters->Write();
//Draw Diff_Signal
Diff_Signal->GetXaxis()->SetTitle("signal");
TString label4;
label4.Form("mean position diff[mm]");
Diff_Signal->GetYaxis()->SetTitle(label.Data());
Diff_Signal->Write();
//Draw Diff_Signal
Diff_ic1->GetXaxis()->SetTitle("ic1");
TString label5;
label5.Form("mean position diff[mm]");
Diff_ic1->GetYaxis()->SetTitle(label.Data());
Diff_ic1->Write();
TObjArray aSlices, bSlices;
Diff_Signal->FitSlicesY(0, 0, -1, 5, "QNR", &aSlices);
for (int i = 0; i < 4; i++)
{
aSlices[i]->Write();
}
Diff_ic1->FitSlicesY(0, 0, -1, 5, "QNR", &bSlices);
for (int i = 0; i < 4; i++)
{
bSlices[i]->Write();
}
mean_ave = beamPosX->GetMean();
sigma_pos = beamPosX->GetStdDev();
sigma_ave = beamFocusX->GetMean();
return;
}
pair<int, int> Beammon::GetSpillTimes(int nSpill)
{
return spillTimes[nSpill];
}
TH2D *Beammon::GetSpillCenters(int nSpill)
{
//return calcTimeCenters(spillTimes[nSpill].first, spillTimes[nSpill].second);
}
TH1D *Beammon::GetSpillSlice(int nSpill)
{
return GetTimeSlice(-1, spillTimes[nSpill].first - 500, spillTimes[nSpill].second + 500, true);
}
TH1D *Beammon::GetSignalDistribution()
{
/*int n = signalDist->GetNbinsX();
for(int i = 1; i<= n; i++){
signalDist->SetBinContent(i,signalDist->GetBinContent(i)*signalDistScale);
}*/
signalDist->Scale(signalDistScale);
return signalDist;
}
TCanvas *Beammon::GetSignalDistributionCanvas()
{
TCanvas *c1 = new TCanvas();
c1->SetName("signalDistributionCanvas");
signalDist->Draw();
signalDist->SetTitle("");
signalDist->GetYaxis()->SetTitle("entries");
TString label;
label.Form("integrated signal [adc counts / %d #mus]", intTime);
signalDist->GetXaxis()->SetTitle(label.Data());
//c1->SaveAs("dummy");
return c1;
}
TH2D *Beammon::GetSimpleView()
{
return timePosHist;
}
TH1D *Beammon::GetBaselineHist()
{
TH1D *tmp = new TH1D("baseline", "baseline", 64, 0.5, 64.5);
for (int n = 1; n <= 64; n++)
{
tmp->SetBinContent(n, baseline[n - 1]);
}
return tmp;
}
void Beammon::Close()
{
delete timePosHist;
fileOut->Close();
}
double Beammon::GetIntSignal()
{
return intSignal;
}
int Beammon::GetBeamOnTime()
{
return beamOnTime;
}
double Beammon::GetSignalStd()
{
return signalDist->GetStdDev(1);
}
bool Beammon::isAboveThreshold(TGraphErrors *slice)
{
double channel, signal;
slice->GetPoint(highestChannelBin, channel, signal);
if (signal > beamOnAdcVal)
return true;
else
return false;
}
bool Beammon::isAboveThreshold(TH1D *slice)
{
double val = slice->GetBinContent(highestChannelBin);
if (val > beamOnAdcVal)
return true;
else
return false;
}
bool Beammon::isAboveThreshold(double val)
{
// old version using the gaussian fit to the measurement
if (val > beamOnAdcVal)
return true;
else
return false;
/*
if(val > baselineIntegral * 1.5 && val < baselineIntegral *100)
return true;
else
return false;*/
}
TH1D *Beammon::icSignalComp()
{
icSignalDist = new TH1D("icSignal", "icSignal", 200, 0., ic1Hist->GetMaximum());
TH1D *icSignalDist2 = new TH1D("icSignal2", "icSignal2", 200, 0., ic2Hist->GetMaximum());
for (int i = 0; i < ic1Hist->GetNbinsX(); i++)
{
if (ic1Hist->GetBinContent(i) != 0)
icSignalDist->Fill(ic1Hist->GetBinContent(i));
if (ic2Hist->GetBinContent(i) != 0)
icSignalDist2->Fill(ic2Hist->GetBinContent(i));
}
icSignalDist2->Write();
return icSignalDist;
}
TCanvas *Beammon::signalOverTime2d()
{
TCanvas *c1 = new TCanvas();
c1->SetName("signalOverTime2d");
TH1D *tmp = GetSignalOverTime();
int binsX = tmp->GetNbinsX();
int max = tmp->GetMaximum();
int min = tmp->GetMinimum();
TH2D *res = new TH2D("sigTime2d", "", binsX / 200, 0., binsX / readOutFreq, 100, min, max);
for (int k = 0; k < binsX; k++)
{
res->Fill(k / readOutFreq, tmp->GetBinContent(k));
}
res->Draw("colz");
res->GetXaxis()->SetTitle("time [s]");
TString label;
label.Form("adc counts / %d #mus", intTime);
res->GetYaxis()->SetTitle(label.Data());
//c1->SaveAs("dummy");
res->Write();
return c1;
}
double Beammon::calcBeamTime()
{
beamOnTime = 0;
for (int i = 0; i < spillTimes.size(); i++)
{
beamOnTime += spillTimes[i].second - spillTimes[i].first;
}
return beamOnTime;
}
TCanvas *Beammon::PrintProfile(int offset)
{
TCanvas *c1 = new TCanvas();
c1->cd();
//TH1D *res = GetChannelSlice(-1,offset);
TGraphErrors *res = new TGraphErrors(64);
TGraphErrors *tmp = GetChannelGraph(-1, offset);
double x, y, errx, erry;
for (int ch = 0; ch < 63; ch++)
{
tmp->GetPoint(ch, x, y);
x *= 0.8;
errx = 0;
erry = tmp->GetErrorY(ch);
res->SetPoint(ch, x, y);
res->SetPointError(ch, errx, erry);
}
res->Fit("gaus", "Q0");
double stepsize = 0.01;
double lim = res->GetFunction("gaus")->GetMaximum() / 2;
double val, x1, x2, pos = 0;
bool below = true;
while (pos <= 64)
{
val = res->GetFunction("gaus")->Eval(pos);
if (below && val >= lim)
{
x1 = pos;
below = false;
}
else if (below == false && val <= lim)
{
x2 = pos;
break;
}
pos += stepsize;
}
fwhm = (x2 - x1) * channelWidth;
res->Draw("AP");
TString sName;
sName.Form("Profile%d", offset);
c1->SetName(sName.Data());
res->SetTitle("");
res->GetXaxis()->SetTitle("position [mm]");
TString label;
label.Form("adc counts / %d #mus", intTime);
res->GetYaxis()->SetTitle(label.Data());
res->GetYaxis()->SetTitleOffset(1.4);
TString str;
str.Form("fwhm = %.3f mm", fwhm);
TText *t = new TText(34, 900, str.Data());
t->Draw("SAME");
delete t;
delete res;
c1->Update();
c1->Modified();
//c1->SaveAs(sName.Data());
return c1;
}
TCanvas *Beammon::FitSnapshot(int offset)
{
//gSyle->SetOptFit(1111);
TCanvas *c1 = new TCanvas();
c1->cd();
TGraphErrors *res = GetChannelGraph(-1, offset);
TString sName;
sName.Form("Snapshot%d", offset);
c1->SetName(sName.Data());
res->SetTitle("");
res->GetXaxis()->SetTitle("channel");
TString label;
label.Form("adc counts / %d #mus", intTime);
res->GetYaxis()->SetTitle(label.Data());
res->GetYaxis()->SetTitleOffset(1.4);
//res->Fit("gaus");
res->Fit("gaus");
double mean = res->GetFunction("gaus")->GetParameter(1);
double sigma = res->GetFunction("gaus")->GetParameter(2);
double constant = res->GetFunction("gaus")->GetParameter(0);
TF1 *func = new TF1("doublegaus", fitfunc, 0, 64, 6); // t
func->SetNpx(1000);
func->SetParameter(0, constant);
func->SetParameter(1, mean);
func->SetParameter(2, sigma);
func->SetParameter(3, constant);
func->SetParameter(4, sigma);
func->SetParameter(5, 0.);
func->SetParLimits(1, 0., 64.);
func->SetParLimits(2, 1., 20.);
func->SetParLimits(4, 5., 100.);
//func->SetParameter(5,0.);
res->Fit("doublegaus");
res->Draw("AP");
c1->SaveAs(sName.Data());
return c1;
}
TCanvas *Beammon::FitAverage()
{
TCanvas *c1 = new TCanvas();
c1->cd();
int n = spillTimes.size();
double yValues[64];
double yErrors[64];
double xValues[64];
double xErrors[64];
int count = 0;
double x, y, errY;
TGraphErrors *tmp;
yValues[0] = 0;
for (int i = 0; i < n; i++)
{
int start = spillTimes[i].first;
int stop = spillTimes[i].second;
for (int k = start; k <= stop; k++)
{
count++;
tmp = GetChannelGraph(k, 0);
for (int bin = 1; bin <= 64; bin++)
{
errY = tmp->GetErrorY(bin);
tmp->GetPoint(bin, x, y);
if (!(y != y || errY != errY))
{
yValues[bin - 1] += y;
yErrors[bin - 1] += errY * errY;
}
else
{
yValues[bin - 1] += 0;
yErrors[bin - 1] += 0;
}
}
delete tmp;
}
}
for (int bin = 0; bin < 64; bin++)
{
xValues[bin] = bin;
xErrors[bin] = 1. / sqrt(12);
yValues[bin] *= 1. / count;
yErrors[bin] = sqrt(yErrors[bin]) / count;
}
/*for (int i = 0; i < 64; i++)
{
printf("channelErrors[%d]: %f ", i, channelErrors[i]);
printf("yValue[%d]: %f ", i, yValues[i]);
printf("yError[%d]: %f \n", i, yErrors[i]);
}*/
TGraph *res = new TGraph(64, xValues, yValues); //, xErrors, yErrors);
c1->SetName("avgGraph");
res->SetTitle("");
res->GetXaxis()->SetTitle("channel");
TString label;
label.Form("adc counts / %d #mus", intTime);
res->GetYaxis()->SetTitle(label.Data());
res->GetYaxis()->SetTitleOffset(1.4);
res->Fit("gaus", "Q0");
double mean = res->GetFunction("gaus")->GetParameter(1);
//mean_ave = mean;
double sigma = res->GetFunction("gaus")->GetParameter(2);
//sigma_ave = sigma;
double constant = res->GetFunction("gaus")->GetParameter(0);
//res->Fit("gaus","","",mean-1.5*sigma,mean+1.5*sigma);
TF1 *func = new TF1("doublegaus", fitfunc, 0, 64, 6); // t
func->SetNpx(1000);
func->SetParameter(0, constant);
func->SetParameter(1, mean);
func->SetParameter(2, sigma);
func->SetParameter(3, constant);
func->SetParameter(4, sigma);
func->SetParameter(5, 0.);
func->SetParLimits(1, 0., 64.);
func->SetParLimits(2, 1., 20.);
func->SetParLimits(4, 5., 100.);
//func->SetParameter(5,0.);
res->Fit("doublegaus");
res->Draw("AP");
//c1->SaveAs("avgFit");
// find fwhm
double stepsize = 0.01;
double lim = res->GetFunction("doublegaus")->GetMaximum() / 2;
double val, x1, x2, pos = 0;
bool below = true;
while (pos <= 64)
{
val = func->Eval(pos);
if (below && val >= lim)
{
x1 = pos;
below = false;
}
else if (below == false && val <= lim)
{
x2 = pos;
break;
}
pos += stepsize;
}
fwhm = (x2 - x1) * channelWidth;
return c1;
}
TCanvas *Beammon::PrintChannelNoise(int channel)
{
TCanvas *c1 = new TCanvas();
c1->cd();
TString name;
name.Form("channelNoise%d", channel);
TString label;
label.Form("adc counts / %d #mus", intTime);
c1->SetName(name.Data());
channelSig[channel]->SetTitle("");
channelSig[channel]->SetName(name.Data());
channelSig[channel]->GetXaxis()->SetTitle(label.Data());
channelSig[channel]->GetYaxis()->SetTitle("entries");
channelSig[channel]->Draw();
gStyle->SetOptStat("mr");
gPad->Modified();
//c1->SaveAs("blub");
return c1;
}
TCanvas *Beammon::IcFitSlicesY()
{
TCanvas *c1 = new TCanvas();
c1->cd();
c1->SetName("icFitSlicesY");
TH1D *tmp = timePosHist->ProjectionY("", 1, 64);
TH1D *means = (TH1D *)sigIcScatterPlot[10 + alignintIc]->ProfileY(); //(TH1D*)objectArray[1];
means->SetLineColor(kBlack);
//sigIcScatterPlot->Draw();
means->Draw("E");
means->GetXaxis()->SetRangeUser(0., means->GetXaxis()->GetXmax());
means->SetMinimum(0.);
means->SetTitle("");
//means->SetName("icFitSlicesY");
TString label1, label2;
label1.Form("mean photodiode signal [particles / %d #mus]", intTime);
label2.Form("ionisation chamber signal [particles / %d #mus]", intTime);
means->GetXaxis()->SetTitle(label2.Data());
means->GetYaxis()->SetTitle(label1.Data());
means->Fit("pol1");
//sigIcScatterPlot->Draw();
gPad->Modified();
//c1->SaveAs("blub");
means->Write();
return c1;
}
TCanvas *Beammon::PrintSigIcGraph()
{
TCanvas *c1 = new TCanvas();
c1->cd();
TString name;
name.Form("sigIcScatterGraph");
TString label, label2;
label.Form("adc counts / %d #mus", intTime);
label2.Form("whatever / %d #mus", intTime);
c1->SetName(name.Data());
sigIcScatterGraph->SetTitle("");
sigIcScatterGraph->SetName(name.Data());
sigIcScatterGraph->GetXaxis()->SetTitle(label2.Data());
sigIcScatterGraph->GetYaxis()->SetTitle(label.Data());
sigIcScatterGraph->Fit("pol1", "W");
sigIcScatterGraph->Draw("AP");
gPad->Modified();
//c1->SaveAs("blub");
return c1;
}
TH1D *Beammon::GetMwPosComp()
{
double bin, mean;
double shift = beamPosX->GetMean() - mwPosX->GetMean();
//printf("beam: %f \n",beamPosTimeGraph->GetMean(2));
int n = beamPosTimeGraph->GetN();
for (int i = 0; i < n; i++)
{
beamPosTimeGraph->GetPoint(i, bin, mean);
if (!mean)
continue;
positionComp->Fill(mean - mw1_posxHist->GetBinContent(bin - alignint) - shift);
}
return positionComp;
}
//
///
/////////////////////////////////////////////////////////////////////////////
TCanvas *Beammon::GetSig_e_ic_Hist()
{
shift = beamPosX->GetMean() - mwPosX->GetMean();
printf("shift = %f", shift);
TCanvas *c1 = new TCanvas();
c1->cd();
c1->SetName("SignalComp_vs_Time");
Double_t norm = ic1Hist->Integral();
ic1Hist->Scale(1. / norm);
ic1Hist->SetLineColor(kRed);
ic1Hist->Draw("HIST L");
Double_t norm0 = sigHist->Integral();
Double_t scaleSig = norm / norm0;
sigHist->Scale(1. / norm0);
sigHist->Draw("HIST SAME L");
TLegend * legend= new TLegend(0.1, 0.7, 0.48, 0.9);
legend->AddEntry(ic1Hist, "ic", "l");
legend->AddEntry(sigHist, "Fibre", "l");
legend->Draw();
TH1D *compsig = new TH1D("compsig", "compsig", 200, -0.5, 0.5);
for (int i = 0; i < nevents; i++)
{
compsig->Fill((ic1Hist->GetBinContent(i) - sigHist->GetBinContent(i)) * 2 / (ic1Hist->GetBinContent(i) + sigHist->GetBinContent(i)));
}
compsig->Write();
TCanvas *c2 = new TCanvas();
c2->cd();
c2->SetName("PositionComp_vs_Time");
TH1D *fibreHistNew = new TH1D("fibreHistNew", "fibreHistNew", nevents, 0., nevents - 1);
for (int i = 0; i < nevents; i++)
{
if (fibreHist->GetBinContent(i) != 0)
{
fibreHistNew->SetBinContent(i, fibreHist->GetBinContent(i) - shift);
}
else
fibreHistNew->SetBinContent(i, 0);
}
fibreHistNew->Draw("HIST L");
mwpcHist->SetLineColor(kRed);
mwpcHist->SetLineColor(kRed);
mwpcHist->Draw("HIST SAME L");
TLegend * legend2 = new TLegend(0.1, 0.7, 0.48, 0.9);
legend2->AddEntry(mwpcHist, "MWPC", "l");
legend2->AddEntry(fibreHistNew, "Fibre", "l");
legend2->Draw();
c2->Write();
return c1;
}
TCanvas *Beammon::GetSig_e_ic_Hist_Weighted()
{
TCanvas *c1 = new TCanvas();
c1->cd();
c1->SetName("SignalComp_vs_Time_all");
Double_t norm = ic1Hist->Integral();
ic1Hist->Scale(1. / norm);
ic1Hist->SetLineColor(kRed);
ic1Hist->SetLineWidth(2);
ic1Hist->Draw("HIST L");
Double_t norm0 = sigHist_weighted->Integral();
sigHist_weighted->Scale(1. / norm0);
sigHist_weighted->Draw("HIST SAME L");
Double_t norm1 = sigHist->Integral();
sigHist->Scale(1. / norm1);
sigHist->SetLineColor(kBlue);
sigHist->Draw("HIST SAME L");
TLegend * legend= new TLegend(0.1, 0.7, 0.48, 0.9);
legend->AddEntry(ic1Hist, "ic", "l");
legend->AddEntry(sigHist, "Fibre_fit", "l");
legend->AddEntry(sigHist_weighted, "Fibre_sum", "l");
legend->Draw();
c1->Write();
TH1D *compsig_weighted = new TH1D("compsig_weighted", "compsig", 200, -0.5, 0.5);
for (int i = 0; i < nevents; i++)
{
compsig_weighted->Fill((ic1Hist->GetBinContent(i) - sigHist_weighted->GetBinContent(i)) * 2 / (ic1Hist->GetBinContent(i) + sigHist_weighted->GetBinContent(i)));
}
compsig_weighted->SetLineColor(kRed);
compsig_weighted->Write();
TH1D *compsig = new TH1D("compsig", "compsig", 200, -0.5, 0.5);
for (int i = 0; i < nevents; i++)
{
compsig->Fill((ic1Hist->GetBinContent(i) - sigHist->GetBinContent(i)) * 2 / (ic1Hist->GetBinContent(i) + sigHist->GetBinContent(i)));
}
c1->cd();
c1->SetName("CompSig_all");
compsig->Draw();
compsig_weighted->Draw("SAME");
TLegend * legend2 = new TLegend(0.1, 0.7, 0.48, 0.9);
legend2->AddEntry(compsig_weighted, "Weighted", "l");
legend2->AddEntry(compsig, "Fit", "l");
legend2->Draw();
c1->Write();
TCanvas *c2 = new TCanvas();
c2->cd();
c2->SetName("PositionComp_vs_Time_all");
TH1D *fibreHistNew_weighted = new TH1D("fibreHistNew_weighted", "fibreHistNew_weighted", nevents, 0., nevents - 1);
for (int i = 0; i < nevents; i++)
{
if (fibreHist_weighted->GetBinContent(i) != 0)
{
fibreHistNew_weighted->SetBinContent(i, fibreHist_weighted->GetBinContent(i) - shift);
}
else
fibreHistNew_weighted->SetBinContent(i, 0);
}
TH1D *fibreHistNew = new TH1D("fibreHistNew", "fibreHistNew", nevents, 0., nevents - 1);
for (int i = 0; i < nevents; i++)
{
if (fibreHist->GetBinContent(i) != 0)
{
fibreHistNew->SetBinContent(i, fibreHist->GetBinContent(i) - shift);
}
else
fibreHistNew->SetBinContent(i, 0);
}
fibreHistNew->Draw("HIST L");
fibreHistNew->SetLineColor(kBlue);
fibreHistNew_weighted->Draw("HIST L SAME");
mwpcHist->SetLineColor(kRed);
mwpcHist->SetLineWidth(2);
mwpcHist->Draw("HIST SAME L");
TLegend * legend3 = new TLegend(0.1, 0.7, 0.48, 0.9);
legend3->AddEntry(mwpcHist, "MWPC", "l");
legend3->AddEntry(fibreHistNew, "Fibre_fit", "l");
legend3->AddEntry(fibreHistNew_weighted, "Fibre_Weighted", "l");
legend3->Draw();
c2->Write();
c2->cd();
c2->Clear();
c2->SetName("PosComp_All");
TH1D *comppos_weighted = new TH1D("comppos_weighted", "comppos", 200, -0.4, 0.4);
for (int i = 0; i < nevents; i++)
{
comppos_weighted->Fill((mwpcHist->GetBinContent(i) - fibreHist_weighted->GetBinContent(i) + shift));
}
comppos_weighted->SetLineColor(kRed);
comppos_weighted->Write();
TH1D *comppos = new TH1D("comppos", "comppos", 200, -0.4, 0.4);
for (int i = 0; i < nevents; i++)
{
comppos->Fill((mwpcHist->GetBinContent(i) - fibreHist->GetBinContent(i) + shift));
}
comppos->Draw();
comppos_weighted->Draw("SAME");
TLegend * legend4 = new TLegend(0.1, 0.7, 0.48, 0.9);
legend4->AddEntry(comppos_weighted, "Weighted", "l");
legend4->AddEntry(comppos, "Fit", "l");
legend4->Draw();
return c2;
}
Reference in New Issue View Git Blame Copy Permalink