detectors/JFCalibration
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
20a6b015baf6f7e79a4e282c38e692dfa0b9765c
JFCalibration/JFMC_CuFluoPeakFit.cpp
T

246 lines
9.0 KiB
C++
Raw Blame History

// file to fit fluo spectrum per pixel
// and save peak position and uncertainty
#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 "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/energyCalibration.h"
#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/energyCalibration.cpp"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"
int main(int argc, char* argv[]) {
jungfrauStyle();
gStyle->SetOptFit(11);
if (argc != 2) {
cout << "Correct usage:" << endl;
cout << "arg 1: specify module number" << endl;
cout << "eg: ./JFMC_FluoCalibFit 006" << endl;
cout << " " << endl;
exit(1);
}
string module_str = argv[1];
char histoname[128];
char savename[128];
TCanvas* c1 = new TCanvas("c1","");
sprintf(savename,"/mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/M%s_CalibAna/Fluo_comb.root", module_str.c_str());
TFile* comb_file = new TFile((const char *)(savename),"READ");
TH1F* peak_fit_pos = new TH1F("peak_fit_pos","",100,250,400);
TH1F* peak_fit_poserr = new TH1F("peak_fit_poserr","",100,0,2);
TH2F* peak_fit_pos_2d = new TH2F("peak_fit_pos_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* peak_fit_poserr_2d = new TH2F("peak_fit_poserr_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH1F* noise_fit_pos = new TH1F("noise_fit_pos","",100,-30,30);
TH1F* noise_fit_poserr = new TH1F("noise_fit_poserr","",100,0,0.1);
TH2F* noise_fit_pos_2d = new TH2F("noise_fit_pos_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* noise_fit_poserr_2d = new TH2F("noise_fit_poserr_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH1F* gain_fit = new TH1F("gain_fit","",100,250,400);
TH1F* gain_fiterr = new TH1F("gain_fiterr","",100,0,2);
TH2F* gain_fit_2d = new TH2F("gain_fit_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
TH2F* gain_fiterr_2d = new TH2F("gain_fiterr_2d","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5);
float noise_pos;
float noise_err;
for (int j=1; j<9; j++) {
cout << "slice " << j << endl;
sprintf(histoname,"adc2d_%i",j);
TH2I* adc2d_j = (TH2I*)(comb_file->Get((const char *)(histoname)));
cout << adc2d_j->GetEntries() << endl;
adc2d_j->Draw("colz");
c1->Update();
for (int i=(65536*(j-1)); i<(65536*(j)); i++) {
if (i%10000==0){cout << "another 10k" << endl;}
TH1D* proj = adc2d_j->ProjectionX("bin1",i-(65536*(j-1))+1,i-(65536*(j-1))+1);
proj->Draw();
if (proj->Integral(100,300)!=0 && proj->Integral(300,600)!=0) {
proj->Rebin(4);
proj->SetStats(kTRUE);
// noise
proj->GetXaxis()->SetRangeUser(-100,100);
proj->Fit("gaus","Q","",-50,50);
proj->Draw();
c1->Update();
TF1 *fit = proj->GetFunction("gaus");
noise_pos = fit->GetParameter(1);
noise_err = fit->GetParError(1);
noise_fit_pos->Fill(noise_pos);
noise_fit_pos_2d->Fill(i%NC,i/NC,noise_pos);
noise_fit_poserr->Fill(noise_err);
noise_fit_poserr_2d->Fill(i%NC,i/NC,noise_err);
if (i > 10000 && i < 10010) {
proj->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
TPaveStats *st0 = (TPaveStats*)proj->FindObject("stats");
st0->SetX1NDC(0.65);
st0->SetX2NDC(0.94);
st0->SetY1NDC(0.75);
st0->SetY2NDC(0.94);
st0->SetBorderSize(0);
sprintf(savename,"plots/M%s/CuFluo/proj_noise_%d.png",module_str.c_str(),i);
c1->SaveAs((const char *)(savename));
}
// peak
proj->GetXaxis()->SetRangeUser(100,400);
Double_t mypar[6];
mypar[0] = 0.0;
mypar[1] = 0.0;
mypar[2] = proj->GetBinCenter(proj->GetMaximumBin());
mypar[3] = 19.;
mypar[4] = 60.;
mypar[5] = 0.18;
Double_t emypar[6];
energyCalibration *thiscalibration = new energyCalibration();
thiscalibration->setScanSign(1);
thiscalibration->setStartParametersKb(mypar);
thiscalibration->fixParameter(0,0.); // no background
thiscalibration->fixParameter(1,0.);
TF1* fittedfun = thiscalibration->fitSpectrum(proj,mypar,emypar);
peak_fit_pos->Fill(mypar[2]);
peak_fit_poserr->Fill(emypar[2]);
peak_fit_pos_2d->Fill(i%NC,i/NC,mypar[2]);
peak_fit_poserr_2d->Fill(i%NC,i/NC,emypar[2]);
if (i > 10000 && i < 10010) {
proj->Draw();
fittedfun->Draw("same");
c1->Update();
proj->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
TPaveStats *st = (TPaveStats*)proj->FindObject("stats");
st->SetX1NDC(0.22);
st->SetX2NDC(0.55);
st->SetY1NDC(0.55);
st->SetY2NDC(0.94);
st->SetBorderSize(0);
sprintf(savename,"plots/M%s/CuFluo/proj_peak_%d.png",module_str.c_str(),i);
c1->SaveAs((const char *)(savename));
}
// gain
gain_fit->Fill(mypar[2] - noise_pos);
gain_fiterr->Fill(sqrt(pow(emypar[2],2) + pow(noise_err,2)));
gain_fit_2d->Fill(i%NC,i/NC,mypar[2] - noise_pos);
gain_fiterr_2d->Fill(i%NC,i/NC,sqrt(pow(emypar[2],2) + pow(noise_err,2)));
delete thiscalibration;
}
}
}
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);
peak_fit_pos_2d->GetXaxis()->SetTitle("Column");
peak_fit_pos_2d->GetYaxis()->SetTitle("Row");
peak_fit_pos_2d->GetYaxis()->SetTitleOffset(0.7);
peak_fit_pos_2d->Draw("colz");
peak_fit_pos_2d->GetZaxis()->SetRangeUser(250,400);
sprintf(savename,"plots/M%s/CuFluo/peak_fit_pos_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
peak_fit_poserr_2d->GetXaxis()->SetTitle("Column");
peak_fit_poserr_2d->GetYaxis()->SetTitle("Row");
peak_fit_poserr_2d->GetYaxis()->SetTitleOffset(0.7);
peak_fit_poserr_2d->Draw("colz");
peak_fit_poserr_2d->GetZaxis()->SetRangeUser(0,2);
sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
noise_fit_pos_2d->GetXaxis()->SetTitle("Column");
noise_fit_pos_2d->GetYaxis()->SetTitle("Row");
noise_fit_pos_2d->GetYaxis()->SetTitleOffset(0.7);
noise_fit_pos_2d->Draw("colz");
noise_fit_pos_2d->GetZaxis()->SetRangeUser(-30,30);
sprintf(savename,"plots/M%s/CuFluo/noise_fit_pos_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
noise_fit_poserr_2d->GetXaxis()->SetTitle("Column");
noise_fit_poserr_2d->GetYaxis()->SetTitle("Row");
noise_fit_poserr_2d->GetYaxis()->SetTitleOffset(0.7);
noise_fit_poserr_2d->Draw("colz");
noise_fit_poserr_2d->GetZaxis()->SetRangeUser(0,0.05);
sprintf(savename,"plots/M%s/CuFluo/noise_fit_poserr_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
gain_fit_2d->GetXaxis()->SetTitle("Column");
gain_fit_2d->GetYaxis()->SetTitle("Row");
gain_fit_2d->GetYaxis()->SetTitleOffset(0.7);
gain_fit_2d->Draw("colz");
gain_fit_2d->GetZaxis()->SetRangeUser(250,400);
sprintf(savename,"plots/M%s/CuFluo/gain_fit_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
gain_fiterr_2d->GetXaxis()->SetTitle("Column");
gain_fiterr_2d->GetYaxis()->SetTitle("Row");
gain_fiterr_2d->GetYaxis()->SetTitleOffset(0.7);
gain_fiterr_2d->Draw("colz");
gain_fiterr_2d->GetZaxis()->SetRangeUser(0,2);
sprintf(savename,"plots/M%s/CuFluo/gain_fiterr_2d.png",module_str.c_str());
mapcanvas->SaveAs((const char *)(savename));
c1->cd();
peak_fit_pos->GetXaxis()->SetTitle("Peak position [ADU]");
peak_fit_pos->Draw();
sprintf(savename,"plots/M%s/CuFluo/peak_fit_pos.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
peak_fit_poserr->GetXaxis()->SetTitle("Peak position uncert [ADU]");
peak_fit_poserr->Draw();
sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
noise_fit_pos->GetXaxis()->SetTitle("Noise position [ADU]");
noise_fit_pos->Draw();
sprintf(savename,"plots/M%s/CuFluo/noise_fit_pos.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
noise_fit_poserr->GetXaxis()->SetTitle("Noise position uncert [ADU]");
noise_fit_poserr->Draw();
sprintf(savename,"plots/M%s/CuFluo/noise_fit_poserr.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
gain_fit->GetXaxis()->SetTitle("Gain G0 [ADU / 8 keV]");
gain_fit->Draw();
sprintf(savename,"plots/M%s/CuFluo/gain_fit.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
gain_fiterr->GetXaxis()->SetTitle("Gain uncert [ADU / 8 keV]");
gain_fiterr->Draw();
sprintf(savename,"plots/M%s/CuFluo/gain_fiterr.png",module_str.c_str());
c1->SaveAs((const char *)(savename));
sprintf(savename,"data/M%s/CuF_gain_M%s.root", module_str.c_str(), module_str.c_str());
TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
gain_fit_2d->Write();
gain_fiterr_2d->Write();
saved_file->Close();
}
Reference in New Issue View Git Blame Copy Permalink