detectors/JFCalibration
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
main
JFCalibration/Fluo_fit.cpp
carulla_m_ 427b525022 BP_analysis.cpp add comments on the script 
BP_analysis_M431 modified
BP_analysis_M431.cpp modified
BP_scan.sh must declare module number as argument and readout speed (full_speed or half_speed) Modified the communication interface and commands to the pulser
BP_scan_both_speeds.sh Added Vb_comp 1220 at the end of the script and Modified the communication interface and commands to the pulser
CS_analysis.cpp add list of new modules for analysis
CS_analysis_JF11 modified
CS_analysis_M431.cpp modified
CS_fit.cpp add new modules in the list
CS_scan.sh must declare module number as argument and readout speed (full_speed or half_speed)
CS_scan_both_speeds.sh Added Vb_comp 1220 at the end of the script
CuFluo_exposure.sh Change HV to 60V and I to 40mA instead of 40V and 60mV+ must declare module number as argument and readout speed (full_speed or half_speed)
CuFluo_exposure_JF11.s Change HV to 60V and I to 40mA instead of 40V and 60mV
CuFluo_exposure_both_speeds.sh Change HV to 60V and I to 40mA instead of 40V and 60m
CuFluo_exposure_sc.sh Change HV to 60V and I to 40mA instead of 40V and 60m. Adjust the period to avoid losing packets
CuFluo_fit_sc.cpp Create a new file to fit Cu and noise peak fit for the storage cells data
Default_pixels_arrays.cpp Creates a default pixels system map given the module numbers of the system
Fluo_analysis_Ti_In.cpp Analysis the data for Fluorecensce target: Ti, Fe, Cu, Se, Mo, or In. The target needs to be passed as an argument.
Fluo_exposure_Ti_In.sh Fluorescence exposure with targets: Ti, Fe, Cu, Se, Mo, or In.
Fluo_fit.cpp Function to fit the fluorescence peaks for the different targets: Ti, Fe, Cu, Se, Mo, or In. The target needs to be passed as an argument.
JFMC_CalibWriter_wBP_sc.cpp change the range of the g0cutmap
Multi_filename_creator.sh Bash file to create 6 different receivers, one for each JFcalib readout board
filename_creator.sh add comments
makefile Add new objects:Fluo_analysis_Ti_In, Fluo_fit, CuFluo_fit_sc, CS_fit, BP_fit, etc...
2026-02-20 16:38:10 +01:00

722 lines
18 KiB
C++
Raw Permalink Blame History

