Added the size check to the range functions and improved plotting.

This commit is contained in:
redford_s
2019-02-08 11:40:43 +01:00
parent 0dd0b92d24
commit f3ba85d343

View File

@ -13,8 +13,11 @@
double highestPointBeforeSwitching(const vector<double> &lower_filter, const vector<double> &higher_filter) {
double highest_point_before_switching;
if (higher_filter.size() > 0) {
// find the highest value in lower_filter that is lower than all entries in higher_filter
double highest_point_before_switching = *min_element(lower_filter.begin(),lower_filter.end())-1.;
highest_point_before_switching = *min_element(lower_filter.begin(),lower_filter.end())-1;
double lowest_entry_in_higher_filter = *min_element(higher_filter.begin(),higher_filter.end());
for(vector<double>::const_iterator it = lower_filter.begin(); it != lower_filter.end(); ++it) {
@ -24,6 +27,9 @@ double highestPointBeforeSwitching(const vector<double> &lower_filter, const vec
}
}
}
} else {
highest_point_before_switching = *max_element(lower_filter.begin(),lower_filter.end());
}
return highest_point_before_switching;
@ -31,8 +37,11 @@ double highestPointBeforeSwitching(const vector<double> &lower_filter, const vec
double lowestPointAfterSwitching(const vector<double> &higher_filter, const vector<double> &lower_filter) {
double lowest_point_after_switching;
if (lower_filter.size() > 0) {
// find the lowest value in higher_filter that is higher than all entries in lower_filter
double lowest_point_after_switching = *max_element(higher_filter.begin(),higher_filter.end())+1;
lowest_point_after_switching = *max_element(higher_filter.begin(),higher_filter.end())+1;
double highest_entry_in_lower_filter = *max_element(lower_filter.begin(),lower_filter.end());
for(vector<double>::const_iterator it = higher_filter.begin(); it != higher_filter.end(); ++it) {
@ -42,6 +51,9 @@ double lowestPointAfterSwitching(const vector<double> &higher_filter, const vect
}
}
}
} else {
lowest_point_after_switching = *min_element(higher_filter.begin(),higher_filter.end());
}
return lowest_point_after_switching;
@ -452,74 +464,36 @@ int main(int argc, char* argv[]) {
r1_adc.push_back(this_g1avg);
r1_adcer.push_back(this_g1avger);
}
}
if (r0_adc.size() > 1 && r1_adc.size() > 1) {
TGraphErrors *grap_g0 = 0;
TGraphErrors *grap_g1 = 0;
double rangemin0 = *min_element(r0_filter.begin(),r0_filter.end());
double rangemax0 = highestPointBeforeSwitching(r0_filter,r1_filter);
range0maxhist->Fill(rangemax0);
range0maxmap->Fill(i%NC,i/NC,rangemax0);
double rangemin1 = lowestPointAfterSwitching(r1_filter,r0_filter);
range1minhist->Fill(rangemin1);
range1minmap->Fill(i%NC,i/NC,rangemin1);
double rangemax1 = *max_element(r1_filter.begin(),r1_filter.end());
if (rangemax0 > rangemin0 && rangemax1 > rangemin1) {
TGraphErrors *grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_filterer[0]),&(r0_adcer[0]));
TGraphErrors *grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_filterer[0]),&(r1_adcer[0]));
TF1 *fit1 = new TF1("fit1","[0]+[1]*x",rangemin0,rangemax0);
TF1 *fit1 = 0;
TF1 *fit2 = 0;
fit1->SetParameter(0, 100.);
fit1->SetParameter(1, 1);
fit1->SetParName(0,"G0 const");
fit1->SetParName(1,"G0 grad");
fit1->SetLineColor(kBlue);
grap_g0->Fit(fit1,"QR","");
double rangemin0 = 0;
double rangemax0 = 0;
double rangemin1 = 0;
double rangemax1 = 0;
g0map->Fill(i%NC,i/NC,fit1->GetParameter(1));
g0ermap->Fill(i%NC,i/NC,fit1->GetParError(1));
// define graphs
if (r0_adc.size() > 1) {
grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_filterer[0]),&(r0_adcer[0]));
grap_g0->SetMarkerStyle(20);
grap_g0->SetMarkerColor(kBlue);
grap_g0->SetLineColor(kBlue);
}
fit2 = new TF1("fit2","[0]+[1]*x",rangemin1,rangemax1);
if (r1_adc.size() > 1) {
grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_filterer[0]),&(r1_adcer[0]));
grap_g1->SetMarkerStyle(20);
grap_g1->SetMarkerColor(kGreen+2);
grap_g1->SetLineColor(kGreen+2);
}
fit2->SetParameter(0, 1000.);
fit2->SetParameter(1, -0.1);
fit2->SetParName(0,"G1 const");
fit2->SetParName(1,"G1 grad");
fit2->SetLineColor(kGreen+2);
grap_g1->Fit(fit2,"QR","");
g1map->Fill(i%NC,i/NC,fit2->GetParameter(1));
g1ermap->Fill(i%NC,i/NC,fit2->GetParError(1));
g0overg1map->Fill(i%NC,i/NC,fit1->GetParameter(1)/fit2->GetParameter(1));
g0overg1hist->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
g0overg1ermap->Fill(i%NC,i/NC,abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));
g0overg1erhist->Fill(abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));
if (isEdge(i)) {
g0overg1hist_isEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isInnerEdge(i)) {
g0overg1hist_isInnerEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isDouble(i)) {
g0overg1hist_isDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isNextToDouble(i)) {
g0overg1hist_isNextToDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isQuad(i)) {
g0overg1hist_isQuad->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isBulk(i)) {
g0overg1hist_isBulk->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
// plot the datapoints
if (r1_adc.size() > 1) {
if ((i >= 58000 && i < 58000+10) || // bulk
(i >= 10 && i < 10+10) || // edge
@ -547,25 +521,80 @@ int main(int argc, char* argv[]) {
mapcanvas->SetLeftMargin(0.13);
mapcanvas->SetRightMargin(0.05);
grap_g0->SetMarkerStyle(20);
grap_g0->SetMarkerColor(kBlue);
grap_g0->SetLineColor(kBlue);
grap_g1->SetMarkerStyle(20);
grap_g1->SetMarkerColor(kGreen+2);
grap_g1->SetLineColor(kGreen+2);
grap_g1->SetMinimum(0);
grap_g1->SetMaximum(17000);
grap_g1->GetXaxis()->SetTitle("Tube current [mA]");
grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
grap_g1->GetYaxis()->SetTitleOffset(0.9);
grap_g1->GetXaxis()->SetLimits(0,85);
grap_g1->SetMinimum(0);
grap_g1->SetMaximum(17000);
grap_g1->Draw("AP");
if (r0_adc.size() > 1) {
grap_g0->Draw("P");
fit1->Draw("same");
fit2->Draw("same");
mapcanvas->Update();
}
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
mapcanvas->SetLeftMargin(0.1);
mapcanvas->SetRightMargin(0.13);
}
}
// define fit ranges and fit
if (r0_adc.size() > 1) {
rangemin0 = *min_element(r0_filter.begin(),r0_filter.end());
rangemax0 = highestPointBeforeSwitching(r0_filter,r1_filter);
range0maxhist->Fill(rangemax0);
range0maxmap->Fill(i%NC,i/NC,rangemax0);
if (rangemax0 > rangemin0) {
fit1 = new TF1("fit1","[0]+[1]*x",rangemin0,rangemax0);
fit1->SetParameter(0, 100.);
fit1->SetParameter(1, 1);
fit1->SetLineColor(kBlue);
fit1->SetParName(0,"G0 const");
fit1->SetParName(1,"G0 grad");
grap_g0->Fit(fit1,"QR","");
g0map->Fill(i%NC,i/NC,fit1->GetParameter(1));
g0ermap->Fill(i%NC,i/NC,fit1->GetParError(1));
if ((i >= 58000 && i < 58000+10) || // bulk
(i >= 10 && i < 10+10) || // edge
(i >= 1024+10 && i < 1024+10+10) || // inner edge
(i >= (256*1024)+10 && i < (256*1024)+10+10) || // double
(i >= (257*1024)+10 && i < (257*1024)+10+10) || // next to double
(i == (255*1024)+255) // quad
) {
string pixel_type = "x";
if (i >= 58000 && i < 58000+10) {
pixel_type = "b";
} else if (i >= 10 && i < 10+10) {
pixel_type = "e";
} else if (i >= 1024+10 && i < 1024+10+10) {
pixel_type = "ie";
} else if (i >= (256*1024)+10 && i < (256*1024)+10+10) {
pixel_type = "d";
} else if (i >= (257*1024)+10 && i < (257*1024)+10+10) {
pixel_type = "ntd";
} else if (i == (255*1024)+255) {
pixel_type = "q";
}
mapcanvas->SetLeftMargin(0.13);
mapcanvas->SetRightMargin(0.05);
grap_g0->SetMinimum(0);
grap_g0->SetMaximum(17000);
grap_g0->GetXaxis()->SetTitle("Tube current [mA]");
grap_g0->GetYaxis()->SetTitle("ADC [ADU]");
grap_g0->GetYaxis()->SetTitleOffset(0.9);
grap_g0->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
grap_g0->Draw("AP");
fit1->Draw("same");
mapcanvas->Update();
TPaveStats *st0 = (TPaveStats*)grap_g0->FindObject("stats");
st0->SetX1NDC(0.2);
st0->SetX2NDC(0.54);
@ -573,32 +602,20 @@ int main(int argc, char* argv[]) {
st0->SetY2NDC(0.37);
st0->SetBorderSize(0);
st0->SetTextColor(kBlue);
TPaveStats *st1 = (TPaveStats*)grap_g1->FindObject("stats");
st1->SetX1NDC(0.6);
st1->SetX2NDC(0.94);
st1->SetY1NDC(0.18);
st1->SetY2NDC(0.37);
st1->SetBorderSize(0);
st1->SetTextColor(kGreen+2);
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
// normalise
vector<double> r0_adc_norm;
vector<double> r1_adc_norm;
for (size_t j = 0; j < r0_adc.size(); j++) {
r0_adc_norm.push_back(r0_adc[j] - fit1->Eval(r0_filter[j]));
}
for (size_t j = 0; j < r1_adc.size(); j++) {
r1_adc_norm.push_back(r1_adc[j] - fit2->Eval(r1_filter[j]));
}
TGraphErrors *grap0_norm = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc_norm[0]),&(r0_filterer[0]),&(r0_adcer[0]));
TGraphErrors *grap1_norm = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_filterer[0]),&(r1_adcer[0]));
grap0_norm->SetMarkerColor(kBlue);
grap0_norm->SetLineColor(kBlue);
TF1* flat_g0 = new TF1("flat_g0","0",rangemin0, rangemax0);
TF1* flat_g1 = new TF1("flat_g1","0",rangemin1, rangemax1);
flat_g0->SetLineColor(kBlue);
TF1* lin_g0_p02pc = new TF1("lin_g0_p02pc","[0]+[1]*x",rangemin0,rangemax0);
lin_g0_p02pc->SetParameter(0,fit1->GetParameter(0)/500.);
@ -620,6 +637,109 @@ int main(int argc, char* argv[]) {
lin_g0_m05pc->SetParameter(1,fit1->GetParameter(1)/-200.);
lin_g0_m05pc->SetLineColor(kOrange+1);
grap0_norm->GetXaxis()->SetTitle("Tube current [mA]");
grap0_norm->GetYaxis()->SetTitle("Normalised ADC [ADU]");
grap0_norm->GetYaxis()->SetTitleOffset(0.9);
grap0_norm->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
grap0_norm->Draw("AP");
flat_g0->Draw("same");
lin_g0_p02pc->Draw("same");
lin_g0_p05pc->Draw("same");
lin_g0_m02pc->Draw("same");
lin_g0_m05pc->Draw("same");
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
grap0_norm->SetMinimum(-100);
grap0_norm->SetMaximum(100);
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_zoom_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
delete grap0_norm;
delete flat_g0;
mapcanvas->SetLeftMargin(0.1);
mapcanvas->SetRightMargin(0.13);
}
}
}
if (r1_adc.size() > 1) {
rangemin1 = lowestPointAfterSwitching(r1_filter,r0_filter);
range1minhist->Fill(rangemin1);
range1minmap->Fill(i%NC,i/NC,rangemin1);
rangemax1 = *max_element(r1_filter.begin(),r1_filter.end());
if (rangemax1 > rangemin1) {
fit2 = new TF1("fit2","[0]+[1]*x",rangemin1,rangemax1);
fit2->SetParameter(0, 1000.);
fit2->SetParameter(1, -0.1);
fit2->SetLineColor(kGreen+2);
fit2->SetParName(0,"G1 const");
fit2->SetParName(1,"G1 grad");
grap_g1->Fit(fit2,"QR","");
g1map->Fill(i%NC,i/NC,fit2->GetParameter(1));
g1ermap->Fill(i%NC,i/NC,fit2->GetParError(1));
if ((i >= 58000 && i < 58000+10) || // bulk
(i >= 10 && i < 10+10) || // edge
(i >= 1024+10 && i < 1024+10+10) || // inner edge
(i >= (256*1024)+10 && i < (256*1024)+10+10) || // double
(i >= (257*1024)+10 && i < (257*1024)+10+10) || // next to double
(i == (255*1024)+255) // quad
) {
string pixel_type = "x";
if (i >= 58000 && i < 58000+10) {
pixel_type = "b";
} else if (i >= 10 && i < 10+10) {
pixel_type = "e";
} else if (i >= 1024+10 && i < 1024+10+10) {
pixel_type = "ie";
} else if (i >= (256*1024)+10 && i < (256*1024)+10+10) {
pixel_type = "d";
} else if (i >= (257*1024)+10 && i < (257*1024)+10+10) {
pixel_type = "ntd";
} else if (i == (255*1024)+255) {
pixel_type = "q";
}
mapcanvas->SetLeftMargin(0.13);
mapcanvas->SetRightMargin(0.05);
grap_g1->SetMinimum(0);
grap_g1->SetMaximum(17000);
grap_g1->GetXaxis()->SetTitle("Tube current [mA]");
grap_g1->GetYaxis()->SetTitle("ADC [ADU]");
grap_g1->GetYaxis()->SetTitleOffset(0.9);
grap_g1->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
grap_g1->Draw("AP");
fit2->Draw("same");
mapcanvas->Update();
TPaveStats *st0 = (TPaveStats*)grap_g1->FindObject("stats");
st0->SetX1NDC(0.6);
st0->SetX2NDC(0.94);
st0->SetY1NDC(0.18);
st0->SetY2NDC(0.37);
st0->SetBorderSize(0);
st0->SetTextColor(kGreen+2);
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
vector<double> r1_adc_norm;
for (size_t j = 0; j < r1_adc.size(); j++) {
r1_adc_norm.push_back(r1_adc[j] - fit2->Eval(r1_filter[j]));
}
TGraphErrors *grap1_norm = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_filterer[0]),&(r1_adcer[0]));
grap1_norm->SetMarkerColor(kGreen+2);
grap1_norm->SetLineColor(kGreen+2);
TF1* flat_g1 = new TF1("flat_g1","0",rangemin1, rangemax1);
flat_g1->SetLineColor(kGreen+2);
TF1* lin_g1_p02pc = new TF1("lin_g1_p02pc","[0]+[1]*x",rangemin1,rangemax1);
lin_g1_p02pc->SetParameter(0,fit2->GetParameter(0)/500.);
lin_g1_p02pc->SetParameter(1,fit2->GetParameter(1)/500.);
@ -640,31 +760,10 @@ int main(int argc, char* argv[]) {
lin_g1_m05pc->SetParameter(1,fit2->GetParameter(1)/-200.);
lin_g1_m05pc->SetLineColor(kOrange+1);
grap0_norm->SetMarkerColor(kBlue);
grap0_norm->SetLineColor(kBlue);
flat_g0->SetLineColor(kBlue);
grap0_norm->GetXaxis()->SetTitle("Tube current [mA]");
grap0_norm->GetYaxis()->SetTitle("Normalised ADC [ADU]");
grap0_norm->GetYaxis()->SetTitleOffset(0.9);
grap0_norm->Draw("AP");
flat_g0->Draw("same");
lin_g0_p02pc->Draw("same");
lin_g0_p05pc->Draw("same");
lin_g0_m02pc->Draw("same");
lin_g0_m05pc->Draw("same");
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
grap0_norm->SetMinimum(-100);
grap0_norm->SetMaximum(100);
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g0norm_zoom_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
grap1_norm->SetMarkerColor(kGreen+2);
grap1_norm->SetLineColor(kGreen+2);
flat_g1->SetLineColor(kGreen+2);
grap1_norm->GetXaxis()->SetTitle("Tube current [mA]");
grap1_norm->GetYaxis()->SetTitle("Normalised ADC [ADU]");
grap1_norm->GetYaxis()->SetTitleOffset(0.9);
grap1_norm->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
grap1_norm->Draw("AP");
flat_g1->Draw("same");
lin_g1_p02pc->Draw("same");
@ -680,15 +779,48 @@ int main(int argc, char* argv[]) {
mapcanvas->SetLeftMargin(0.1);
mapcanvas->SetRightMargin(0.13);
delete grap1_norm;
delete flat_g1;
}
}
}
// get ratio measurements
if (rangemax0 > rangemin0 && rangemax1 > rangemin1) {
g0overg1map->Fill(i%NC,i/NC,fit1->GetParameter(1)/fit2->GetParameter(1));
g0overg1hist->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
g0overg1ermap->Fill(i%NC,i/NC,abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));
g0overg1erhist->Fill(abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2)));
if (isEdge(i)) {
g0overg1hist_isEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isInnerEdge(i)) {
g0overg1hist_isInnerEdge->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isDouble(i)) {
g0overg1hist_isDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isNextToDouble(i)) {
g0overg1hist_isNextToDouble->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isQuad(i)) {
g0overg1hist_isQuad->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
if (isBulk(i)) {
g0overg1hist_isBulk->Fill(fit1->GetParameter(1)/fit2->GetParameter(1));
}
}
delete fit1;
delete fit2;
delete grap_g0;
delete grap_g1;
}
}
}
}
range1minhist->GetXaxis()->SetTitle("Start of gain 1 fit [ADU]");
range1minhist->Draw();