detectors/JFCalibration
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
63dd95f1c51a8e3bbdcc09962fbf6ba4bb87f6b9
JFCalibration/JFMC_CalibWriter.cpp

320 lines
14 KiB
C++
Raw Blame History

// file to combine the three sets of results
// to write the calibration constants for the junfgrau module
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonHeader.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauCommonFunctions.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauFile.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPixelMask.C"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPedestal.C"
#include "TFile.h"
struct GainMaps
{
static double g0VALs[NCH]; // declaration, incomplete type
static double g1VALs[NCH];
static double g2VALs[NCH];
};
double GainMaps::g0VALs[NCH]; // definition, complete type
double GainMaps::g1VALs[NCH];
double GainMaps::g2VALs[NCH];
int main(int argc, char* argv[]) {
jungfrauStyle();
if (argc != 2) {
cout << "Correct usage:" << endl;
cout << "arg 1: specify module number" << endl;
cout << "eg: ./JFMC_CalibWriter 006" << endl;
cout << " " << endl;
exit(1);
}
string module_str = argv[1];
char savename[128];
// Fluo dataset
sprintf(savename,"data/M%s/CuF_peakpos_M%s.root", module_str.c_str(), module_str.c_str());
TFile* FL_file = new TFile((char*)savename,"READ");
TH2F* FL_peak_map = (TH2F*)FL_file->Get("peak_fit_pos_2d");
TH2F* FL_peaker_map = (TH2F*)FL_file->Get("peak_fit_poserr_2d");
// Direct beam dataset
sprintf(savename,"data/M%s/DB_ratio_M%s.root", module_str.c_str(), module_str.c_str());
TFile* DB_file = new TFile((char*)savename,"READ");
TH2F* DB_ratio_map = (TH2F*)DB_file->Get("g0overg1map");
TH2F* DB_ratioer_map = (TH2F*)DB_file->Get("g0overg1ermap");
// Current source scan dataset
sprintf(savename,"data/M%s/CS_ratio_M%s.root", module_str.c_str(), module_str.c_str());
TFile* CS_file = new TFile((char*)savename,"READ");
TH2F* CS_ratio_map = (TH2F*)CS_file->Get("g1overg2_map");
TH2F* CS_ratioer_map = (TH2F*)CS_file->Get("g1overg2er_map");
// Now write binary maps
// units of ADU/keV
// one output file per module with 3 maps in
GainMaps *mapsObject = new GainMaps;
TH1F* g0histall = new TH1F("g0histall","",100,20,60);
TH1F* g1histall = new TH1F("g1histall","",100,-4,0);
TH1F* g2histall = new TH1F("g2histall","",100,-0.5,0);
TH1F* g0ferhistall = new TH1F("g0ferhistall","",100,0,0.005);
TH1F* g1ferhistall = new TH1F("g1ferhistall","",100,0,0.05);
TH1F* g2ferhistall = new TH1F("g2ferhistall","",100,0,0.05);
TH1F* g0histcut = new TH1F("g0histcut","",100,20,60);
TH1F* g1histcut = new TH1F("g1histcut","",100,-4,0);
TH1F* g2histcut = new TH1F("g2histcut","",100,-0.5,0);
TH2F* g0mapall = new TH2F("g0mapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g1mapall = new TH2F("g1mapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g2mapall = new TH2F("g2mapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g0fermapall = new TH2F("g0fermapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g1fermapall = new TH2F("g1fermapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g2fermapall = new TH2F("g2fermapall","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g0mapcut = new TH2F("g0mapcut","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g1mapcut = new TH2F("g1mapcut","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* g2mapcut = new TH2F("g2mapcut","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* defaultvalmap = new TH2F("defaultvalmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* maskmap = new TH2F("maskmap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
sprintf(savename,"data/M%s/gainMaps_M%s.bin", module_str.c_str(), module_str.c_str());
fstream outfile;
outfile.open(savename, ios::binary | ios::out);
for (int i = 0; i < NCH; i++) {
// first make initial maps and histograms
// check for masked pixels
if (FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0 &&
DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0 &&
CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0) {
mapsObject->g0VALs[i] = FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1) / 8.0;
mapsObject->g1VALs[i] = mapsObject->g0VALs[i] / DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1);
mapsObject->g2VALs[i] = mapsObject->g1VALs[i] / CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1);
} else {
// 0 in masked pixels: calculated energy will be 0?
mapsObject->g0VALs[i] = 0.;
mapsObject->g1VALs[i] = 0.;
mapsObject->g2VALs[i] = 0.;
maskmap->Fill(i%NC,i/NC,1);
}
g0histall->Fill(mapsObject->g0VALs[i]);
g1histall->Fill(mapsObject->g1VALs[i]);
g2histall->Fill(mapsObject->g2VALs[i]);
g0mapall->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g0VALs[i]);
g1mapall->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g1VALs[i]);
g2mapall->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g2VALs[i]);
if (FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0) {
g0fermapall->SetBinContent((i%NC)+1,(i/NC)+1,FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1));
g0ferhistall->Fill(FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1));
}
if (FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0 && DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0) {
g1fermapall->SetBinContent((i%NC)+1,(i/NC)+1, sqrt(pow((FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((DB_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2)));
g1ferhistall->Fill(sqrt(pow((FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((DB_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2)));
}
if (FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0 &&
DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0 &&
CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1) != 0) {
g2fermapall->SetBinContent((i%NC)+1,(i/NC)+1, sqrt(pow((FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((DB_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((CS_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2)));
g2ferhistall->Fill(sqrt(pow((FL_peaker_map->GetBinContent((i%NC)+1,(i/NC)+1)/FL_peak_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((DB_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/DB_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2) + pow((CS_ratioer_map->GetBinContent((i%NC)+1,(i/NC)+1)/CS_ratio_map->GetBinContent((i%NC)+1,(i/NC)+1)),2)));
}
}
float defaultG0 = g0histall->GetBinCenter(g0histall->GetMaximumBin());
float defaultG1 = g1histall->GetBinCenter(g1histall->GetMaximumBin());
float defaultG2 = g2histall->GetBinCenter(g2histall->GetMaximumBin());
cout << "default values for module " << module_str << " " << defaultG0 << " " << defaultG1 << " " << defaultG2 << endl;
for (int i = 0; i < NCH; i++) {
// crazy value check: replace with default
// think of a way to get the limits better
if (mapsObject->g0VALs[i] < 30 || mapsObject->g0VALs[i] > 48 ||
mapsObject->g1VALs[i] < -2.0 || mapsObject->g1VALs[i] > -0.5 ||
mapsObject->g2VALs[i] < -0.2 || mapsObject->g2VALs[i] > -0.02) {
mapsObject->g0VALs[i] = defaultG0;
mapsObject->g1VALs[i] = defaultG1;
mapsObject->g2VALs[i] = defaultG2;
defaultvalmap->Fill(i%NC,i/NC,1);
}
g0histcut->Fill(mapsObject->g0VALs[i]);
g1histcut->Fill(mapsObject->g1VALs[i]);
g2histcut->Fill(mapsObject->g2VALs[i]);
g0mapcut->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g0VALs[i]);
g1mapcut->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g1VALs[i]);
g2mapcut->SetBinContent((i%NC)+1,(i/NC)+1,mapsObject->g2VALs[i]);
}
outfile.write((char*)mapsObject->g0VALs, sizeof(mapsObject->g0VALs));
outfile.write((char*)mapsObject->g1VALs, sizeof(mapsObject->g1VALs));
outfile.write((char*)mapsObject->g2VALs, sizeof(mapsObject->g2VALs));
outfile.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);
g0mapall->GetXaxis()->SetTitle("Column");
g0mapall->GetYaxis()->SetTitle("Row");
g0mapall->GetYaxis()->SetTitleOffset(0.7);
g0mapall->Draw("colz");
g0mapall->GetZaxis()->SetRangeUser(30,50);
sprintf(savename,"plots/M%s/gain0_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0fermapall->GetXaxis()->SetTitle("Column");
g0fermapall->GetYaxis()->SetTitle("Row");
g0fermapall->GetYaxis()->SetTitleOffset(0.7);
g0fermapall->Draw("colz");
g0fermapall->GetZaxis()->SetRangeUser(0.,0.005);
sprintf(savename,"plots/M%s/gain0fer_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g1mapall->GetXaxis()->SetTitle("Column");
g1mapall->GetYaxis()->SetTitle("Row");
g1mapall->GetYaxis()->SetTitleOffset(0.7);
g1mapall->Draw("colz");
g1mapall->GetZaxis()->SetRangeUser(-2,-1);
sprintf(savename,"plots/M%s/gain1_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g1fermapall->GetXaxis()->SetTitle("Column");
g1fermapall->GetYaxis()->SetTitle("Row");
g1fermapall->GetYaxis()->SetTitleOffset(0.7);
g1fermapall->Draw("colz");
g1fermapall->GetZaxis()->SetRangeUser(0,0.1);
sprintf(savename,"plots/M%s/gain1fer_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g2mapall->GetXaxis()->SetTitle("Column");
g2mapall->GetYaxis()->SetTitle("Row");
g2mapall->GetYaxis()->SetTitleOffset(0.7);
g2mapall->Draw("colz");
g2mapall->GetZaxis()->SetRangeUser(-0.15,-0.05);
sprintf(savename,"plots/M%s/gain2_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g2fermapall->GetXaxis()->SetTitle("Column");
g2fermapall->GetYaxis()->SetTitle("Row");
g2fermapall->GetYaxis()->SetTitleOffset(0.7);
g2fermapall->Draw("colz");
g2fermapall->GetZaxis()->SetRangeUser(0,0.1);
sprintf(savename,"plots/M%s/gain2fer_all_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g0mapcut->GetXaxis()->SetTitle("Column");
g0mapcut->GetYaxis()->SetTitle("Row");
g0mapcut->GetYaxis()->SetTitleOffset(0.7);
g0mapcut->Draw("colz");
g0mapcut->GetZaxis()->SetRangeUser(30,50);
sprintf(savename,"plots/M%s/gain0_cut_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g1mapcut->GetXaxis()->SetTitle("Column");
g1mapcut->GetYaxis()->SetTitle("Row");
g1mapcut->GetYaxis()->SetTitleOffset(0.7);
g1mapcut->Draw("colz");
g1mapcut->GetZaxis()->SetRangeUser(-2,-1);
sprintf(savename,"plots/M%s/gain1_cut_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
g2mapcut->GetXaxis()->SetTitle("Column");
g2mapcut->GetYaxis()->SetTitle("Row");
g2mapcut->GetYaxis()->SetTitleOffset(0.7);
g2mapcut->Draw("colz");
g2mapcut->GetZaxis()->SetRangeUser(-0.15,-0.05);
sprintf(savename,"plots/M%s/gain2_cut_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
defaultvalmap->GetXaxis()->SetTitle("Column");
defaultvalmap->GetYaxis()->SetTitle("Row");
defaultvalmap->GetYaxis()->SetTitleOffset(0.7);
defaultvalmap->Draw("colz");
sprintf(savename,"plots/M%s/defaultvalmap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
maskmap->GetXaxis()->SetTitle("Column");
maskmap->GetYaxis()->SetTitle("Row");
maskmap->GetYaxis()->SetTitleOffset(0.7);
maskmap->Draw("colz");
sprintf(savename,"plots/M%s/maskmap_M%s.png", module_str.c_str(), module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
TCanvas* c1 = new TCanvas("c1","");
g0histcut->SetLineColor(kRed);
g1histcut->SetLineColor(kRed);
g2histcut->SetLineColor(kRed);
g0histcut->GetXaxis()->SetTitle("Gain G0 [ADU/keV]");
g0histcut->Draw();
g0histall->Draw("same");
sprintf(savename,"plots/M%s/g0hist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy();
sprintf(savename,"plots/M%s/g0hist_log_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy(0);
g1histcut->GetXaxis()->SetTitle("Gain G1 [ADU/keV]");
g1histcut->Draw();
g1histall->Draw("same");
sprintf(savename,"plots/M%s/g1hist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy();
sprintf(savename,"plots/M%s/g1hist_log_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy(0);
g2histcut->GetXaxis()->SetTitle("Gain G2 [ADU/keV]");
g2histcut->Draw();
g2histall->Draw("same");
sprintf(savename,"plots/M%s/g2hist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy();
sprintf(savename,"plots/M%s/g2hist_log_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
c1->SetLogy(0);
g0ferhistall->GetXaxis()->SetTitle("Gain G0 frac uncert");
g0ferhistall->Draw();
sprintf(savename,"plots/M%s/g0ferhist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
g1ferhistall->GetXaxis()->SetTitle("Gain G1 frac uncert");
g1ferhistall->Draw();
sprintf(savename,"plots/M%s/g1ferhist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
g2ferhistall->GetXaxis()->SetTitle("Gain G2 frac uncert");
g2ferhistall->Draw();
sprintf(savename,"plots/M%s/g2ferhist_M%s.png", module_str.c_str(), module_str.c_str());
c1->SaveAs((const char *)(savename));
}
Reference in New Issue View Git Blame Copy Permalink