// to analyse the backplane pulsing data per module
// changes by VH 210906: to eliminate hardcoded absolute paths, uses location of the analysis root files as additional input argument (accordingly changed in filename_creator.sh)
#include "TApplication.h"
#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 "sls_detector_calibration/energyCalibration.h"
#include "sls_detector_calibration/energyCalibration.cpp"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TF1.h"
#include "TFile.h"
#include "TPaveStats.h"
#include "TLegend.h"
#include "TPaveText.h"
#include <sys/stat.h>
//#include <sstream>
//#define NB_ENABLE 1
//void nonblock(int state);
//#define NB_DISABLE 0
//TApplication* rootapp;
//TCanvas *A2;
//TCanvas *A3;
//TCanvas *A4;
//TCanvas *A5;
//TCanvas *A6;
//TPaveStats *st0;
//void PlotCanvas(void);
double checkRangeMaxForAmplifierPlateau(double range_max) {
// check that the range maximum is no more than 6.4 V
// to avoid non-linearity coming from amplifier plateau
if (range_max > 6400) {
return 6400;
} else {
return range_max;
}
}
bool isHGX=false;
int main(int argc, char* argv[]) {
//nonblock(NB_ENABLE);
cout <<"opening the rootapp" <<endl;
TApplication rootapp("example",&argc, argv);
jungfrauStyle();
//gROOT->SetBatch(1);
gStyle->SetOptFit(11);
/*
if (argc != 3) {
cout << "Correct usage:" << endl;
cout << "arg 1: specify module number" << endl;
cout << "arg 2: specify data location" << endl;
cout << "arg 3: specify column (x)" << endl;
cout << "arg 4: specify row (y)" << endl;
exit(1);
}
*/ //uncomment for SR
if (argc != 6) {
cout << "Correct usage:" << endl;
cout << "arg 1: specify module number" << endl;
cout << "arg 2: specify fluorescence target" << endl;
cout << "arg 3: gain settings" << endl;
cout << "arg 4: specify pixel x position" << endl;
cout << "arg 5: specify pixel y position" << endl;
exit(1);
} //uncomment for VH 210906
string module_str = argv[1];
string gain_str = argv[3];
string target= argv[2];
string str2 =("HG0G1G2");
string str3 =("HGOG1G2"); //filename creator had this bug
int column = atoi(argv[4]);
int row = atoi(argv[5]);
int pixel = column+row*1024;
char data_loc[256];
sprintf(data_loc,"/mnt/sls_det_storage/jungfrau_calib/jungfrau_ana_sophie/M%s_%s_CalibAna/", module_str.c_str(),target.c_str());
cout << data_loc << endl;
std::string folder_path(data_loc);
if (folder_path.find(str2) != string::npos) isHGX=true;
if (folder_path.find(str3) != string::npos) isHGX=true;
bool isJF11=false;
if (gain_str == "HG0JF11") {
gain_str = "HG0";
isJF11=true;
}
// cout << data_loc.find(str2)<<" " << string::npos << " " << str2 << " " << data_loc <<endl;
if (isHGX) {
cout << " HG0->HG1->HG2 sequence - dynamicHG0" <<endl;
// plotfolder_str="BackplanePulsing_HG0G1G2";
}
else {
cout << " G0->G1->G2 sequence - dynamicG0" <<endl;
// plotfolder_str="BackplanePulsing";
}
//char savename[128]; //uncomment for SR
//char filename[128]; //uncomment for SR
char filename[256]; //uncomment for VH 210902
int low_ADU_peak = 0;
int high_ADU_peak = 0;
if (target=="Ti") {
if (gain_str == "HG0") {
low_ADU_peak = 350;
high_ADU_peak = 600;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="Fe") {
if (gain_str == "HG0") {
low_ADU_peak = 500;
high_ADU_peak = 800;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="Cu") {
if (gain_str == "HG0") {
low_ADU_peak = 700;
high_ADU_peak = 1000;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="Se") {
if (gain_str == "HG0") {
low_ADU_peak = 1050;
high_ADU_peak = 1350;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="Mo") {
if (gain_str == "HG0") {
low_ADU_peak = 1650;
high_ADU_peak = 2150;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="Ag") {
if (gain_str == "HG0") {
low_ADU_peak = 2100;
high_ADU_peak = 2650;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
if (target=="In") {
if (gain_str == "HG0") {
low_ADU_peak = 2300;
high_ADU_peak = 2800;
} else if (gain_str == "G0") {
low_ADU_peak = 250;
high_ADU_peak = 400;
}
}
int low_bin_noise = 100;
int high_bin_noise = 350;
int low_bin_peak = 0;
int high_bin_peak = 0;
int nrebin = 4;
if (target=="Ti") {
if (gain_str == "HG0") {
low_bin_peak = 450;
high_bin_peak = 750;
nrebin = 15;
} else if (gain_str == "G0") {
low_bin_peak = 250;
high_bin_peak = 400;
}
}
if (target=="Fe") {
if (gain_str == "HG0") {
low_bin_peak = 558;
high_bin_peak = 955;
nrebin = 8;
} else if (gain_str == "G0") {
low_bin_peak = 240;
high_bin_peak = 518;
}
}
if (target=="Cu") {
if (gain_str == "HG0") {
low_bin_peak = 701;
high_bin_peak = 1200;
nrebin = 8;
} else if (gain_str == "G0") {
low_bin_peak = 301;
high_bin_peak = 651;
}
}
if (target=="Se") {
if (gain_str == "HG0") {
low_bin_peak = 900;
high_bin_peak = 1500;
nrebin = 8;
} else if (gain_str == "G0") {
low_bin_peak = 420;
high_bin_peak = 908;
}
}
if (target=="Mo") {
if (gain_str == "HG0") {
low_bin_peak = 1400;
high_bin_peak = 2400;
nrebin = 8;
} else if (gain_str == "G0") {
low_bin_peak = 654;
high_bin_peak = 1414;
}
}
if (target=="Ag") {
if (gain_str == "HG0") {
low_bin_peak = 1900;
high_bin_peak = 2900;
nrebin = 15;
} else if (gain_str == "G0") {
low_bin_peak = 829;
high_bin_peak = 1793;
}
}
if (target=="In") {
if (gain_str == "HG0") {
low_bin_peak = 2000;
high_bin_peak = 3000;
nrebin = 15;
} else if (gain_str == "G0") {
low_bin_peak = 906;
high_bin_peak = 1959;
}
}
TCanvas *A2 = new TCanvas("A2","Plot Noise",150,10,800,400);
TCanvas *A3 = new TCanvas("A3","Plot Pedestal",150,10,800,400);
A2->SetLeftMargin(0.1);
A2->SetRightMargin(0.13);
A2->SetTopMargin(0.08);
A2->SetBottomMargin(0.15);
A3->SetLeftMargin(0.1);
A3->SetRightMargin(0.13);
A3->SetTopMargin(0.08);
A3->SetBottomMargin(0.15);
//file name
sprintf(filename,"%sCuFluo_%s_file0to21.root",folder_path.c_str(), gain_str.c_str());
cout << "Loading file " << filename << endl;
TFile* f = new TFile((const char *)(filename),"READ");
int chip;
int CH = row*1024+column;
if (CH < (65536*1)) {
chip = 1;
} else if (CH < (65536*2)) {
chip = 2;
} else if (CH < (65536*3)) {
chip = 3;
} else if (CH < (65536*4)) {
chip = 4;
} else if (CH < (65536*5)) {
chip = 5;
} else if (CH < (65536*6)) {
chip = 6;
} else if (CH < (65536*7)) {
chip = 7;
} else if (CH < (65536*8)){
chip = 8;
}
cout<< "Chip "<< chip << " Channel number " << CH << endl;
TH2I* hist0=(TH2I*)f->Get(Form("adc2d_%d",chip));
cout << "Creating histogram for " << Form("adc2d_%d",chip) << endl;
TH1D* proj = hist0->ProjectionX("bin1",CH-(65536*(chip-1))+1,CH-(65536*(chip-1))+1);
cout << "Data for pixel "<< pixel << " is loaded" << endl;
if (proj->Integral(low_bin_noise,high_bin_noise)!=0 && proj->Integral(low_bin_peak,high_bin_peak)!=0) {
A2->cd();
cout << "Canvas Noise" << endl;
// noise
TH1D *proj_noise = dynamic_cast<TH1D*>(proj->Rebin(nrebin,"proj_noise"));
proj_noise->SetStats(kTRUE);
proj_noise->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_noise),proj->GetBinLowEdge(high_bin_noise+1));
proj_noise->Fit("gaus","Q");
TF1 *fit = proj_noise->GetFunction("gaus");
proj_noise->Draw();
A2->Update();
proj_noise->GetXaxis()->SetTitle("Noise RMS [ADU]");
proj_noise->GetXaxis()->SetRangeUser(-200,200);
fit->SetParNames("N_{#gamma}", "Peak pos", "Noise RMS");
TPaveStats *st0 = (TPaveStats*)proj_noise->FindObject("stats");
st0->SetX1NDC(0.53);
st0->SetX2NDC(0.94);
st0->SetY1NDC(0.75);
st0->SetY2NDC(0.94);
st0->SetBorderSize(0);
st0->SetTextSize(0.04);
A2->Modified();
A2->Update();
A3->cd();
// peak
TH1D *proj_peak = dynamic_cast<TH1D*>(proj->Rebin(nrebin,"proj_peak"));
proj_peak->SetStats(kTRUE);
proj_peak->GetXaxis()->SetRangeUser(proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
Double_t mypar[8];
mypar[0] = 0.0;
mypar[1] = 0.0;
if (target=="Ti") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 24.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.12;
}
// Kb/Ka position ratio
mypar[6] = 1.0;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.2;
}
}
if (target=="Fe") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 27.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.2;
}
// Kb/Ka position ratio
mypar[6] = 1.11;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.12;
}
}
if (target=="Cu") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 29.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.14;
}
// Kb/Ka position ratio
mypar[6] = 1.11; // mypar[6] = 1.12;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.14;
}
}
if (target=="Se") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 29.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.14;
}
// Kb/Ka position ratio
mypar[6] = 1.11;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.14;
}
}
if (target=="Mo") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 29.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.14;
}
// Kb/Ka position ratio
mypar[6] = 1.12;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.14;
}
}
if (target=="Ag") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 29.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.14;
}
// Kb/Ka position ratio
mypar[6] = 1.10;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.14;
}
}
if (target=="In") {
// Peak position
mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin());
// Noise RMS
if (gain_str == "G0") {
mypar[3] = 16.;
} else if (gain_str == "HG0") {
mypar[3] = 29.;
}
// Peak hight
mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin());
// Charge sharing
if (gain_str == "G0") {
mypar[5] = 0.17;
} else if (gain_str == "HG0") {
mypar[5] = 0.14;
}
// Kb/Ka position ratio
mypar[6] = 1.13;
// Kb fraction
if (gain_str == "G0") {
mypar[7] = 0.12;
} else if (gain_str == "HG0") {
mypar[7] = 0.14;
}
}
Double_t emypar[8];
energyCalibration *thiscalibration = new energyCalibration();
thiscalibration->setScanSign(1);
thiscalibration->setStartParametersKb(mypar);
thiscalibration->fixParameter(0,0.); // no background
thiscalibration->fixParameter(1,0.);
TF1* fittedfun = thiscalibration->fitSpectrumKb(proj_peak,mypar,emypar);
TF1 *gaus_Ka = new TF1("gaus_Ka","gaus",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
gaus_Ka->SetParameters(mypar[4],mypar[2],mypar[3]);
gaus_Ka->SetLineColor(kBlue);
TF1 *erfc_Ka = new TF1("erfc_Ka","[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
erfc_Ka->SetParameters(mypar[4]*mypar[5], 1, mypar[2], mypar[3]);
erfc_Ka->SetLineColor(kOrange);
TF1 *gaus_Kb = new TF1("gaus_Kb","gaus",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
gaus_Kb->SetParameters(mypar[4]*mypar[7],mypar[6]*mypar[2],mypar[3]);
gaus_Kb->SetLineColor(kGreen+2);
TF1 *erfc_Kb = new TF1("erfc_Kb","[0]/2.*(TMath::Erfc(([1]*(x-[2])/[3])/(TMath::Sqrt(2.))))",proj->GetBinLowEdge(low_bin_peak),proj->GetBinLowEdge(high_bin_peak+1));
erfc_Kb->SetParameters(mypar[4]*mypar[7]*mypar[5], 1, mypar[6]*mypar[2], mypar[3]);
erfc_Kb->SetLineColor(kOrange+7);
proj_peak->Draw();
A3->Update();
erfc_Kb->Draw("same");
erfc_Ka->Draw("same");
gaus_Kb->Draw("same");
gaus_Ka->Draw("same");
fittedfun->Draw("same");
A3->Update();
proj_peak->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
fittedfun->SetParNames("Bkg height", "Bkg grad", "K_{#alpha} pos", "Noise RMS", "K_{#alpha} height", "CS", "K_{#beta}/K_{#alpha} pos", "K_{#beta} frac");
TPaveStats *st = (TPaveStats*)proj_peak->FindObject("stats");
st->SetX1NDC(0.15);
st->SetX2NDC(0.55);
st->SetY1NDC(0.7);
st->SetY2NDC(0.94);
st->SetBorderSize(0);
st->SetTextSize(0.04);
A3->Modified();
A3->Update();
} else {
std::cout << "Masked pixel";
}
rootapp.Run();
return 0;
}
Reference in New Issue View Git Blame Copy Permalink