detectors/JFCalibration
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
master
JFCalibration/DB_analysis.cpp
redford_s 5c3e77b091 Removed cut on adc error from array filling. 
With the new range finding, this is now taken into account.
Points with negative error will not be fit, since there will be an entry in the second array.
Not having these points in the first array would change the range for the second array.
2019-02-20 12:35:05 +01:00

893 lines
34 KiB
C++
Raw Permalink Blame History

#include "../sls_detector_calibration/jungfrauCommonHeader.h"
#include "../sls_detector_calibration/jungfrauCommonFunctions.h"
#include "../sls_detector_calibration/jungfrauFile.C"
#include "../sls_detector_calibration/jungfrauPixelMask.C"
#include "../sls_detector_calibration/jungfrauPedestal.C"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"
#include "TLegend.h"
int main(int argc, char* argv[]) {
jungfrauStyle();
gROOT->SetBatch(1);
gStyle->SetOptFit(11);
if (argc != 5) {
cout << "Correct usage:" << endl;
cout << "arg 1: specify module number" << endl;
cout << "arg 2: specify data location" << endl;
cout << "arg 3: specify pede file" << endl;
cout << "arg 4: specify data file" << endl;
cout << " " << endl;
exit(1);
}
string module_str = argv[1];
string data_loc = argv[2];
string pede_file = argv[3];
string data_file = argv[4];
jungfrauFile *thisfile = new jungfrauFile();
jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
bool pixel_mask [NCH];
jungfrauPedestal *pedestalObject = new jungfrauPedestal();
pedestalObject->pedestalSetNFrames(1000); // using 1000 frames, rolling window
TH2F* pedestalsG0 = new TH2F("pedestalsG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* pedestalsG1 = new TH2F("pedestalsG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* pedestalsG2 = new TH2F("pedestalsG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* pedeRMSG0 = new TH2F("pedeRMSG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* pedeRMSG1 = new TH2F("pedeRMSG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* pedeRMSG2 = new TH2F("pedeRMSG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* no_g0g1_calib_map = new TH2F("no_g0g1_calib_map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
char savename[128];
int createHistoFile = 0;
int createPerStepHistos = 0;
// open pede file
sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str());
thisfile->open((char*)savename, 0);
// calculate pixel mask
pixelMaskObject->initialisePixelMask(pixel_mask);
int framecounter = 0;
while (thisfile->readNextFrame()) {
framecounter++;
}
thisfile->rewind();
cout << "read " << framecounter << " frames" << endl;
for (int i = 0; i < 1000; i++) {
thisfile->readNextFrame();
pixelMaskObject->maskIfGainNot(0, thisfile->getFrameDataHandle(), pixel_mask);
}
for (int i = 0; i < 1000; i++) {
thisfile->readNextFrame();
pixelMaskObject->maskIfGainNot(1, thisfile->getFrameDataHandle(), pixel_mask);
}
for (int i = 0; i < 999; i++) {
thisfile->readNextFrame();
pixelMaskObject->maskIfGainNot(3, thisfile->getFrameDataHandle(), pixel_mask);
}
thisfile->rewind();
sprintf(savename,"plots/M%s/DirectBeam/pixelmask.png", module_str.c_str());
pixelMaskObject->plotPixelMask(pixel_mask,savename);
cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
// caluclate pedestals
for (int i = 0; i < 1000; i++) {
thisfile->readNextFrame();
pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
}
for (int i = 0; i < NCH; i++) {
if (pixel_mask[i] == true) {
pedestalsG0->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
pedeRMSG0->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
}
}
pedestalObject->pedestalClear();
for (int i = 0; i < 1000; i++) {
thisfile->readNextFrame();
pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
}
for (int i = 0; i < NCH; i++) {
if (pixel_mask[i] == true) {
pedestalsG1->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
pedeRMSG1->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
}
}
pedestalObject->pedestalClear();
for (int i = 0; i < 999; i++) {
thisfile->readNextFrame();
pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle());
}
for (int i = 0; i < NCH; i++) {
if (pixel_mask[i] == true) {
pedestalsG2->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i));
pedeRMSG2->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i));
}
}
pedestalObject->pedestalClear();
thisfile->close();
TCanvas *mapcanvas = new TCanvas("mapcanvas","",150,10,800,400);
mapcanvas->SetLeftMargin(0.1);
mapcanvas->SetRightMargin(0.13);
mapcanvas->SetTopMargin(0.08);
mapcanvas->SetBottomMargin(0.15);
TCanvas* c1 = new TCanvas("c1","");
mapcanvas->cd();
pedestalsG0->GetXaxis()->SetTitle("Column");
pedestalsG0->GetYaxis()->SetTitle("Row");
pedestalsG0->GetYaxis()->SetTitleOffset(0.7);
pedestalsG0->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeG0_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
pedestalsG1->GetXaxis()->SetTitle("Column");
pedestalsG1->GetYaxis()->SetTitle("Row");
pedestalsG1->GetYaxis()->SetTitleOffset(0.7);
pedestalsG1->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeG1_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
pedestalsG2->GetXaxis()->SetTitle("Column");
pedestalsG2->GetYaxis()->SetTitle("Row");
pedestalsG2->GetYaxis()->SetTitleOffset(0.7);
pedestalsG2->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeG2_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
pedeRMSG0->GetXaxis()->SetTitle("Column");
pedeRMSG0->GetYaxis()->SetTitle("Row");
pedeRMSG0->GetYaxis()->SetTitleOffset(0.7);
pedeRMSG0->GetZaxis()->SetRangeUser(0,25);
pedeRMSG0->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG0_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
pedeRMSG1->GetXaxis()->SetTitle("Column");
pedeRMSG1->GetYaxis()->SetTitle("Row");
pedeRMSG1->GetYaxis()->SetTitleOffset(0.7);
pedeRMSG1->GetZaxis()->SetRangeUser(0,15);
pedeRMSG1->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG1_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
pedeRMSG2->GetXaxis()->SetTitle("Column");
pedeRMSG2->GetYaxis()->SetTitle("Row");
pedeRMSG2->GetYaxis()->SetTitleOffset(0.7);
pedeRMSG2->GetZaxis()->SetRangeUser(0,15);
pedeRMSG2->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG2_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
jungfrauPedestal *pedestalObjectG0 = new jungfrauPedestal();
pedestalObjectG0->pedestalSetNFrames(1000); // using 1000 frames, rolling window
jungfrauPedestal *pedestalObjectG1 = new jungfrauPedestal();
pedestalObjectG1->pedestalSetNFrames(1000); // using 1000 frames, rolling window
int nsteps = 0;
if (module_str == "094" || module_str == "066") {
nsteps = 79;
} else {
nsteps = 40;
}
// 2D hist for every point, to save
// mean, error and number of fills
// can be loaded from a file
TH2F* avg_adc_g0_map [nsteps];
TH2F* avg_adc_g1_map [nsteps];
TH2F* avg_adcer_g0_map [nsteps];
TH2F* avg_adcer_g1_map [nsteps];
// creating the histograms
if (createHistoFile) {
for (int j = 0; j < nsteps; j++) {
sprintf(savename,"avg_adc_g0_map_%d", j);
avg_adc_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
sprintf(savename,"avg_adc_g1_map_%d", j);
avg_adc_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
sprintf(savename,"avg_adcer_g0_map_%d", j);
avg_adcer_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
sprintf(savename,"avg_adcer_g1_map_%d", j);
avg_adcer_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
}
int nfiles = ceil(nsteps/10.);
// total frame number check
int totalFrames = 0;
for (int filei = 0; filei < nfiles; filei++) {
sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str());
thisfile->open((char*)savename, filei);
while(thisfile->readNextFrame()) {
totalFrames++;
}
thisfile->close();
}
if (totalFrames == nsteps*1000) {
cout << "expected number of frames " << totalFrames << endl;
} else {
cout << "unexpected number of frames " << totalFrames << endl;
cout << "code expected nsteps*1000 = " << nsteps*1000 << " frames" << endl;
cout << "exiting to make sure you notice" << endl;
exit(1);
}
int framecounter = 0;
for (int filei = 0; filei < nfiles; filei++) {
// open data file
sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str());
thisfile->open((char*)savename, filei);
while (thisfile->readNextFrame()) {
if (module_str == "028") {
if (framecounter%1000 < 100) { // shutter skip
framecounter++;
continue;
}
}
pedestalObjectG0->addG0FrameToPedestalCalculation(thisfile->getFrameDataHandle());
pedestalObjectG1->addG1FrameToPedestalCalculation(thisfile->getFrameDataHandle());
framecounter++;
if (framecounter%1000 == 0) {
for (int i = 0; i < NCH; i++) {
avg_adc_g0_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG0->pedestalOfChannel(i));
avg_adcer_g0_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG0->semOfChannel(i));
avg_adc_g1_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG1->pedestalOfChannel(i));
avg_adcer_g1_map[framecounter/1000 -1]->SetBinContent((i%NC)+1,(i/NC)+1,pedestalObjectG1->semOfChannel(i));
}
pedestalObjectG0->pedestalClear();
pedestalObjectG1->pedestalClear();
if (createPerStepHistos == 1) {
avg_adc_g0_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
avg_adc_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
avg_adc_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
avg_adc_g0_map[framecounter/1000 -1]->Draw("colz");
avg_adc_g0_map[framecounter/1000 -1]->SetMinimum(0);
sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adc_g0_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
avg_adc_g1_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
avg_adc_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
avg_adc_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
avg_adc_g1_map[framecounter/1000 -1]->Draw("colz");
avg_adc_g1_map[framecounter/1000 -1]->SetMinimum(0);
sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adc_g1_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
avg_adcer_g0_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
avg_adcer_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
avg_adcer_g0_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
avg_adcer_g0_map[framecounter/1000 -1]->Draw("colz");
avg_adcer_g0_map[framecounter/1000 -1]->SetMinimum(0);
sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adcer_g0_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
avg_adcer_g1_map[framecounter/1000 -1]->GetXaxis()->SetTitle("Column");
avg_adcer_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitle("Row");
avg_adcer_g1_map[framecounter/1000 -1]->GetYaxis()->SetTitleOffset(0.7);
avg_adcer_g1_map[framecounter/1000 -1]->Draw("colz");
avg_adcer_g1_map[framecounter/1000 -1]->SetMinimum(0);
sprintf(savename,"plots/M%s/DirectBeam/perStep/avg_adcer_g1_map_%d_M%s.png", module_str.c_str(), ((framecounter/1000 -1)+1)*2, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
}
}
}
thisfile->close();
} // end of files
// save histograms
sprintf(savename,"/mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/M%s_CalibAna/DB_histos_M%s.root", module_str.c_str(), module_str.c_str());
TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
for (int i = 0; i < nsteps; i++) {
avg_adc_g0_map[i]->Write();
avg_adc_g1_map[i]->Write();
avg_adcer_g0_map[i]->Write();
avg_adcer_g1_map[i]->Write();
}
saved_file->Close();
} else {
// load histos
cout << "LOADING HISTOS" << endl;
sprintf(savename,"/mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/M%s_CalibAna/DB_histos_M%s.root", module_str.c_str(), module_str.c_str());
TFile* saved_file = new TFile((const char *)(savename),"READ");
for (int i = 0; i < nsteps; i++) {
sprintf(savename,"avg_adc_g0_map_%d", i);
avg_adc_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
sprintf(savename,"avg_adc_g1_map_%d", i);
avg_adc_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
sprintf(savename,"avg_adcer_g0_map_%d", i);
avg_adcer_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
sprintf(savename,"avg_adcer_g1_map_%d", i);
avg_adcer_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename));
}
}
TH1F* g0overg1hist = new TH1F("g0overg1hist","",100,-34,-24);
TH2F* g0map = new TH2F("g0map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g1map = new TH2F("g1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g0overg1map = new TH2F("g0overg1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g0ermap = new TH2F("g0ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g1ermap = new TH2F("g1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH1F* g0overg1erhist = new TH1F("g0overg1erhist","",100,0,2);
TH2F* g0overg1ermap = new TH2F("g0overg1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH1F* g0overg1hist_isEdge = new TH1F("g0overg1hist_isEdge","",100,-34,-24);
TH1F* g0overg1hist_isInnerEdge = new TH1F("g0overg1hist_isInnerEdge","",100,-34,-24);
TH1F* g0overg1hist_isDouble = new TH1F("g0overg1hist_isDouble","",100,-34,-24);
TH1F* g0overg1hist_isNextToDouble = new TH1F("g0overg1hist_isNextToDouble","",100,-34,-24);
TH1F* g0overg1hist_isQuad = new TH1F("g0overg1hist_isQuad","",100,-34,-24);
TH1F* g0overg1hist_isBulk = new TH1F("g0overg1hist_isBulk","",100,-34,-24);
TH1F* range0maxhist = new TH1F("range0maxhist","",100,-10.5,89.5);
TH1F* range1minhist = new TH1F("range1minhist","",100,-10.5,89.5);
TH2F* range0maxmap = new TH2F("range0maxmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* range1minmap = new TH2F("range1minmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
mapcanvas->cd();
mapcanvas->SetLeftMargin(0.13);
mapcanvas->SetRightMargin(0.05);
for (int i = 0; i < NCH; i++) {
vector<double> r0_adc;
vector<double> r0_adcer;
vector<double> r0_filter;
vector<double> r0_ferr;
vector<double> r1_adc;
vector<double> r1_adcer;
vector<double> r1_filter;
vector<double> r1_ferr;
for (int j = 0; j < nsteps; j++) {
int current = 0;
if (nsteps == 40) {
current = (j+1)*2;
} else if (nsteps == 79) {
current = j+2;
}
double this_g0avg = avg_adc_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
double this_g1avg = avg_adc_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
double this_g0avger = avg_adcer_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
double this_g1avger = avg_adcer_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1);
if (this_g0avg != 0) {
r0_filter.push_back(current);
r0_ferr.push_back(0);
r0_adc.push_back(this_g0avg);
r0_adcer.push_back(this_g0avger);
}
if (this_g1avg != 0) {
r1_filter.push_back(current);
r1_ferr.push_back(0);
r1_adc.push_back(this_g1avg);
r1_adcer.push_back(this_g1avger);
}
}
TGraphErrors *grap_g0 = 0;
TGraphErrors *grap_g1 = 0;
TF1 *fit_g0 = 0;
TF1 *fit_g1 = 0;
double rangemin_g0 = 0;
double rangemax_g0 = 0;
double rangemin_g1 = 0;
double rangemax_g1 = 0;
// define graphs
if (r0_adc.size() > 1) {
grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_ferr[0]),&(r0_adcer[0]));
grap_g0->SetMarkerStyle(20);
grap_g0->SetMarkerColor(kBlue);
grap_g0->SetLineColor(kBlue);
}
if (r1_adc.size() > 1) {
grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_ferr[0]),&(r1_adcer[0]));
grap_g1->SetMarkerStyle(20);
grap_g1->SetMarkerColor(kGreen+2);
grap_g1->SetLineColor(kGreen+2);
}
// plot the datapoints
if (r1_adc.size() > 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";
}
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->Draw("AP");
if (r0_adc.size() > 1) {
grap_g0->Draw("P");
}
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));
}
}
// define fit ranges and fit
if (r0_adc.size() > 1) {
rangemin_g0 = *min_element(r0_filter.begin(),r0_filter.end());
rangemax_g0 = highestPointBeforeSwitching(r0_filter,r1_filter);
range0maxhist->Fill(rangemax_g0);
range0maxmap->Fill(i%NC,i/NC,rangemax_g0);
if (rangemax_g0 > rangemin_g0) {
fit_g0 = new TF1("fit_g0","[0]+[1]*x",rangemin_g0,rangemax_g0);
fit_g0->SetParameter(0, 100.);
fit_g0->SetParameter(1, 1);
fit_g0->SetLineColor(kBlue);
fit_g0->SetParName(0,"G0 const");
fit_g0->SetParName(1,"G0 grad");
grap_g0->Fit(fit_g0,"QR","");
g0map->Fill(i%NC,i/NC,fit_g0->GetParameter(1));
g0ermap->Fill(i%NC,i/NC,fit_g0->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";
}
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");
fit_g0->Draw("same");
mapcanvas->Update();
TPaveStats *st0 = (TPaveStats*)grap_g0->FindObject("stats");
st0->SetX1NDC(0.2);
st0->SetX2NDC(0.54);
st0->SetY1NDC(0.18);
st0->SetY2NDC(0.37);
st0->SetBorderSize(0);
st0->SetTextColor(kBlue);
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));
vector<double> r0_adc_norm;
for (size_t j = 0; j < r0_adc.size(); j++) {
r0_adc_norm.push_back(r0_adc[j] - fit_g0->Eval(r0_filter[j]));
}
TGraphErrors *norm_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc_norm[0]),&(r0_ferr[0]),&(r0_adcer[0]));
norm_g0->SetMarkerColor(kBlue);
norm_g0->SetLineColor(kBlue);
TF1* flat_g0 = new TF1("flat_g0","0",rangemin_g0, rangemax_g0);
flat_g0->SetLineColor(kBlue);
TF1* lin_g0_p02pc = new TF1("lin_g0_p02pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
lin_g0_p02pc->SetParameter(0,fit_g0->GetParameter(0)/500.);
lin_g0_p02pc->SetParameter(1,fit_g0->GetParameter(1)/500.);
lin_g0_p02pc->SetLineColor(kOrange);
TF1* lin_g0_p05pc = new TF1("lin_g0_p05pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
lin_g0_p05pc->SetParameter(0,fit_g0->GetParameter(0)/200.);
lin_g0_p05pc->SetParameter(1,fit_g0->GetParameter(1)/200.);
lin_g0_p05pc->SetLineColor(kOrange+1);
TF1* lin_g0_m02pc = new TF1("lin_g0_m02pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
lin_g0_m02pc->SetParameter(0,fit_g0->GetParameter(0)/-500.);
lin_g0_m02pc->SetParameter(1,fit_g0->GetParameter(1)/-500.);
lin_g0_m02pc->SetLineColor(kOrange);
TF1* lin_g0_m05pc = new TF1("lin_g0_m05pc","[0]+[1]*x",rangemin_g0,rangemax_g0);
lin_g0_m05pc->SetParameter(0,fit_g0->GetParameter(0)/-200.);
lin_g0_m05pc->SetParameter(1,fit_g0->GetParameter(1)/-200.);
lin_g0_m05pc->SetLineColor(kOrange+1);
norm_g0->GetXaxis()->SetTitle("Tube current [mA]");
norm_g0->GetYaxis()->SetTitle("Normalised ADC [ADU]");
norm_g0->GetYaxis()->SetTitleOffset(0.9);
norm_g0->GetXaxis()->SetLimits(*min_element(r0_filter.begin(),r0_filter.end()),*max_element(r0_filter.begin(),r0_filter.end()));
norm_g0->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));
norm_g0->SetMinimum(-100);
norm_g0->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 norm_g0;
delete flat_g0;
}
}
}
if (r1_adc.size() > 1) {
rangemin_g1 = lowestPointAfterSwitching(r1_filter,r0_filter);
range1minhist->Fill(rangemin_g1);
range1minmap->Fill(i%NC,i/NC,rangemin_g1);
rangemax_g1 = *max_element(r1_filter.begin(),r1_filter.end());
if (rangemax_g1 > rangemin_g1) {
fit_g1 = new TF1("fit_g1","[0]+[1]*x",rangemin_g1,rangemax_g1);
fit_g1->SetParameter(0, 1000.);
fit_g1->SetParameter(1, -0.1);
fit_g1->SetLineColor(kGreen+2);
fit_g1->SetParName(0,"G1 const");
fit_g1->SetParName(1,"G1 grad");
grap_g1->Fit(fit_g1,"QR","");
g1map->Fill(i%NC,i/NC,fit_g1->GetParameter(1));
g1ermap->Fill(i%NC,i/NC,fit_g1->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";
}
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");
fit_g1->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] - fit_g1->Eval(r1_filter[j]));
}
TGraphErrors *norm_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc_norm[0]),&(r1_ferr[0]),&(r1_adcer[0]));
norm_g1->SetMarkerColor(kGreen+2);
norm_g1->SetLineColor(kGreen+2);
TF1* flat_g1 = new TF1("flat_g1","0",rangemin_g1, rangemax_g1);
flat_g1->SetLineColor(kGreen+2);
TF1* lin_g1_p02pc = new TF1("lin_g1_p02pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
lin_g1_p02pc->SetParameter(0,fit_g1->GetParameter(0)/500.);
lin_g1_p02pc->SetParameter(1,fit_g1->GetParameter(1)/500.);
lin_g1_p02pc->SetLineColor(kOrange);
TF1* lin_g1_p05pc = new TF1("lin_g1_p05pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
lin_g1_p05pc->SetParameter(0,fit_g1->GetParameter(0)/200.);
lin_g1_p05pc->SetParameter(1,fit_g1->GetParameter(1)/200.);
lin_g1_p05pc->SetLineColor(kOrange+1);
TF1* lin_g1_m02pc = new TF1("lin_g1_m02pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
lin_g1_m02pc->SetParameter(0,fit_g1->GetParameter(0)/-500.);
lin_g1_m02pc->SetParameter(1,fit_g1->GetParameter(1)/-500.);
lin_g1_m02pc->SetLineColor(kOrange);
TF1* lin_g1_m05pc = new TF1("lin_g1_m05pc","[0]+[1]*x",rangemin_g1,rangemax_g1);
lin_g1_m05pc->SetParameter(0,fit_g1->GetParameter(0)/-200.);
lin_g1_m05pc->SetParameter(1,fit_g1->GetParameter(1)/-200.);
lin_g1_m05pc->SetLineColor(kOrange+1);
norm_g1->GetXaxis()->SetTitle("Tube current [mA]");
norm_g1->GetYaxis()->SetTitle("Normalised ADC [ADU]");
norm_g1->GetYaxis()->SetTitleOffset(0.9);
norm_g1->GetXaxis()->SetLimits(*min_element(r1_filter.begin(),r1_filter.end()),*max_element(r1_filter.begin(),r1_filter.end()));
norm_g1->Draw("AP");
flat_g1->Draw("same");
lin_g1_p02pc->Draw("same");
lin_g1_p05pc->Draw("same");
lin_g1_m02pc->Draw("same");
lin_g1_m05pc->Draw("same");
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1norm_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
norm_g1->SetMinimum(-50);
norm_g1->SetMaximum(50);
sprintf(savename,"plots/M%s/DirectBeam/graps_%s_%d_g1norm_zoom_M%s.png", module_str.c_str(), pixel_type.c_str(), i, module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
delete norm_g1;
delete flat_g1;
}
}
}
// get ratio measurements
if (rangemax_g0 > rangemin_g0 && rangemax_g1 > rangemin_g1) {
double this_g0overg1 = fit_g0->GetParameter(1)/fit_g1->GetParameter(1);
double this_g0overg1er = abs(fit_g0->GetParameter(1)/fit_g1->GetParameter(1))*sqrt(pow((fit_g0->GetParError(1)/fit_g0->GetParameter(1)),2) + pow((fit_g1->GetParError(1)/fit_g1->GetParameter(1)),2));
g0overg1map->Fill(i%NC,i/NC,this_g0overg1);
g0overg1hist->Fill(this_g0overg1);
g0overg1ermap->Fill(i%NC,i/NC,this_g0overg1er);
g0overg1erhist->Fill(this_g0overg1er);
if (isEdge(i)) {
g0overg1hist_isEdge->Fill(this_g0overg1);
}
if (isInnerEdge(i)) {
g0overg1hist_isInnerEdge->Fill(this_g0overg1);
}
if (isDouble(i)) {
g0overg1hist_isDouble->Fill(this_g0overg1);
}
if (isNextToDouble(i)) {
g0overg1hist_isNextToDouble->Fill(this_g0overg1);
}
if (isQuad(i)) {
g0overg1hist_isQuad->Fill(this_g0overg1);
}
if (isBulk(i)) {
g0overg1hist_isBulk->Fill(this_g0overg1);
}
} else {
no_g0g1_calib_map->Fill(i%NC,i/NC,1);
}
delete fit_g0;
delete fit_g1;
delete grap_g0;
delete grap_g1;
}
mapcanvas->SetLeftMargin(0.1);
mapcanvas->SetRightMargin(0.13);
range1minhist->GetXaxis()->SetTitle("Start of gain 1 fit [ADU]");
range1minhist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/range1minhist_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
range0maxhist->GetXaxis()->SetTitle("End of gain 0 fit [ADU]");
range0maxhist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/range0maxhist_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
range1minmap->GetXaxis()->SetTitle("Column");
range1minmap->GetYaxis()->SetTitle("Row");
range1minmap->GetYaxis()->SetTitleOffset(0.7);
range1minmap->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/range1minmap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
range0maxmap->GetXaxis()->SetTitle("Column");
range0maxmap->GetYaxis()->SetTitle("Row");
range0maxmap->GetYaxis()->SetTitleOffset(0.7);
range0maxmap->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/range0maxmap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g1map->GetXaxis()->SetTitle("Column");
g1map->GetYaxis()->SetTitle("Row");
g1map->GetYaxis()->SetTitleOffset(0.7);
g1map->Draw("colz");
g1map->GetZaxis()->SetRangeUser(-40,0);
sprintf(savename,"plots/M%s/DirectBeam/g1map_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0map->GetXaxis()->SetTitle("Column");
g0map->GetYaxis()->SetTitle("Row");
g0map->GetYaxis()->SetTitleOffset(0.7);
g0map->Draw("colz");
g0map->GetZaxis()->SetRangeUser(200,800);
sprintf(savename,"plots/M%s/DirectBeam/g0map_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g1ermap->GetXaxis()->SetTitle("Column");
g1ermap->GetYaxis()->SetTitle("Row");
g1ermap->GetYaxis()->SetTitleOffset(0.7);
g1ermap->Draw("colz");
g1ermap->GetZaxis()->SetRangeUser(0,0.7);
sprintf(savename,"plots/M%s/DirectBeam/g1ermap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0ermap->GetXaxis()->SetTitle("Column");
g0ermap->GetYaxis()->SetTitle("Row");
g0ermap->GetYaxis()->SetTitleOffset(0.7);
g0ermap->Draw("colz");
g0ermap->GetZaxis()->SetRangeUser(2,8);
sprintf(savename,"plots/M%s/DirectBeam/g0ermap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0overg1map->GetXaxis()->SetTitle("Column");
g0overg1map->GetYaxis()->SetTitle("Row");
g0overg1map->GetYaxis()->SetTitleOffset(0.7);
g0overg1map->Draw("colz");
g0overg1map->GetZaxis()->SetRangeUser(-35,-20);
sprintf(savename,"plots/M%s/DirectBeam/g0overg1map_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0overg1ermap->GetXaxis()->SetTitle("Column");
g0overg1ermap->GetYaxis()->SetTitle("Row");
g0overg1ermap->GetYaxis()->SetTitleOffset(0.7);
g0overg1ermap->Draw("colz");
g0overg1ermap->GetZaxis()->SetRangeUser(0,1);
sprintf(savename,"plots/M%s/DirectBeam/g0overg1ermap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
no_g0g1_calib_map->GetXaxis()->SetTitle("Column");
no_g0g1_calib_map->GetYaxis()->SetTitle("Row");
no_g0g1_calib_map->GetYaxis()->SetTitleOffset(0.7);
no_g0g1_calib_map->Draw("colz");
sprintf(savename,"plots/M%s/DirectBeam/no_g0g1_calib_map_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
c1->cd();
g0overg1hist->GetXaxis()->SetTitle("G0 / G1");
g0overg1hist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
g0overg1erhist->GetXaxis()->SetTitle("#sigma_{G0 / G1}");
g0overg1erhist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/g0overg1erhist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
g0overg1hist_isEdge->SetLineColor(kBlue);
g0overg1hist_isInnerEdge->SetLineColor(kCyan);
g0overg1hist_isDouble->SetLineColor(kGreen+2);
g0overg1hist_isNextToDouble->SetLineColor(kRed);
g0overg1hist_isQuad->SetLineColor(kOrange);
g0overg1hist_isEdge->Scale(1./g0overg1hist_isEdge->GetEntries());
g0overg1hist_isInnerEdge->Scale(1./g0overg1hist_isInnerEdge->GetEntries());
g0overg1hist_isDouble->Scale(1./g0overg1hist_isDouble->GetEntries());
g0overg1hist_isNextToDouble->Scale(1./g0overg1hist_isNextToDouble->GetEntries());
g0overg1hist_isQuad->Scale(1./g0overg1hist_isQuad->GetEntries());
g0overg1hist_isBulk->Scale(1./g0overg1hist_isBulk->GetEntries());
TLegend *leg = new TLegend(0.62,0.6,0.93,0.93);
leg->AddEntry(g0overg1hist_isBulk, "Normal", "l");
leg->AddEntry(g0overg1hist_isDouble, "Double", "l");
leg->AddEntry(g0overg1hist_isNextToDouble, "Next to D", "l");
leg->AddEntry(g0overg1hist_isEdge, "Edge", "l");
leg->AddEntry(g0overg1hist_isInnerEdge, "Inner E", "l");
g0overg1hist_isDouble->GetXaxis()->SetTitle("G0 / G1");
g0overg1hist_isDouble->GetYaxis()->SetTitle("Normalised");
g0overg1hist_isDouble->GetYaxis()->SetTitleOffset(1.3);
g0overg1hist_isDouble->SetMinimum(0.0);
g0overg1hist_isDouble->SetMaximum(0.1);
g0overg1hist_isDouble->Draw();
g0overg1hist_isEdge->Draw("same");
g0overg1hist_isInnerEdge->Draw("same");
g0overg1hist_isNextToDouble->Draw("same");
g0overg1hist_isBulk->Draw("same");
leg->Draw("same");
sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist_perType_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
sprintf(savename,"data/M%s/DB_ratio_M%s.root", module_str.c_str(), module_str.c_str());
TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE");
g0overg1map->Write();
g0overg1ermap->Write();
saved_file2->Close();
}
Reference in New Issue View Git Blame Copy Permalink