From 8ca78cb101996fd0b15d959f47f15035b7859cfe Mon Sep 17 00:00:00 2001 From: hinger_v Date: Thu, 19 May 2022 10:29:43 +0200 Subject: [PATCH] Split SC analysis in two files handling raw analysis and fitting, each script handles a single storage cell at a time allowing multi-threading --- CuFluo_analysis_sc_singlethread_data.cpp | 377 ++++++++++ CuFluo_analysis_sc_singlethread_fits.cpp | 846 +++++++++++++++++++++++ makefile | 8 +- 3 files changed, 1230 insertions(+), 1 deletion(-) create mode 100644 CuFluo_analysis_sc_singlethread_data.cpp create mode 100644 CuFluo_analysis_sc_singlethread_fits.cpp diff --git a/CuFluo_analysis_sc_singlethread_data.cpp b/CuFluo_analysis_sc_singlethread_data.cpp new file mode 100644 index 0000000..39985ee --- /dev/null +++ b/CuFluo_analysis_sc_singlethread_data.cpp @@ -0,0 +1,377 @@ +// file to calculate pedestal correction of fluo data for all storage cells, only G0 +// make correction and save spectrum per pixel +// then fit fluo spectrum per pixel +// and save peak position and uncertainty +// 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 "../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 "TObjArray.h" + +#include +#include + +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]; //uncomment for VH 210906 + + bool save_to_afs = false; + + int createHistoFile = 1; + + //char histoname[128]; + //char savename[128]; + //char histoname[256]; // VH 210902 + char savename[256]; // VH 210902 + char plotfolder[256]; + char rootdatafolder[256]; + int filen = 352; + int pedefilen = 5; + int this_storagecell = stoi(argv[7]); + + // 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); + // /mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/Mxxx_CalibAna + //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna", module_str.c_str()); //uncomment for SR + sprintf(savename,"%s", anadata_loc.c_str()); //uncomment for VH 210906 + 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); + + //TCanvas* c1 = new TCanvas("c1",""); + + /*****************************************************************************************************************/ + /* 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 pedestals16_G0_start( NCH ); //I suppose, this is to track pedestal shifting over the course of data taking + vector 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() ); //VH: note, this adds a double slash in the filepath + 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 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; + + //cout << "Initialize vectors of histograms..." << endl; + + //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 ); + + //cout << "Done." << endl; + + 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()); //VH: note, this adds a double slash in the filepath + 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,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna/CuFluo_%s_file0to%d.root", module_str.c_str(), gain_str.c_str(), filen-1); + sprintf(savename,"%s/CuFluo_%s_sc%d_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), this_storagecell, filen-1); //uncomment for VH 210906 + 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 i creatHistoFile + + +} diff --git a/CuFluo_analysis_sc_singlethread_fits.cpp b/CuFluo_analysis_sc_singlethread_fits.cpp new file mode 100644 index 0000000..93165f8 --- /dev/null +++ b/CuFluo_analysis_sc_singlethread_fits.cpp @@ -0,0 +1,846 @@ +// file to calculate pedestal correction of fluo data for all storage cells, only G0 +// make correction and save spectrum per pixel +// then fit fluo spectrum per pixel +// and save peak position and uncertainty +// 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 "../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 "TObjArray.h" + +#include +#include +#include +#include + +int main(int argc, char* argv[]) { + + jungfrauStyle(); + gStyle->SetOptFit(11); + + /* + if (argc != 6) { + 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" << endl; + cout << "arg 4: specify pede file" << endl; + cout << "arg 5: specify data file" << endl; + cout << " " << endl; + exit(1); + } + */ //uncomment for SR + + 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); + } //uncomment for VH 210906 + + 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]; //uncomment for VH 210906 + + bool save_to_afs = false; + + int createHistoFile = 1; + + //char histoname[128]; + //char savename[128]; + char histoname[256]; // VH 210902 + char savename[256]; // VH 210902 + char plotfolder[256]; + char rootdatafolder[256]; + int filen = 352; + int pedefilen = 5; + int this_storagecell = stoi(argv[7]); + + // 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); + // /mnt/pcmoench_jungfrau_data/jungfrau_ana_sophie/Mxxx_CalibAna + //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna", module_str.c_str()); //uncomment for SR + sprintf(savename,"%s", anadata_loc.c_str()); //uncomment for VH 210906 + 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); + + TCanvas* c1 = new TCanvas("c1",""); + + /*****************************************************************************************************************/ + /* 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); + + 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() ); //VH: note, this adds a double slash in the filepath + 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 + + cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl; + + } // end if i creatHistoFile + + // We have saved the data into root-file by first calling CuFluo_analysis_sc_singlethread_data.cpp. Now, we read from that same root-file and perform the fitting. + + c1->cd(); + + //sprintf(savename,"/mnt/sls_det_storage/jungfrau_data1/jungfrau_ana_sophie/M%s_CalibAna/CuFluo_%s_file0to%d.root", module_str.c_str(), gain_str.c_str(), filen-1); + sprintf(savename,"%s/CuFluo_%s_sc%d_file0to%d.root", anadata_loc.c_str(), gain_str.c_str(), this_storagecell, filen-1); //uncomment for VH 210906 + TFile* comb_file = new TFile((const char *)(savename),"READ"); + + pixelMaskObject->initialisePixelMask(pixel_mask); + + int low_ADU_peak = 0; + int high_ADU_peak = 0; + if (gain_str == "HG0") { + + + low_ADU_peak = 700; + high_ADU_peak = 900; + if (isJF11) { + low_ADU_peak = 850; + high_ADU_peak = 1350; + + } + + } else if (gain_str == "G0") { + low_ADU_peak = 200; // was 250 + high_ADU_peak = 400; + } + + //declare + TH1F* fit_par3; + TH1F* fit_par4; + TH1F* fit_par5; + TH1F* fit_par6; + TH1F* fit_par7; + + TH2F* fit_par3_2d; + TH2F* fit_par4_2d; + TH2F* fit_par5_2d; + TH2F* fit_par6_2d; + TH2F* fit_par7_2d; + + TH1F* peak_fit_pos; + TH1F* peak_fit_poserr; + TH2F* peak_fit_pos_2d; + TH2F* peak_fit_poserr_2d; + + TH1F* noise_fit_pos; + TH1F* noise_fit_poserr; + TH2F* noise_fit_pos_2d; + TH2F* noise_fit_poserr_2d; + + TH1F* gain_fit; + TH1F* gain_fiterr; + TH2F* gain_fit_2d; + TH2F* gain_fiterr_2d; + + TH1F* gain_fit_isEdge; + TH1F* gain_fit_isInnerEdge; + TH1F* gain_fit_isDouble; + TH1F* gain_fit_isNextToDouble; + TH1F* gain_fit_isQuad; + TH1F* gain_fit_isBulk; + + TH2F* gain_ADUper1keV_2d; + TH2F* gainerr_ADUper1keV_2d; + + //initialize + Char_t *fithistoname = new Char_t[50]; + snprintf( fithistoname, 50, "fit_par3_sc%d", this_storagecell ); + fit_par3 = new TH1F( fithistoname, "", 100, 0, 50 ); + snprintf( fithistoname, 50, "fit_par4_sc%d", this_storagecell ); + fit_par4 = new TH1F( fithistoname, "", 100, 0, 500 ); + snprintf( fithistoname, 50, "fit_par5_sc%d", this_storagecell ); + fit_par5 = new TH1F( fithistoname, "", 100, 0, 0.5 ); + snprintf( fithistoname, 50, "fit_par6_sc%d", this_storagecell ); + fit_par6 = new TH1F( fithistoname, "", 100, 1.05, 1.3 ); // was 1.05, 1.25 + snprintf( fithistoname, 50, "fit_par7_sc%d", this_storagecell ); + fit_par7 = new TH1F( fithistoname, "", 100, 0, 0.4 ); + + snprintf( fithistoname, 50, "fit_par3_2d_sc%d", this_storagecell ); + fit_par3_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "fit_par4_2d_sc%d", this_storagecell ); + fit_par4_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "fit_par5_2d_sc%d", this_storagecell ); + fit_par5_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "fit_par6_2d_sc%d", this_storagecell ); + fit_par6_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "fit_par7_2d_sc%d", this_storagecell ); + fit_par7_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + + snprintf( fithistoname, 50, "peak_fit_pos_sc%d", this_storagecell ); + peak_fit_pos = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "peak_fit_poserr_sc%d", this_storagecell ); + peak_fit_poserr = new TH1F( fithistoname, "", 100, 0, 2 ); + snprintf( fithistoname, 50, "peak_fit_pos_2d_sc%d", this_storagecell ); + peak_fit_pos_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "peak_fit_poserr_2d_sc%d", this_storagecell ); + peak_fit_poserr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + + snprintf( fithistoname, 50, "noise_fit_pos_sc%d", this_storagecell ); + noise_fit_pos = new TH1F( fithistoname, "", 100, -10, 10 ); + snprintf( fithistoname, 50, "noise_fit_poserr_sc%d", this_storagecell ); + noise_fit_poserr = new TH1F( fithistoname, "", 100, 0, 0.15 ); // was 0, 0.1 + snprintf( fithistoname, 50, "noise_fit_pos_2d_sc%d", this_storagecell ); + noise_fit_pos_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "noise_fit_poserr_2d_sc%d", this_storagecell ); + noise_fit_poserr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + + snprintf( fithistoname, 50, "gain_fit_sc%d", this_storagecell ); + gain_fit = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fiterr_sc%d", this_storagecell ); + gain_fiterr = new TH1F( fithistoname, "", 100, 0, 2 ); + snprintf( fithistoname, 50, "gain_fit_2d_sc%d", this_storagecell ); + gain_fit_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "gain_fiterr_2d_sc%d", this_storagecell ); + gain_fiterr_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + + snprintf( fithistoname, 50, "gain_fit_isEdge_sc%d", this_storagecell ); + gain_fit_isEdge = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fit_isInnerEdge_sc%d", this_storagecell ); + gain_fit_isInnerEdge = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fit_isDouble_sc%d", this_storagecell ); + gain_fit_isDouble = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fit_isNextToDouble_sc%d", this_storagecell ); + gain_fit_isNextToDouble = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fit_isQuad_sc%d", this_storagecell ); + gain_fit_isQuad = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + snprintf( fithistoname, 50, "gain_fit_isBulk_sc%d", this_storagecell ); + gain_fit_isBulk = new TH1F( fithistoname, "", 100, low_ADU_peak, high_ADU_peak ); + + snprintf( fithistoname, 50, "gain_ADUper1keV_2d_sc%d", this_storagecell ); + gain_ADUper1keV_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + snprintf( fithistoname, 50, "gainerr_ADUper1keV_2d_sc%d", this_storagecell ); + gainerr_ADUper1keV_2d = new TH2F( fithistoname, "", NC, -0.5, NC-0.5, NR, -0.5, NR-0.5 ); + + int low_bin_noise = 101; + int high_bin_noise = 301; + int low_bin_peak = 0; + int high_bin_peak = 0; + if (gain_str == "HG0") { + low_bin_peak = 701; + high_bin_peak = 1200; + if (isJF11) { + low_bin_peak = 801; + high_bin_peak = 1451; + } + + } else if (gain_str == "G0") { + low_bin_peak = 301; + high_bin_peak = 651; + } + + for ( int j = 1; j < 9; ++j ) { + cout << "SC = " << this_storagecell << " slice " << j << endl; + sprintf( histoname, "adc2d_%d_sc%d", j, this_storagecell ); + 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); + + if ( proj->Integral( low_bin_noise, high_bin_noise ) != 0 && proj->Integral( low_bin_peak, high_bin_peak ) != 0 ) { + + // noise + TH1D *proj_noise = dynamic_cast(proj->Rebin(4,"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"); + + noise_fit_pos->Fill(fit->GetParameter(1)); + noise_fit_pos_2d->Fill(i%NC,i/NC,fit->GetParameter(1)); + noise_fit_poserr->Fill(fit->GetParError(1)); + noise_fit_poserr_2d->Fill(i%NC,i/NC,fit->GetParError(1)); + + // peak + TH1D *proj_peak = dynamic_cast(proj->Rebin(4,"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; + mypar[2] = proj_peak->GetBinCenter(proj_peak->GetMaximumBin()); + if (gain_str == "G0") { + mypar[3] = 16.; + } else if (gain_str == "HG0") { + mypar[3] = 29.; + } + mypar[4] = proj_peak->GetBinContent(proj_peak->GetMaximumBin()); + if (gain_str == "G0") { + mypar[5] = 0.17; + } else if (gain_str == "HG0") { + mypar[5] = 0.14; + } + mypar[6] = 1.12; + 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 ); + + fit_par3->Fill( mypar[3] ); + fit_par4->Fill( mypar[4] ); + fit_par5->Fill( mypar[5] ); + fit_par6->Fill( mypar[6] ); + fit_par7->Fill( mypar[7] ); + + fit_par3_2d->Fill( i%NC, i/NC, mypar[3] ); + fit_par4_2d->Fill( i%NC, i/NC, mypar[4] ); + fit_par5_2d->Fill( i%NC, i/NC, mypar[5] ); + fit_par6_2d->Fill( i%NC, i/NC, mypar[6] ); + fit_par7_2d->Fill( i%NC, i/NC, mypar[7] ); + + 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 >= 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"; + } + + proj_noise->Draw(); + c1->Update(); + proj_noise->GetXaxis()->SetTitle( "Pedestal corrected ADC [ADU]" ); + proj_noise->GetXaxis()->SetRangeUser( -100, 150 ); + 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); + sprintf( savename, "%s/noise_%s_sc%d_%d_%s_M%s.png", plotfolder, pixel_type.c_str(), this_storagecell, i, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + 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(); + erfc_Kb->Draw("same"); + erfc_Ka->Draw("same"); + gaus_Kb->Draw("same"); + gaus_Ka->Draw("same"); + fittedfun->Draw("same"); + c1->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.22); + st->SetX2NDC(0.62); + st->SetY1NDC(0.7); + st->SetY2NDC(0.94); + st->SetBorderSize(0); + st->SetTextSize(0.04); + sprintf( savename, "%s/peak_%s_sc%d_%d_%s_M%s.png", plotfolder, pixel_type.c_str(), this_storagecell, i, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + } + + // gain + gain_fit->Fill( mypar[2] - fit->GetParameter(1) ); + gain_fiterr->Fill( sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) ); + gain_fit_2d->Fill( i%NC, i/NC, mypar[2] - fit->GetParameter(1) ); + gain_fiterr_2d->Fill( i%NC, i/NC, sqrt( pow(emypar[2],2) + pow(fit->GetParError(1),2) ) ); + gain_ADUper1keV_2d->Fill( i%NC, i/NC, ( mypar[2] - fit->GetParameter(1) ) / 8.0 ); + gainerr_ADUper1keV_2d->Fill( i%NC, i/NC, sqrt( pow( emypar[2],2 ) + pow( fit->GetParError(1),2) ) / 8.0 ); + + if (isEdge(i)) { + gain_fit_isEdge->Fill( mypar[2] - fit->GetParameter(1) ); + } + if (isInnerEdge(i)) { + gain_fit_isInnerEdge->Fill( mypar[2] - fit->GetParameter(1) ); + } + if (isDouble(i)) { + gain_fit_isDouble->Fill( mypar[2] - fit->GetParameter(1) ); + } + if (isNextToDouble(i)) { + gain_fit_isNextToDouble->Fill( mypar[2] - fit->GetParameter(1) ); + } + if (isQuad(i)) { + gain_fit_isQuad->Fill( mypar[2] - fit->GetParameter(1) ); + } + if (isBulk(i)) { + gain_fit_isBulk->Fill( mypar[2] - fit->GetParameter(1) ); + } + + delete thiscalibration; + delete proj_noise; + delete proj_peak; + delete proj; + + } else { + pixel_mask[i] = false; + } //end if can be fitted + + } // end for loop pixel numbers (i) + + } // end for loop sensor slices (j) + + sprintf( savename, "%s/pixelmask_afterfit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + pixelMaskObject->plotPixelMask( pixel_mask, savename ); + + TPaveText *pave = new TPaveText( 0.86, 0.95, 0.91, 0.98, "blNDC" ); + pave->SetBorderSize(0); + pave->SetFillStyle(0); + pave->SetTextSize(0.06); + pave->SetTextAlign(32); + + mapcanvas->cd(); + + fit_par3_2d->GetXaxis()->SetTitle("Column"); + fit_par3_2d->GetYaxis()->SetTitle("Row"); + fit_par3_2d->GetYaxis()->SetTitleOffset(0.7); + fit_par3_2d->Draw("colz"); + fit_par3_2d->GetZaxis()->SetRangeUser(0,50); + sprintf( savename, "%s/fit_par3_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + fit_par3_2d->Delete(); + + fit_par4_2d->GetXaxis()->SetTitle("Column"); + fit_par4_2d->GetYaxis()->SetTitle("Row"); + fit_par4_2d->GetYaxis()->SetTitleOffset(0.7); + fit_par4_2d->Draw("colz"); + fit_par4_2d->GetZaxis()->SetRangeUser(0,500); + sprintf( savename, "%s/fit_par4_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + fit_par4_2d->Delete(); + + fit_par5_2d->GetXaxis()->SetTitle("Column"); + fit_par5_2d->GetYaxis()->SetTitle("Row"); + fit_par5_2d->GetYaxis()->SetTitleOffset(0.7); + fit_par5_2d->Draw("colz"); + fit_par5_2d->GetZaxis()->SetRangeUser(0,0.5); + sprintf( savename, "%s/fit_par5_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + fit_par5_2d->Delete(); + + fit_par6_2d->GetXaxis()->SetTitle("Column"); + fit_par6_2d->GetYaxis()->SetTitle("Row"); + fit_par6_2d->GetYaxis()->SetTitleOffset(0.7); + fit_par6_2d->Draw("colz"); + fit_par6_2d->GetZaxis()->SetRangeUser(1.0,1.3); //was (1.0,1.25) + sprintf( savename, "%s/fit_par6_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + fit_par6_2d->Delete(); + + fit_par7_2d->GetXaxis()->SetTitle("Column"); + fit_par7_2d->GetYaxis()->SetTitle("Row"); + fit_par7_2d->GetYaxis()->SetTitleOffset(0.7); + fit_par7_2d->Draw("colz"); + fit_par7_2d->GetZaxis()->SetRangeUser(0.,0.4); + sprintf( savename, "%s/fit_par7_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + 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( low_ADU_peak, high_ADU_peak ); + sprintf( savename, "%s/peak_fit_pos_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + peak_fit_pos_2d->Delete(); + + 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, "%s/peak_fit_poserr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + peak_fit_poserr_2d->Delete(); + + 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(-5,5); + sprintf( savename, "%s/noise_fit_pos_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + noise_fit_pos_2d->Delete(); + + 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"); + if (gain_str == "HG0") { + noise_fit_poserr_2d->GetZaxis()->SetRangeUser(0,0.1); + } else if (gain_str == "G0") { + noise_fit_poserr_2d->GetZaxis()->SetRangeUser(0,0.1); // was 0,0.05 + } + sprintf( savename, "%s/noise_fit_poserr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + noise_fit_poserr_2d->Delete(); + + gain_fit_2d->GetXaxis()->SetTitle("Column"); + gain_fit_2d->GetYaxis()->SetTitle("Row"); + gain_fit_2d->GetYaxis()->SetTitleOffset(0.7); + gain_fit_2d->Draw("colz"); + sprintf( savename, "%s [ADU/8 keV]", gain_str.c_str() ); + pave->AddText( (const char *)(savename) ); + pave->Draw(); + gain_fit_2d->GetZaxis()->SetRangeUser( low_ADU_peak, high_ADU_peak ); + sprintf( savename, "%s/gain_fit_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + gain_fit_2d->Delete(); + + 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, "%s/gain_fiterr_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + gain_fiterr_2d->Delete(); + + gain_ADUper1keV_2d->GetXaxis()->SetTitle("Column"); + gain_ADUper1keV_2d->GetYaxis()->SetTitle("Row"); + gain_ADUper1keV_2d->GetYaxis()->SetTitleOffset(0.7); + gain_ADUper1keV_2d->Draw("colz"); + if (gain_str == "HG0") { + gain_ADUper1keV_2d->GetZaxis()->SetRangeUser(80,120); + } else if (gain_str == "G0") { + gain_ADUper1keV_2d->GetZaxis()->SetRangeUser(25,45); //was 35,50 + } + sprintf( savename, "%s/gain_ADUper1keV_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + //gain_ADUper1keV_2d->Delete(); + + gainerr_ADUper1keV_2d->GetXaxis()->SetTitle("Column"); + gainerr_ADUper1keV_2d->GetYaxis()->SetTitle("Row"); + gainerr_ADUper1keV_2d->GetYaxis()->SetTitleOffset(0.7); + gainerr_ADUper1keV_2d->Draw("colz"); + if (gain_str == "HG0") { + gainerr_ADUper1keV_2d->GetZaxis()->SetRangeUser(0,0.5); + } else if (gain_str == "G0") { + gainerr_ADUper1keV_2d->GetZaxis()->SetRangeUser(0,0.25); + } + sprintf( savename, "%s/gainerr_ADUper1keV_2d_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + mapcanvas->SaveAs( (const char *)(savename) ); + + //gainerr_ADUper1keV_2d[this_storagecell]->Delete(); + + c1->cd(); + + fit_par3->GetXaxis()->SetTitle("Fit par 3"); + fit_par3->Draw(); + sprintf( savename, "%s/fit_par3_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + fit_par3->Delete(); + + fit_par4->GetXaxis()->SetTitle("Fit par 4"); + fit_par4->Draw(); + sprintf( savename, "%s/fit_par4_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + fit_par4->Delete(); + + fit_par5->GetXaxis()->SetTitle("Fit par 5"); + fit_par5->Draw(); + sprintf( savename, "%s/fit_par5_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + fit_par5->Delete(); + + fit_par6->GetXaxis()->SetTitle("Fit par 6"); + fit_par6->Draw(); + sprintf( savename, "%s/fit_par6_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + fit_par6->Delete(); + + fit_par7->GetXaxis()->SetTitle("Fit par 7"); + fit_par7->Draw(); + sprintf( savename, "%s/fit_par7_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + peak_fit_pos->GetXaxis()->SetTitle("Peak position [ADU]"); + peak_fit_pos->Draw(); + sprintf( savename, "%s/peak_fit_pos_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + peak_fit_pos->Delete(); + + peak_fit_poserr->GetXaxis()->SetTitle("Peak position uncert [ADU]"); + peak_fit_poserr->Draw(); + sprintf( savename, "%s/peak_fit_poserr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + peak_fit_poserr->Delete(); + + noise_fit_pos->GetXaxis()->SetTitle("Noise position [ADU]"); + noise_fit_pos->Draw(); + sprintf( savename, "%s/noise_fit_pos_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + noise_fit_pos->Delete(); + + noise_fit_poserr->GetXaxis()->SetTitle("Noise position uncert [ADU]"); + noise_fit_poserr->Draw(); + sprintf( savename, "%s/noise_fit_poserr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + noise_fit_poserr->Delete(); + + sprintf( savename, "Gain %s [ADU / 8 keV]", gain_str.c_str() ); + gain_fit->GetXaxis()->SetTitle( (const char *)(savename) ); + gain_fit->Draw(); + sprintf( savename, "%s/gain_fit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + //gain_fit->Delete(); + + gain_fit->GetXaxis()->SetRangeUser( low_ADU_peak+30, high_ADU_peak ); + gain_fit->Fit("gaus"); + gain_fit->Draw(); + c1->Update(); + TPaveText *pave2 = new TPaveText( 0.6, 0.8, 0.94, 0.94, "blNDC" ); + pave2->SetBorderSize(0); + pave2->SetFillStyle(0); + pave2->SetTextSize(0.04); + pave2->SetTextAlign(32); + TF1* gain_fit_gaus = gain_fit->GetFunction("gaus"); + sprintf( savename, "Mean %0.2f #pm %0.2f", gain_fit_gaus->GetParameter(1), gain_fit_gaus->GetParError(1) ); + pave2->AddText( (const char *)(savename) ); + sprintf( savename, "Sigma %0.2f #pm %0.2f", gain_fit_gaus->GetParameter(2), gain_fit_gaus->GetParError(2) ); + pave2->AddText( (const char *)(savename) ); + pave2->Draw(); + gain_fit->SetStats(kFALSE); + sprintf( savename, "%s/gain_fit_fit_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + gain_fit->Delete(); + + gain_fiterr->GetXaxis()->SetTitle("Gain uncert [ADU / 8 keV]"); + gain_fiterr->Draw(); + sprintf( savename, "%s/gain_fiterr_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + gain_fiterr->Delete(); + + gain_fit_isEdge->SetLineColor(kBlue); + gain_fit_isInnerEdge->SetLineColor(kCyan); + gain_fit_isDouble->SetLineColor(kGreen+2); + gain_fit_isNextToDouble->SetLineColor(kRed); + gain_fit_isQuad->SetLineColor(kOrange); + + gain_fit_isEdge->Scale( 1./gain_fit_isEdge->GetEntries() ); + gain_fit_isInnerEdge->Scale( 1./gain_fit_isInnerEdge->GetEntries() ); + gain_fit_isDouble->Scale( 1./gain_fit_isDouble->GetEntries() ); + gain_fit_isNextToDouble->Scale( 1./gain_fit_isNextToDouble->GetEntries() ); + gain_fit_isQuad->Scale( 1./gain_fit_isQuad->GetEntries() ); + gain_fit_isBulk->Scale( 1./gain_fit_isBulk->GetEntries() ); + + TLegend *leg = new TLegend( 0.62, 0.6, 0.93, 0.93 ); + leg->AddEntry( gain_fit_isBulk, "Normal", "l" ); + leg->AddEntry( gain_fit_isDouble, "Double", "l" ); + leg->AddEntry( gain_fit_isNextToDouble, "Next to D", "l" ); + leg->AddEntry( gain_fit_isEdge, "Edge", "l" ); + leg->AddEntry( gain_fit_isInnerEdge, "Inner E", "l" ); + + sprintf( savename, "Gain %s [ADU / 8 keV]", gain_str.c_str() ); + gain_fit_isDouble->GetXaxis()->SetTitle( (const char *)(savename) ); + gain_fit_isDouble->GetYaxis()->SetTitle("Normalised"); + gain_fit_isDouble->GetYaxis()->SetTitleOffset(1.3); + gain_fit_isDouble->SetMinimum(0.0); + gain_fit_isDouble->SetMaximum(0.16); + gain_fit_isDouble->Draw(); + gain_fit_isEdge->Draw("same"); + gain_fit_isInnerEdge->Draw("same"); + gain_fit_isNextToDouble->Draw("same"); + gain_fit_isBulk->Draw("same"); + leg->Draw("same"); + sprintf( savename, "%s/gain_fit_perType_sc%d_%s_M%s.png", plotfolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + c1->SaveAs( (const char *)(savename) ); + + gain_fit_isDouble->Delete(); + gain_fit_isEdge->Delete(); + gain_fit_isInnerEdge->Delete(); + gain_fit_isNextToDouble->Delete(); + gain_fit_isBulk->Delete(); + + sprintf( savename, "%s/CuFluo_gain_sc%d_%s_M%s.root", rootdatafolder, this_storagecell, gain_str.c_str(), module_str.c_str() ); + TFile* saved_file = new TFile( (const char *)(savename), "RECREATE" ); + gain_ADUper1keV_2d->Write(); + gainerr_ADUper1keV_2d->Write(); + saved_file->Close(); + + gain_ADUper1keV_2d->Delete(); + gainerr_ADUper1keV_2d->Delete(); + +} diff --git a/makefile b/makefile index 6682304..7fef141 100644 --- a/makefile +++ b/makefile @@ -3,7 +3,13 @@ CuFluo_analysis: CuFluo_analysis.cpp g++ -Wall -O3 -m64 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic CuFluo_analysis.cpp -o CuFluo_analysis CuFluo_analysis_sc: CuFluo_analysis_sc.cpp - g++ -Wall -O3 -m64 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic CuFluo_analysis_sc.cpp -o CuFluo_analysis_sc + g++ -Wall -O3 -m64 -std=c++11 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic CuFluo_analysis_sc.cpp -o CuFluo_analysis_sc + +CuFluo_analysis_sc_singlethread_data: CuFluo_analysis_sc_singlethread_data.cpp + g++ -Wall -O3 -m64 -std=c++11 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic CuFluo_analysis_sc_singlethread_data.cpp -o CuFluo_analysis_sc_singlethread_data + +CuFluo_analysis_sc_singlethread_fits: CuFluo_analysis_sc_singlethread_fits.cpp + g++ -Wall -O3 -m64 -std=c++11 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic CuFluo_analysis_sc_singlethread_fits.cpp -o CuFluo_analysis_sc_singlethread_fits DB_analysis: DB_analysis.cpp g++ -Wall -O3 -m64 -I$(ROOTSYS)/include -L$(ROOTSYS)/lib -lGui -lCore -lCint -lRIO -lNet -lHist -lGraf -lGraf3d -lGpad -lTree -lRint -lPostscript -lMatrix -lPhysics -lMathCore -lThread -pthread -lm -ldl -rdynamic DB_analysis.cpp -o DB_analysis