364 lines
13 KiB
C++
364 lines
13 KiB
C++
// calculate pedestal correction of fluo data for a given storage cell, only G0
|
|
// (can be called multiple times to handle multi-storage-cell data, each time with different Sc number as argument, i.e. supports manual multi-threading on the same data set)
|
|
// make correction and save spectrum per pixel
|
|
// for fitting call CuFluo_analysis_sc_singlethread_fits (with correct SC number) afterward
|
|
|
|
#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 <vector>
|
|
|
|
int main(int argc, char* argv[]) {
|
|
|
|
jungfrauStyle();
|
|
gStyle->SetOptFit(11);
|
|
|
|
if (argc != 8) {
|
|
cout << "Correct usage:" << endl;
|
|
cout << "arg 1: specify module number" << endl;
|
|
cout << "arg 2: specify HG0 or G0" << endl;
|
|
cout << "arg 3: specify data location folder" << endl;
|
|
cout << "arg 4: specify pede file prefix" << endl;
|
|
cout << "arg 5: specify data file prefix" << endl;
|
|
cout << "arg 6: specify location of analysis root files" << endl;
|
|
cout << "arg 7: specify storage cell number" << endl;
|
|
cout << " " << endl;
|
|
exit(1);
|
|
}
|
|
|
|
string module_str = argv[1];
|
|
string gain_str = argv[2];
|
|
bool isJF11=false;
|
|
if (gain_str == "HG0JF11") {
|
|
gain_str = "HG0";
|
|
isJF11=true;
|
|
}
|
|
|
|
string data_loc = argv[3];
|
|
string pede_file = argv[4];
|
|
string data_file = argv[5];
|
|
string anadata_loc = argv[6];
|
|
|
|
bool save_to_afs = false;
|
|
|
|
int createHistoFile = 1;
|
|
|
|
char savename[256];
|
|
char plotfolder[256];
|
|
char rootdatafolder[256];
|
|
int filen = 352;
|
|
int pedefilen = 5;
|
|
int this_storagecell = stoi(argv[7]); // needs compiler flag -std=c++11
|
|
|
|
// create necessary directories with permissions drwxrwxr-x
|
|
// data/Mxxx
|
|
sprintf(savename,"data/M%s", module_str.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx
|
|
sprintf(savename,"plots/M%s", module_str.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo
|
|
sprintf(savename,"plots/M%s/CuFluo", module_str.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo/G0
|
|
sprintf(savename,"plots/M%s/CuFluo/G0", module_str.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo/HG0
|
|
sprintf(savename,"plots/M%s/CuFluo/HG0", module_str.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
sprintf(savename,"%s", anadata_loc.c_str());
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
|
|
if ( save_to_afs ) {
|
|
|
|
cout << "Saving to afs." << endl;
|
|
|
|
sprintf( plotfolder, "plots/M%s/CuFluo/%s", module_str.c_str(), gain_str.c_str() );
|
|
sprintf( rootdatafolder, "data/M%s", module_str.c_str() );
|
|
|
|
} else {
|
|
|
|
cout << "Not saving to afs." << endl;
|
|
|
|
// data/Mxxx
|
|
sprintf( savename,"%s/data", anadata_loc.c_str() );
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx
|
|
sprintf( savename,"%s/plots", anadata_loc.c_str() );
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo
|
|
sprintf( savename,"%s/plots/CuFluo", anadata_loc.c_str() );
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo/G0
|
|
sprintf( savename,"%s/plots/CuFluo/G0", anadata_loc.c_str() );
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
// plots/Mxxx/CuFluo/HG0
|
|
sprintf( savename,"%s/plots/CuFluo/HG0", anadata_loc.c_str() );
|
|
mkdir(savename, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
|
|
|
|
sprintf( plotfolder, "%s/plots/CuFluo/%s", anadata_loc.c_str(), gain_str.c_str() );
|
|
sprintf( rootdatafolder, "%s/data", anadata_loc.c_str() );
|
|
|
|
}
|
|
|
|
jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask();
|
|
bool pixel_mask [NCH];
|
|
pixelMaskObject->initialisePixelMask(pixel_mask);
|
|
|
|
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);
|
|
|
|
/*****************************************************************************************************************/
|
|
/* NOTE: Since arrays, as they were used before, tend to cause memory leaks if they grow too large on the stack, */
|
|
/* it might be smarter to realize all array objects with std::vector instead. */
|
|
/*****************************************************************************************************************/
|
|
|
|
if (createHistoFile) {
|
|
|
|
jungfrauFile *thisfile = new jungfrauFile();
|
|
|
|
jungfrauPedestal* pedestalObject_SC = new jungfrauPedestal();
|
|
pedestalObject_SC->pedestalSetNFrames(100);
|
|
|
|
vector<uint16_t> pedestals16_G0_start( NCH ); //Track pedestal shifting over the course of data taking
|
|
vector<double> pedeRMS16_G0( NCH );
|
|
|
|
//cout << "Last element of pedeRMS16_G0 = " << pedeRMS16_G0[NCH-1] << endl;
|
|
|
|
//double *pedeRMSary = pedeRMS16_G0.data();
|
|
|
|
//cout << "Last element of pedeRMSary = " << pedeRMSary[NCH-1] << endl;
|
|
|
|
{ //scope to limit the following 2 histograms
|
|
TH2F* pedestalsG0;
|
|
TH2F* pedeRMSG0;
|
|
|
|
Char_t *pedehistoname = new Char_t[50];
|
|
snprintf( pedehistoname, 50, "pedestalsG0_sc%d", this_storagecell );
|
|
pedestalsG0 = new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
|
|
snprintf( pedehistoname, 50, "pedeRMSG0_sc%d", this_storagecell );
|
|
pedeRMSG0 = new TH2F( pedehistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
|
|
|
|
for (int pedefilei = 0; pedefilei < pedefilen; ++pedefilei) {
|
|
|
|
// open pede file
|
|
sprintf( savename, "%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str() );
|
|
thisfile->open( (char*)savename, pedefilei );
|
|
|
|
//count events in file
|
|
int nevents = 0;
|
|
while ( thisfile->readNextFrame() ) {
|
|
nevents++;
|
|
}
|
|
thisfile->rewind();
|
|
cout << "read " << nevents << " events" << endl;
|
|
|
|
while ( thisfile->readNextFrame() ) {
|
|
// calculate pixel mask
|
|
pixelMaskObject->maskIfGainNot( 0, thisfile->getFrameDataHandle(), pixel_mask );
|
|
// caluclate pedestals
|
|
if( thisfile->currentSCnumber() == (uint64_t)this_storagecell ) {
|
|
pedestalObject_SC->addFrameToPedestalCalculation( thisfile->getFrameDataHandle() );
|
|
}
|
|
}
|
|
|
|
thisfile->close();
|
|
|
|
} // end of loops over files
|
|
|
|
for (int i = 0; i < NCH; ++i) {
|
|
if (pixel_mask[i] == true) {
|
|
pedestalsG0->Fill( i%NC, i/NC, pedestalObject_SC->pedestalOfChannel(i) );
|
|
pedeRMSG0->Fill( i%NC, i/NC, pedestalObject_SC->rmsOfChannel(i) );
|
|
}
|
|
}
|
|
pedestalObject_SC->pedestalData((uint16_t*)(pedestals16_G0_start.data())); //this loads the pedestal data into the array pedestals16_G0_start
|
|
pedestalObject_SC->pedestalRMSData(pedeRMS16_G0.data()); //same here
|
|
pedestalObject_SC->pedestalResetUpdates();
|
|
//pedestalObject[sci]->pedestalClear(); //I don't need to clear if I only calculate pedestal once.
|
|
|
|
sprintf( savename, "%s/pixelmask_%s_M%s.png", plotfolder, gain_str.c_str(), module_str.c_str() );
|
|
pixelMaskObject->plotPixelMask( pixel_mask, savename );
|
|
cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl;
|
|
|
|
mapcanvas->cd();
|
|
|
|
pedestalsG0->GetXaxis()->SetTitle("Column");
|
|
pedestalsG0->GetYaxis()->SetTitle("Row");
|
|
pedestalsG0->GetYaxis()->SetTitleOffset(0.7);
|
|
pedestalsG0->Draw("colz");
|
|
sprintf(savename,"%s/pedeG0_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str());
|
|
mapcanvas->SaveAs((const char *)(savename));
|
|
|
|
pedestalsG0->Delete(); //not needed hereafter (try to save memory)
|
|
|
|
pedeRMSG0->GetXaxis()->SetTitle("Column");
|
|
pedeRMSG0->GetYaxis()->SetTitle("Row");
|
|
pedeRMSG0->GetYaxis()->SetTitleOffset(0.7);
|
|
pedeRMSG0->GetZaxis()->SetRangeUser(0,30);
|
|
pedeRMSG0->Draw("colz");
|
|
sprintf(savename,"%s/pedeRMSG0_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str());
|
|
mapcanvas->SaveAs((const char *)(savename));
|
|
|
|
pedeRMSG0->Delete(); //not needed hereafter (try to save memory)
|
|
|
|
//Make sure histograms are deleted
|
|
|
|
} // end of scope to limit pedestalsG0 and pedeRMSG0
|
|
|
|
int adc2d_nbin= 1200;
|
|
if (isJF11) adc2d_nbin= 1600;
|
|
|
|
vector<TH2I*> adc2d; //declare
|
|
adc2d.reserve(8);
|
|
|
|
cout << "Initializing adc2d vector" << endl;
|
|
|
|
//initialize
|
|
for ( int adci = 0; adci < 8; ++adci ) {
|
|
Char_t *_histoname = new Char_t[50];
|
|
snprintf( _histoname, 50, "adc2d_%d_sc%d", adci + 1, this_storagecell );
|
|
adc2d.push_back( new TH2I( _histoname, "", adc2d_nbin, -200-0.5, adc2d_nbin-200-0.5, 65536, ( 65536*adci-0.5 ), ( 65536*( adci+1 )-0.5 ) ) );
|
|
}
|
|
|
|
cout << "Done." << endl;
|
|
|
|
//declare
|
|
TH1D* adcpc_spec;
|
|
TH2F* pede_updates;
|
|
TH2F* pede_diff;
|
|
|
|
//initialize
|
|
Char_t *_histoname = new Char_t[50];
|
|
snprintf( _histoname, 50, "adcpc_spec_sc%d", this_storagecell );
|
|
adcpc_spec = new TH1D( _histoname, "", 300, 0, 3000 ); //spectrum for every sc, every 10000 frames
|
|
snprintf( _histoname, 50, "pede_updates_sc%d", this_storagecell );
|
|
pede_updates = new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
|
|
snprintf( _histoname, 50, "pede_diff_sc%d", this_storagecell );
|
|
pede_diff = new TH2F( _histoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 );
|
|
|
|
for (int filei = 0; filei < filen; 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()) {
|
|
|
|
uint16_t* imagedptr = thisfile->getFrameDataHandle();
|
|
uint64_t scnumber = thisfile->currentSCnumber();
|
|
|
|
if ( scnumber == (uint64_t)this_storagecell ) {
|
|
pedestalObject_SC->addG0FrameToPedestalCalculationWThreshold( imagedptr, pedestalObject_SC, pedeRMS16_G0.data() );
|
|
|
|
for (int i = 0; i < NCH; i++) {
|
|
|
|
uint16_t gain = (imagedptr[i]&0xc000) >> 14;
|
|
|
|
if (gain == 0) {
|
|
|
|
int adcpc = (imagedptr[i]&0x3fff) - pedestalObject_SC->pedestalOfChannel(i);
|
|
adcpc_spec->Fill(adcpc);
|
|
|
|
if (i < (65536*1)) {
|
|
adc2d[0]->Fill(adcpc,i);
|
|
} else if (i < (65536*2)) {
|
|
adc2d[1]->Fill(adcpc,i);
|
|
} else if (i < (65536*3)) {
|
|
adc2d[2]->Fill(adcpc,i);
|
|
} else if (i < (65536*4)) {
|
|
adc2d[3]->Fill(adcpc,i);
|
|
} else if (i < (65536*5)) {
|
|
adc2d[4]->Fill(adcpc,i);
|
|
} else if (i < (65536*6)) {
|
|
adc2d[5]->Fill(adcpc,i);
|
|
} else if (i < (65536*7)) {
|
|
adc2d[6]->Fill(adcpc,i);
|
|
} else if (i < (65536*8)) {
|
|
adc2d[7]->Fill(adcpc,i);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
}
|
|
thisfile->close();
|
|
|
|
if ( filei%16 == 0 ) { // every 16th file equals every 10000 frames for each storage cell
|
|
|
|
adcpc_spec->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]");
|
|
adcpc_spec->Draw();
|
|
mapcanvas->SetLogy();
|
|
sprintf(savename,"%s/adcpc_spec_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
|
|
mapcanvas->SaveAs((const char *)(savename));
|
|
mapcanvas->SetLogy(0);
|
|
adcpc_spec->Reset();
|
|
|
|
pede_updates->Reset();
|
|
pede_diff->Reset();
|
|
for (int i = 0; i < NCH; i++) {
|
|
pede_updates->Fill( i%NC, i/NC, pedestalObject_SC->pedestalUpdates(i) );
|
|
pede_diff->Fill( i%NC, i/NC, pedestalObject_SC->pedestalOfChannel(i) - pedestals16_G0_start[i] );
|
|
pedestals16_G0_start[i] = pedestalObject_SC->pedestalOfChannel(i);
|
|
}
|
|
|
|
pede_updates->GetXaxis()->SetTitle("Column");
|
|
pede_updates->GetYaxis()->SetTitle("Row");
|
|
pede_updates->GetYaxis()->SetTitleOffset(0.7);
|
|
pede_updates->GetZaxis()->SetRangeUser(0,10000);
|
|
pede_updates->Draw("colz");
|
|
sprintf(savename,"%s/pede_updates_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
|
|
mapcanvas->SaveAs((const char *)(savename));
|
|
|
|
pede_diff->GetXaxis()->SetTitle("Column");
|
|
pede_diff->GetYaxis()->SetTitle("Row");
|
|
pede_diff->GetYaxis()->SetTitleOffset(0.7);
|
|
pede_diff->GetZaxis()->SetRangeUser(-40,40);
|
|
pede_diff->Draw("colz");
|
|
sprintf(savename,"%s/pede_diff_sc%d_slice%d_%s_M%s.png", plotfolder, this_storagecell, ( filei/16 ), gain_str.c_str(), module_str.c_str());
|
|
mapcanvas->SaveAs((const char *)(savename));
|
|
|
|
pedestalObject_SC->pedestalResetUpdates();
|
|
|
|
|
|
} //end of if filei%16
|
|
|
|
} // end of file loop
|
|
|
|
sprintf(savename,"%s/CuFluo_%s_sc%d_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), this_storagecell, filen-1);
|
|
TFile* saved_file = new TFile((const char *)(savename),"RECREATE");
|
|
for ( int adci = 0; adci < 8; ++adci ) {
|
|
adc2d[adci]->Write();
|
|
|
|
//In theory, I believe, adc2d can now also be destroyed since it has already been written to file
|
|
adc2d[adci]->Delete();
|
|
}
|
|
//All pede updates have been written, so they can be destroyed.
|
|
pede_updates->Delete();
|
|
pede_diff->Delete();
|
|
|
|
saved_file->Close();
|
|
|
|
} // end if creatHistoFile
|
|
|
|
|
|
}
|