commit cecd8e0ab8d4a7d30a3e60115c96623a6823f789 Author: Sophie Redford Date: Mon Nov 28 14:28:10 2016 +0100 First commit for JUNGFRAU module calibration version 1. diff --git a/JFMC_CalibWriter.cpp b/JFMC_CalibWriter.cpp new file mode 100644 index 0000000..1939cb2 --- /dev/null +++ b/JFMC_CalibWriter.cpp @@ -0,0 +1,296 @@ +// 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* 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)); + } + + 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))); + } + + 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))); + } + + } + + 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)); + + 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()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->SetLogy(); + sprintf(savename,"plots/M%s/g0hist_log_M%s.png", module_str.c_str(), module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->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()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->SetLogy(); + sprintf(savename,"plots/M%s/g1hist_log_M%s.png", module_str.c_str(), module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->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()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->SetLogy(); + sprintf(savename,"plots/M%s/g2hist_log_M%s.png", module_str.c_str(), module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->SetLogy(0); + +} diff --git a/JFMC_CuFluoPeak.cpp b/JFMC_CuFluoPeak.cpp new file mode 100644 index 0000000..eaef147 --- /dev/null +++ b/JFMC_CuFluoPeak.cpp @@ -0,0 +1,275 @@ +// file to calculate pedestal correction of fluo data +// make correction and save spectrum per pixel + +#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 "TGraph.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "TFile.h" + + +int main(int argc, char* argv[]) { + jungfrauStyle(); + + if (argc != 5) { + cout << "Correct usage:" << endl; + cout << "arg 1: specify module number" << endl; + cout << "arg 2: specify data location" << endl; + cout << "arg 3: specify pede file" << endl; + cout << "arg 4: specify data file" << endl; + cout << "eg: ./JFMC_FluoCalibFit 006 /data_pool/Module_006_161116 10us_500Hz_QS_allgain_pede 10us_500Hz_QS_G0_Cu" << endl; + cout << " " << endl; + exit(1); + } + + string module_str = argv[1]; + string data_loc = argv[2]; + string pede_file = argv[3]; + string data_file = argv[4]; + + jungfrauFile *thisfile = new jungfrauFile(); + + jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask(); + static int pixel_mask [NCH]; + + jungfrauPedestal *pedestalObject = new jungfrauPedestal(); + pedestalObject->pedestalSetNFrames(1000); // using 1000 frames, rolling window + static uint16_t pedestals16_G0[NCH]; + + TH2F* pedestalsG0 = new TH2F("pedestalsG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedestalsG1 = new TH2F("pedestalsG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedestalsG2 = new TH2F("pedestalsG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH2F* pedeRMSG0 = new TH2F("pedeRMSG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedeRMSG1 = new TH2F("pedeRMSG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedeRMSG2 = new TH2F("pedeRMSG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH1D* adcpc_spec = new TH1D("adcpc_spec","",300,0,3000); + + TH2I* adc2d_1 = new TH2I("adc2d_1","",1000,-200-0.5,800-0.5,65536,(65536*0-0.5),(65536*1-0.5)); + TH2I* adc2d_2 = new TH2I("adc2d_2","",1000,-200-0.5,800-0.5,65536,(65536*1-0.5),(65536*2-0.5)); + TH2I* adc2d_3 = new TH2I("adc2d_3","",1000,-200-0.5,800-0.5,65536,(65536*2-0.5),(65536*3-0.5)); + TH2I* adc2d_4 = new TH2I("adc2d_4","",1000,-200-0.5,800-0.5,65536,(65536*3-0.5),(65536*4-0.5)); + TH2I* adc2d_5 = new TH2I("adc2d_5","",1000,-200-0.5,800-0.5,65536,(65536*4-0.5),(65536*5-0.5)); + TH2I* adc2d_6 = new TH2I("adc2d_6","",1000,-200-0.5,800-0.5,65536,(65536*5-0.5),(65536*6-0.5)); + TH2I* adc2d_7 = new TH2I("adc2d_7","",1000,-200-0.5,800-0.5,65536,(65536*6-0.5),(65536*7-0.5)); + TH2I* adc2d_8 = new TH2I("adc2d_8","",1000,-200-0.5,800-0.5,65536,(65536*7-0.5),(65536*8-0.5)); + + char savename[128]; + + 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); + + // open pede file + sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str()); + thisfile->open((char*)savename, 0); + + // calculate pixel mask + pixelMaskObject->initialisePixelMask(pixel_mask); + if (module_str == "006") { + pixelMaskObject->maskChip(6, pixel_mask); + pixelMaskObject->maskSupercolumn(4,1, pixel_mask); + pixelMaskObject->maskSupercolumn(4,2, pixel_mask); + } + + int nevents = 0; + while (thisfile->readNextFrame()) { + nevents++; + } + thisfile->rewind(); + cout << "read " << nevents << " events" << endl; + + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(0, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(1, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + for (int i = 0; i < 999; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(3, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + thisfile->rewind(); + + sprintf(savename,"plots/M%s/CuFluo/pixelmask.png", module_str.c_str()); + pixelMaskObject->plotPixelMask(pixel_mask,savename); + cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl; + + // caluclate pedestals + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG0->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG0->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalData((uint16_t*)(&pedestals16_G0)); + pedestalObject->pedestalClear(); + + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG1->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG1->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalClear(); + + for (int i = 0; i < 999; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG2->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG2->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalClear(); + + thisfile->close(); + + pedestalsG0->GetXaxis()->SetTitle("Column"); + pedestalsG0->GetYaxis()->SetTitle("Row"); + pedestalsG0->GetYaxis()->SetTitleOffset(0.7); + pedestalsG0->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeG0.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedestalsG1->GetXaxis()->SetTitle("Column"); + pedestalsG1->GetYaxis()->SetTitle("Row"); + pedestalsG1->GetYaxis()->SetTitleOffset(0.7); + pedestalsG1->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeG1.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedestalsG2->GetXaxis()->SetTitle("Column"); + pedestalsG2->GetYaxis()->SetTitle("Row"); + pedestalsG2->GetYaxis()->SetTitleOffset(0.7); + pedestalsG2->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeG2.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG0->GetXaxis()->SetTitle("Column"); + pedeRMSG0->GetYaxis()->SetTitle("Row"); + pedeRMSG0->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG0->GetZaxis()->SetRangeUser(0,30); + pedeRMSG0->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeRMSG0.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG1->GetXaxis()->SetTitle("Column"); + pedeRMSG1->GetYaxis()->SetTitle("Row"); + pedeRMSG1->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG1->GetZaxis()->SetRangeUser(0,100); + pedeRMSG1->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeRMSG1.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG2->GetXaxis()->SetTitle("Column"); + pedeRMSG2->GetYaxis()->SetTitle("Row"); + pedeRMSG2->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG2->GetZaxis()->SetRangeUser(0,300); + pedeRMSG2->Draw("colz"); + sprintf(savename,"plots/M%s/CuFluo/pedeRMSG2.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedestalsG0->Reset(); + pedestalsG1->Reset(); + pedestalsG2->Reset(); + pedeRMSG0->Reset(); + pedeRMSG1->Reset(); + pedeRMSG2->Reset(); + + for (int filen = 11; filen < 20; filen++) { + + // open data file + sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str()); + thisfile->open((char*)savename, filen); + + while (thisfile->readNextFrame()) { + + uint16_t* imagedptr = thisfile->getFrameDataHandle(); + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + + uint16_t gain = (imagedptr[i]&0xc000) >> 14; + + if (gain == 0) { + + int adcpc = (imagedptr[i]&0x3fff) - pedestals16_G0[i]; + adcpc_spec->Fill(adcpc); + + if (i < (65536*1)) { + adc2d_1->Fill(adcpc,i); + } else if (i < (65536*2)) { + adc2d_2->Fill(adcpc,i); + } else if (i < (65536*3)) { + adc2d_3->Fill(adcpc,i); + } else if (i < (65536*4)) { + adc2d_4->Fill(adcpc,i); + } else if (i < (65536*5)) { + adc2d_5->Fill(adcpc,i); + } else if (i < (65536*6)) { + adc2d_6->Fill(adcpc,i); + } else if (i < (65536*7)) { + adc2d_7->Fill(adcpc,i); + } else if (i < (65536*8)) { + adc2d_8->Fill(adcpc,i); + } + } + } + } + } + thisfile->close(); + + sprintf(savename,"%s_ana_M%s/Fluo_file%d.root", data_loc.c_str(), module_str.c_str(), filen); + TFile* saved_file = new TFile((const char *)(savename),"RECREATE"); + adc2d_1->Write(); + adc2d_2->Write(); + adc2d_3->Write(); + adc2d_4->Write(); + adc2d_5->Write(); + adc2d_6->Write(); + adc2d_7->Write(); + adc2d_8->Write(); + saved_file->Close(); + + adcpc_spec->GetXaxis()->SetTitle("Pedestal corrected ADC [ADU]"); + adcpc_spec->Draw(); + mapcanvas->SetLogy(); + sprintf(savename,"plots/M%s/CuFluo/adcpc_spec_file%d.png",module_str.c_str(), filen); + mapcanvas->SaveAs((const char *)(savename)); + mapcanvas->SetLogy(0); + + adc2d_1->Reset(); + adc2d_2->Reset(); + adc2d_3->Reset(); + adc2d_4->Reset(); + adc2d_5->Reset(); + adc2d_6->Reset(); + adc2d_7->Reset(); + adc2d_8->Reset(); + adcpc_spec->Reset(); + + } // end of file loop + +} diff --git a/JFMC_CuFluoPeakFit.cpp b/JFMC_CuFluoPeakFit.cpp new file mode 100644 index 0000000..322ec74 --- /dev/null +++ b/JFMC_CuFluoPeakFit.cpp @@ -0,0 +1,152 @@ +// 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 != 3) { + cout << "Correct usage:" << endl; + cout << "arg 1: specify module number" << endl; + cout << "arg 2: specify data location" << endl; + cout << "eg: ./JFMC_FluoCalibFit 006 /data_pool/Module_006_161116" << endl; + cout << " " << endl; + exit(1); + } + + string module_str = argv[1]; + string data_loc = argv[2]; + + char histoname[128]; + char savename[128]; + TCanvas* c1 = new TCanvas("c1",""); + c1->SetLogy(); + + sprintf(savename,"%s_ana_M%s/Fluo_comb.root", data_loc.c_str(), module_str.c_str()); + TFile* comb_file = new TFile((const char *)(savename),"READ"); + + TH1F* peak_fit_pos = new TH1F("peak_fit_pos","",100,200,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); + + 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(300,600)!=0) { + + proj->Rebin(4); + proj->GetXaxis()->SetRangeUser(100,400); + proj->SetStats(kTRUE); + + //Double_t mypar[8]; + Double_t mypar[6]; + mypar[0] = 0.0; + mypar[1] = 0.0; + mypar[2] = proj->GetBinCenter(proj->GetMaximumBin()); + //mypar[6] = 1.127; + //mypar[7] = 0.2; + mypar[3] = 19.; + mypar[4] = 60.; + mypar[5] = 0.18; + + //Double_t emypar[8]; + Double_t emypar[6]; + energyCalibration *thiscalibration = new energyCalibration(); + thiscalibration->setScanSign(1); // new + thiscalibration->setStartParametersKb(mypar); + thiscalibration->fixParameter(0,0.); // no background + thiscalibration->fixParameter(1,0.); // new + //TF1* fittedfun = thiscalibration->fitSpectrumKb(proj,mypar,emypar); + 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(); + 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_test%d.png",module_str.c_str(),i); + c1->SaveAs((const char *)(savename)); + } + + 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(200,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,3); + sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr_2d.png",module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + peak_fit_pos->GetXaxis()->SetTitle("Peak position [ADC]"); + peak_fit_pos->Draw(); + sprintf(savename,"plots/M%s/CuFluo/peak_fit_pos.png",module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + peak_fit_poserr->GetXaxis()->SetTitle("Peak position uncert [ADC]"); + peak_fit_poserr->Draw(); + sprintf(savename,"plots/M%s/CuFluo/peak_fit_poserr.png",module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + sprintf(savename,"data/M%s/CuF_peakpos_M%s.root", module_str.c_str(), module_str.c_str()); + TFile* saved_file = new TFile((const char *)(savename),"RECREATE"); + peak_fit_pos_2d->Write(); + peak_fit_poserr_2d->Write(); + saved_file->Close(); + +} diff --git a/JFMC_CurrentSourceScan.cpp b/JFMC_CurrentSourceScan.cpp new file mode 100644 index 0000000..1be4806 --- /dev/null +++ b/JFMC_CurrentSourceScan.cpp @@ -0,0 +1,742 @@ +// file to take current source scan files, analyse and plot them. + +#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/jungfrauPedestal.C" +#include "/afs/psi.ch/project/mythen/sophie/sls_detector_calibration/jungfrauPixelMask.C" + +#include "TGraph.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "TFile.h" + +#include +#include + +int main(int argc, char* argv[]) { + jungfrauStyle(); + + if (argc != 4) { + cout << "Correct usage:" << endl; + cout << "arg 1: specify module number" << endl; + cout << "arg 2: specify data location" << endl; + cout << "arg 3: specify data file" << endl; + cout << "eg: ./JFMC_CurrentSourceScan 006 ..." << endl; + cout << " " << endl; + exit(1); + } + + string module_str = argv[1]; + string data_loc = argv[2]; + string data_file = argv[3]; + + jungfrauFile *thisfile = new jungfrauFile(); + char savename[128]; + sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str()); + thisfile->open((char*)savename, 0); + + TH2F *gainmap_all = new TH2F("gainmap_all","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *adcmap_all = new TH2F("adcmap_all","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *gainmap_sel = new TH2F("gainmap_sel","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *adcmap_sel = new TH2F("adcmap_sel","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH1F *gain_avg = new TH1F("gain_avg","",100,-1,4); + TH2F *gainmap_avg = new TH2F("gainmap_avg","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH1F *adc_avg_g1 = new TH1F("adc_avg_g1","",100,0,17000); + TH2F *adcmap_avg_g1 = new TH2F("adcmap_avg_g1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH1F *adc_avg_g2 = new TH1F("adc_avg_g2","",100,0,17000); + TH2F *adcmap_avg_g2 = new TH2F("adcmap_avg_g2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + jungfrauPedestal *pedestalObject = new jungfrauPedestal(); + pedestalObject->pedestalSetNFrames(640); + + jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask(); + static int pixel_mask [NCH]; + pixelMaskObject->initialisePixelMask(pixel_mask); + if (module_str == "032") { + pixelMaskObject->maskChip(2,pixel_mask); + pixelMaskObject->maskChip(3,pixel_mask); + pixelMaskObject->maskChip(4,pixel_mask); + } else if (module_str == "006") { + pixelMaskObject->maskChip(6, pixel_mask); + pixelMaskObject->maskSupercolumn(4,1, pixel_mask); + pixelMaskObject->maskSupercolumn(4,2, pixel_mask); + } + + int frame_counter = 0; + + static TH1I *adc_histos_g1 [NCH]; + static TH1I *adc_histos_g2 [NCH]; + static TH1I *gain_histos [NCH]; + + for(int i = 0; i < NCH; i++) { + ostringstream histogramNameStreamg1; + histogramNameStreamg1 << "adc_histos_g1_" << i; + adc_histos_g1[i] = new TH1I(histogramNameStreamg1.str().c_str(),"",100,0,17000); + ostringstream histogramNameStreamg2; + histogramNameStreamg2 << "adc_histos_g2_" << i; + adc_histos_g2[i] = new TH1I(histogramNameStreamg2.str().c_str(),"",100,0,17000); + ostringstream histogramNameStream2; + histogramNameStream2 << "gain_histos_" << i; + gain_histos[i] = new TH1I(histogramNameStream2.str().c_str(),"",100,-1,4); + } + + double filter[35]; + for (int i = 0; i < 9; i++) { + filter[i] = 0.05+(i*0.025); + } + for (int i = 0; i < 9; i++) { + filter[i+9] = 0.5+(i*0.25); + } + for (int i = 0; i < 17; i++) { + filter[i+9+9] = 5.+(i*2.5); + } + + static double adcmeans_g1[NCH][35]; + static double adcmeans_g2[NCH][35]; + static double adcmeanerrs_g1[NCH][35]; + static double adcmeanerrs_g2[NCH][35]; + + float pede_ene[1]; + pede_ene[0] = 0.01; + + float pede_val1[NCH]; + float pede_val2[NCH]; + + 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",""); + + for (int j = 0; j < 640; j++) { + frame_counter++; + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + cout << pedestalObject->pedestalOfChannel(0) << endl; + cout << pedestalObject->rmsOfChannel(0) << endl; + pedestalObject->pedestalClear(); + + for (int j = 0; j < 640; j++) { + frame_counter++; + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + pede_val1[i] = pedestalObject->pedestalOfChannel(i); + } + cout << pedestalObject->pedestalOfChannel(0) << endl; + cout << pedestalObject->rmsOfChannel(0) << endl; + pedestalObject->pedestalClear(); + + for (int j = 0; j < 640; j++) { + frame_counter++; + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + pede_val2[i] = pedestalObject->pedestalOfChannel(i); + } + cout << pedestalObject->pedestalOfChannel(0) << endl; + cout << pedestalObject->rmsOfChannel(0) << endl; + pedestalObject->pedestalClear(); + + int j = 0; + int f = 1; + int f_max = 0; + if (module_str == "032" || module_str == "021" || module_str == "022") { + f_max = 3; + } else if (module_str == "006" || module_str == "008") { + f_max = 4; + } + + while (f < f_max) { + while (thisfile->readNextFrame()) { + frame_counter++; + j++; + + uint16_t* imagedptr = thisfile->getFrameDataHandle(); + for (int i = 0; i < NCH; i++) { + + if ((j-1) < 10 || (j > 8075 && j < 8085)) { + if (pixel_mask[i] == 1) { + uint16_t adc = imagedptr[i]&0x3fff; + uint16_t gain = (imagedptr[i]&0xc000) >> 14; + adcmap_all->Fill(i%NC,i/NC,adc); + gainmap_all->Fill(i%NC,i/NC,gain); + } + } + + if (((i/NC <= 255) && (i%64 == (j-1)%64)) || + ((i/NC >= 256) &&((-1*(i-524287))%64 == (j-1)%64))) { + + if (pixel_mask[i] == 1) { + uint16_t adc = imagedptr[i]&0x3fff; + uint16_t gain = (imagedptr[i]&0xc000) >> 14; + + gain_histos[i]->Fill(gain); + if (gain == 1) { + adc_histos_g1[i]->Fill(adc); + } else if (gain == 3) { + adc_histos_g2[i]->Fill(adc); + } + + if ((j-1) < 10 || (j > 8075 && j < 8085)) { + adcmap_sel->Fill(i%NC,i/NC,adc); + gainmap_sel->Fill(i%NC,i/NC,gain); + } + + } + } + } + if ((j-1) < 10 || (j > 8075 && j < 8085)) { + mapcanvas->cd(); + adcmap_sel->GetXaxis()->SetTitle("Column"); + adcmap_sel->GetYaxis()->SetTitle("Row"); + adcmap_sel->GetYaxis()->SetTitleOffset(0.7); + adcmap_sel->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/adcmap_%i_sel.png", module_str.c_str(), j-1); + mapcanvas->SaveAs((const char *)(savename)); + adcmap_all->GetXaxis()->SetTitle("Column"); + adcmap_all->GetYaxis()->SetTitle("Row"); + adcmap_all->GetYaxis()->SetTitleOffset(0.7); + adcmap_all->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/adcmap_%i_all.png", module_str.c_str(), j-1); + mapcanvas->SaveAs((const char *)(savename)); + + gainmap_sel->GetXaxis()->SetTitle("Column"); + gainmap_sel->GetYaxis()->SetTitle("Row"); + gainmap_sel->GetYaxis()->SetTitleOffset(0.7); + gainmap_sel->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/gainmap_%i_sel.png", module_str.c_str(), j-1); + mapcanvas->SaveAs((const char *)(savename)); + gainmap_all->GetXaxis()->SetTitle("Column"); + gainmap_all->GetYaxis()->SetTitle("Row"); + gainmap_all->GetYaxis()->SetTitleOffset(0.7); + gainmap_all->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/gainmap_%i_all.png", module_str.c_str(), j-1); + mapcanvas->SaveAs((const char *)(savename)); + + adcmap_sel->Reset(); + adcmap_all->Reset(); + gainmap_sel->Reset(); + gainmap_all->Reset(); + } + + // catch the last frame (bad ctrl-c exit) + int last_frame = 0; + if (module_str == "032" || module_str == "021" || module_str == "022") { + last_frame = 19199; + } else if (module_str == "006" || module_str == "008") { + last_frame = 24319; + } + + if (j%640 == 0 || j == last_frame) { + if (j == last_frame) { + j = last_frame+1; + } + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + adcmeans_g1[i][j/640-1] = adc_histos_g1[i]->GetMean(); + adcmeans_g2[i][j/640-1] = adc_histos_g2[i]->GetMean(); + adcmeanerrs_g1[i][j/640-1] = adc_histos_g1[i]->GetMeanError(); + adcmeanerrs_g2[i][j/640-1] = adc_histos_g2[i]->GetMeanError(); + + if (gain_histos[i]->GetEntries() > 0) { + gain_avg->Fill(gain_histos[i]->GetMean()); + gainmap_avg->Fill(i%NC,i/NC,gain_histos[i]->GetMean()); + } + gain_histos[i]->Reset(); + if (adc_histos_g1[i]->GetEntries() > 0) { + adc_avg_g1->Fill(adc_histos_g1[i]->GetMean()); + adcmap_avg_g1->Fill(i%NC,i/NC,adc_histos_g1[i]->GetMean()); + } + adc_histos_g1[i]->Reset(); + if (adc_histos_g2[i]->GetEntries() > 0) { + adc_avg_g2->Fill(adc_histos_g2[i]->GetMean()); + adcmap_avg_g2->Fill(i%NC,i/NC,adc_histos_g2[i]->GetMean()); + } + adc_histos_g2[i]->Reset(); + } + } + c1->cd(); + gain_avg->GetXaxis()->SetTitle("Average gain"); + gain_avg->GetYaxis()->SetTitle("Pixels"); + gain_avg->GetXaxis()->SetTitleOffset(1.1); + gain_avg->GetYaxis()->SetTitleOffset(1.5); + gain_avg->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/perStep/gain_avg_%i.png", module_str.c_str(), j/640); + c1->SaveAs((const char *)(savename)); + mapcanvas->cd(); + gainmap_avg->GetXaxis()->SetTitle("Column"); + gainmap_avg->GetYaxis()->SetTitle("Row"); + gainmap_avg->GetYaxis()->SetTitleOffset(0.7); + gainmap_avg->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/perStep/gainmap_avg_%i.png", module_str.c_str(), j/640); + mapcanvas->SaveAs((const char *)(savename)); + + c1->cd(); + adc_avg_g1->GetXaxis()->SetTitle("Average G1 ADC"); + adc_avg_g1->GetYaxis()->SetTitle("Pixels"); + adc_avg_g1->GetXaxis()->SetTitleOffset(1.1); + adc_avg_g1->GetYaxis()->SetTitleOffset(1.5); + adc_avg_g1->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/perStep/adc_avg_g1_%i.png", module_str.c_str(), j/640); + c1->SaveAs((const char *)(savename)); + mapcanvas->cd(); + adcmap_avg_g1->GetXaxis()->SetTitle("Column"); + adcmap_avg_g1->GetYaxis()->SetTitle("Row"); + adcmap_avg_g1->GetYaxis()->SetTitleOffset(0.7); + adcmap_avg_g1->Draw("colz"); + mapcanvas->Update(); + sprintf(savename,"plots/M%s/CurrentSource/perStep/adcmap_avg_g1_%i.png", module_str.c_str(), j/640); + mapcanvas->SaveAs((const char *)(savename)); + + c1->cd(); + adc_avg_g2->GetXaxis()->SetTitle("Average G2 ADC"); + adc_avg_g2->GetYaxis()->SetTitle("Pixels"); + adc_avg_g2->GetXaxis()->SetTitleOffset(1.1); + adc_avg_g2->GetYaxis()->SetTitleOffset(1.5); + adc_avg_g2->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/perStep/adc_avg_g2_%i.png", module_str.c_str(), j/640); + c1->SaveAs((const char *)(savename)); + mapcanvas->cd(); + adcmap_avg_g2->GetXaxis()->SetTitle("Column"); + adcmap_avg_g2->GetYaxis()->SetTitle("Row"); + adcmap_avg_g2->GetYaxis()->SetTitleOffset(0.7); + adcmap_avg_g2->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/perStep/adcmap_avg_g2_%i.png", module_str.c_str(), j/640); + mapcanvas->SaveAs((const char *)(savename)); + gain_avg->Reset(); + gainmap_avg->Reset(); + adc_avg_g1->Reset(); + adcmap_avg_g1->Reset(); + adc_avg_g2->Reset(); + adcmap_avg_g2->Reset(); + } + + } // end of while + + thisfile->close(); + sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), data_file.c_str()); + thisfile->open((char*)savename, f); + f++; + + } // end of file while + thisfile->close(); + + for (int i = 0; i < NCH; i++) { + delete gain_histos[i]; + delete adc_histos_g1[i]; + delete adc_histos_g2[i]; + } + + cout << "frame_counter " << frame_counter << endl; + + TH1F *g1_avg; + if (module_str == "032") { + g1_avg = new TH1F("g1_avg","",100,-12000,-4000); + } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008") { + g1_avg = new TH1F("g1_avg","",100,-3000,0); + } + TH1F *g2_avg; + if (module_str == "032") { + g2_avg = new TH1F("g2_avg","",100,-1000,-200); + } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008") { + g2_avg = new TH1F("g2_avg","",100,-300,0); + } + + TH2F *g1map_avg = new TH2F("g1map_avg","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *g2map_avg = new TH2F("g2map_avg","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH1F *g1overg2 = new TH1F("g1overg2","",100,10,15); + TH1F *g1overg2er = new TH1F("g1overg2er","",100,0,2); + + TH2F *g1overg2_map = new TH2F("g1overg2_map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *g1overg2er_map = new TH2F("g1overg2er_map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH1F *g1_fit_er = new TH1F("g1_fit_er","",100,0,50); + TH1F *g2_fit_er = new TH1F("g2_fit_er","",100,0,10); + + TH2F *g1map_fit_er = new TH2F("g1map_fit_er","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *g2map_fit_er = new TH2F("g2map_fit_er","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH1F *diffFromLinG1 = new TH1F("diffFromLinG1","",100,-0.05,0.05); + TH1F *diffFromLinG2 = new TH1F("diffFromLinG2","",100,-0.5,0.5); + + TH2F *diffFromLinG1map = new TH2F("diffFromLinG1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F *diffFromLinG2map = new TH2F("diffFromLinG2map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + vector r0_adc; + vector r0_filter; + vector r0_adcerr; + vector r0_ferr; + + vector r1_adc; + vector r1_filter; + vector r1_adcerr; + vector r1_ferr; + + vector r3_adc; + vector r3_filter; + vector r3_adcerr; + vector r3_ferr; + + // filter out points at zero + int m_max = 0; + if (module_str == "032" || module_str == "021" || module_str == "022") { + m_max = 27; + } else if (module_str == "006" || module_str == "008") { + m_max = 35; + } + + for (int m = 0; m < m_max; m++) { + if (adcmeans_g1[i][m] != 0) { + r1_adc.push_back(adcmeans_g1[i][m]); + r1_filter.push_back(filter[m]); + r1_adcerr.push_back(adcmeanerrs_g1[i][m]); + r1_ferr.push_back(0.); + } + if (adcmeans_g2[i][m] != 0) { + r3_adc.push_back(adcmeans_g2[i][m]); + r3_filter.push_back(filter[m]); + r3_adcerr.push_back(adcmeanerrs_g2[i][m]); + r3_ferr.push_back(0.); + } + } + + float this_pede_val1[1] = {pede_val1[i]}; + float this_pede_val2[1] = {pede_val2[i]}; + + TGraph *grap_p1 = new TGraph(1,pede_ene,this_pede_val1); + TGraph *grap_p2 = new TGraph(1,pede_ene,this_pede_val2); + + TGraphErrors *grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_ferr[0]),&(r1_adcerr[0])); + TGraphErrors *grap_g2 = new TGraphErrors(r3_adc.size(),&(r3_filter[0]),&(r3_adc[0]),&(r3_ferr[0]),&(r3_adcerr[0])); + + TF1 *fit2; + TF1 *fit3; + + if (r1_adc.size() > 1) { + double rangemin2 = 0; + double rangemax2 = 0; + if (module_str == "032") { + rangemin2 = 0.75; + rangemax2 = *max_element(r1_filter.begin(),r1_filter.end()); + } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008") { + rangemin2 = 0.075; + rangemax2 = 2.5; + } + fit2 = new TF1("fit2","[0]+[1]*x",rangemin2,rangemax2); + fit2->SetParameter(0, 10000.); + fit2->SetParameter(1, -0.1); + fit2->SetLineColor(kGreen+2); + grap_g1->Fit(fit2,"QR",""); + + g1_avg->Fill(fit2->GetParameter(1)); + g1map_avg->Fill(i%NC,i/NC,fit2->GetParameter(1)); + g1_fit_er->Fill(fit2->GetParError(1)); + g1map_fit_er->Fill(i%NC,i/NC,fit2->GetParError(1)); + + double g1residsum = 0.; + double g1residn = 0.; + for (size_t it = 0; it < r1_adc.size(); it++) { + if (r1_filter[it] >= rangemin2) { + diffFromLinG1->Fill((r1_adc[it] - fit2->Eval(r1_filter[it])) / r1_adc[it]); + g1residsum = g1residsum + ((r1_adc[it] - fit2->Eval(r1_filter[it])) / r1_adc[it]); + g1residn = g1residn + 1.; + } + } + diffFromLinG1map->Fill(i%NC,i/NC,g1residsum/g1residn); + } + + if (r3_adc.size() > 1) { + double rangemin3 = 0; + double rangemax3 = 0; + if (module_str == "032" || module_str == "022" || module_str == "021") { + rangemin3 = 5; + rangemax3 = *max_element(r3_filter.begin(),r3_filter.end()); + } else if (module_str == "006" || module_str == "008") { + rangemin3 = 7.5; + rangemax3 = *max_element(r3_filter.begin(),r3_filter.end()); + } + fit3 = new TF1("fit3","[0]+[1]*x",rangemin3, rangemax3); + fit3->SetParameter(0, 10000.); + fit3->SetParameter(1, -0.01); + fit3->SetLineColor(kRed); + grap_g2->Fit(fit3,"QR+",""); + + g2_avg->Fill(fit3->GetParameter(1)); + g2map_avg->Fill(i%NC,i/NC,fit3->GetParameter(1)); + g2_fit_er->Fill(fit3->GetParError(1)); + g2map_fit_er->Fill(i%NC,i/NC,fit3->GetParError(1)); + + double g2residsum = 0.; + double g2residn = 0.; + for (size_t it = 0; it < r3_adc.size(); it++) { + diffFromLinG2->Fill((r3_adc[it] - fit3->Eval(r3_filter[it])) / r3_adc[it]); + g2residsum = g2residsum + ((r3_adc[it] - fit3->Eval(r3_filter[it])) / r3_adc[it]); + g2residn = g2residn + 1.; + } + diffFromLinG2map->Fill(i%NC,i/NC,g2residsum/g2residn); + } + + if (r1_adc.size() > 1 && r3_adc.size() > 1) { + g1overg2->Fill(fit2->GetParameter(1) / fit3->GetParameter(1)); + g1overg2er->Fill(abs(fit2->GetParameter(1)/fit3->GetParameter(1))*sqrt(pow((fit2->GetParError(1)/fit2->GetParameter(1)),2) + pow((fit3->GetParError(1)/fit3->GetParameter(1)),2))); + g1overg2_map->Fill(i%NC,i/NC,fit2->GetParameter(1) / fit3->GetParameter(1)); + g1overg2er_map->Fill(i%NC,i/NC,abs(fit2->GetParameter(1)/fit3->GetParameter(1))*sqrt(pow((fit2->GetParError(1)/fit2->GetParameter(1)),2) + pow((fit3->GetParError(1)/fit3->GetParameter(1)),2))); + } + + if (i > 100000 && i < 100010) { + c1->cd(); + grap_p1->SetMarkerColor(kGreen+2); + grap_p2->SetMarkerColor(kRed); + grap_p1->SetMarkerStyle(20); + grap_p2->SetMarkerStyle(20); + grap_g1->SetMarkerStyle(20); + grap_g2->SetMarkerStyle(20); + grap_g1->SetMarkerColor(kGreen+2); + grap_g2->SetMarkerColor(kRed); + grap_g1->SetLineColor(kGreen+2); + grap_g2->SetLineColor(kRed); + grap_g1->SetMinimum(0); + grap_g1->SetMaximum(16000); + grap_g2->SetMinimum(0); + grap_g2->SetMaximum(16000); + grap_g1->GetXaxis()->SetTitle("Integration time [#mus]"); + grap_g1->GetYaxis()->SetTitle("ADC"); + grap_g1->GetXaxis()->SetTitleOffset(1.1); + grap_g1->GetYaxis()->SetTitleOffset(1.5); + grap_g2->GetXaxis()->SetTitle("Integration time [#mus]"); + grap_g2->GetYaxis()->SetTitle("ADC"); + grap_g2->GetXaxis()->SetTitleOffset(1.1); + grap_g2->GetYaxis()->SetTitleOffset(1.5); + + TF1 *fit2_e; + TF1 *fit3_e; + + if (r1_adc.size() > 1) { + fit2_e = new TF1("fit2_e","[0]+[1]*x",0.009,1.1E3); + fit2_e->SetParameters(fit2->GetParameter(0),fit2->GetParameter(1)); + fit2_e->SetLineColor(kGreen+2); + fit2_e->SetLineStyle(2); + grap_g1->Draw("AP"); + grap_p1->Draw("P"); + fit2_e->Draw("same"); + sprintf(savename,"plots/M%s/CurrentSource/pixel%i_g1.png", module_str.c_str(), i); + c1->SaveAs((const char *)(savename)); + } + if (r3_adc.size() > 1) { + fit3_e = new TF1("fit3_e","[0]+[1]*x",0.009,1.1E3); + fit3_e->SetParameters(fit3->GetParameter(0),fit3->GetParameter(1)); + fit3_e->SetLineColor(kRed); + fit3_e->SetLineStyle(2); + grap_g2->GetXaxis()->SetLimits(0,50); + grap_g2->Draw("AP"); + grap_p2->Draw("P"); + fit3_e->Draw("same"); + sprintf(savename,"plots/M%s/CurrentSource/pixel%i_g2.png", module_str.c_str(), i); + c1->SaveAs((const char *)(savename)); + } + + grap_g1->Draw("AP"); + grap_g1->GetXaxis()->SetLimits(0.009,1.1E2); + + fit2_e->Draw("same"); + if (r3_adc.size() > 1) { + grap_g2->Draw("P"); + fit3_e->Draw("same"); + } + grap_p1->Draw("P"); + grap_p2->Draw("P"); + c1->SetLogx(); + sprintf(savename,"plots/M%s/CurrentSource/pixel%i.png", module_str.c_str(), i); + c1->SaveAs((const char *)(savename)); + c1->SetLogx(0); + } + + delete grap_g1; + delete grap_g2; + if (r1_adc.size() > 1) { + delete fit2; + } + if (r3_adc.size() > 1) { + delete fit3; + } + } + } + + c1->cd(); + + g1_avg->GetXaxis()->SetTitle("G1 gradient"); + g1_avg->GetYaxis()->SetTitle("Pixels"); + g1_avg->GetXaxis()->SetTitleOffset(1.1); + g1_avg->GetYaxis()->SetTitleOffset(1.5); + g1_avg->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g1_avg.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + g2_avg->GetXaxis()->SetTitle("G2 gradient"); + g2_avg->GetYaxis()->SetTitle("Pixels"); + g2_avg->GetXaxis()->SetTitleOffset(1.1); + g2_avg->GetYaxis()->SetTitleOffset(1.5); + g2_avg->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g2_avg.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + g1overg2->GetXaxis()->SetTitle("G1/G2 gradient"); + g1overg2->GetYaxis()->SetTitle("Pixels"); + g1overg2->GetXaxis()->SetTitleOffset(1.1); + g1overg2->GetYaxis()->SetTitleOffset(1.5); + g1overg2->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g1overg2.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + g1overg2er->GetXaxis()->SetTitle("G1/G2 uncertainty"); + g1overg2er->GetYaxis()->SetTitle("Pixels"); + g1overg2er->GetXaxis()->SetTitleOffset(1.1); + g1overg2er->GetYaxis()->SetTitleOffset(1.5); + g1overg2er->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g1overg2er.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + g1_fit_er->GetXaxis()->SetTitle("Fitted G1 uncert"); + g1_fit_er->GetYaxis()->SetTitle("Pixels"); + g1_fit_er->GetXaxis()->SetTitleOffset(1.1); + g1_fit_er->GetYaxis()->SetTitleOffset(1.5); + g1_fit_er->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g1er.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + g2_fit_er->GetXaxis()->SetTitle("Fitted G2 uncert"); + g2_fit_er->GetYaxis()->SetTitle("Pixels"); + g2_fit_er->GetXaxis()->SetTitleOffset(1.1); + g2_fit_er->GetYaxis()->SetTitleOffset(1.5); + g2_fit_er->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_g2er.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + diffFromLinG1->GetXaxis()->SetTitle("Normalised residual G1"); + diffFromLinG1->GetYaxis()->SetTitle(""); + diffFromLinG1->GetXaxis()->SetTitleOffset(1.1); + diffFromLinG1->GetYaxis()->SetTitleOffset(1.5); + diffFromLinG1->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_diffFromLinG1.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + diffFromLinG2->GetXaxis()->SetTitle("Normalised residual G2"); + diffFromLinG2->GetYaxis()->SetTitle(""); + diffFromLinG2->GetXaxis()->SetTitleOffset(1.1); + diffFromLinG2->GetYaxis()->SetTitleOffset(1.5); + diffFromLinG2->Draw(); + sprintf(savename,"plots/M%s/CurrentSource/module_diffFromLinG2.png", module_str.c_str()); + c1->SaveAs((const char *)(savename)); + + cout << "G1 underflow " << diffFromLinG1->GetBinContent(0) << endl; + cout << "G2 underflow " << diffFromLinG2->GetBinContent(0) << endl; + cout << "G1 overflow " << diffFromLinG1->GetBinContent(101) << endl; + cout << "G2 overflow " << diffFromLinG2->GetBinContent(101) << endl; + + mapcanvas->cd(); + + g1map_avg->GetXaxis()->SetTitle("Column"); + g1map_avg->GetYaxis()->SetTitle("Row"); + g1map_avg->GetYaxis()->SetTitleOffset(0.7); + if (module_str == "032") { + g1map_avg->GetZaxis()->SetRangeUser(-12000,-4000); + } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008") { + g1map_avg->GetZaxis()->SetRangeUser(-3000,0); + } + g1map_avg->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/module_g1map_avg.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g2map_avg->GetXaxis()->SetTitle("Column"); + g2map_avg->GetYaxis()->SetTitle("Row"); + g2map_avg->GetYaxis()->SetTitleOffset(0.7); + if (module_str == "032") { + g2map_avg->GetZaxis()->SetRangeUser(-800,-400); + } else if (module_str == "022" || module_str == "021" || module_str == "006" || module_str == "008") { + g2map_avg->GetZaxis()->SetRangeUser(-300,0); + } + g2map_avg->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/module_g2map_avg.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g1overg2_map->GetXaxis()->SetTitle("Column"); + g1overg2_map->GetYaxis()->SetTitle("Row"); + g1overg2_map->GetYaxis()->SetTitleOffset(0.7); + g1overg2_map->GetZaxis()->SetRangeUser(10,15); + g1overg2_map->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/module_g1overg2_map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g1overg2_map->GetXaxis()->SetRangeUser(0,256); + g1overg2_map->GetYaxis()->SetRangeUser(240,280); + g1overg2_map->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/module_g1overg2_map_zoom.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g1overg2er_map->GetXaxis()->SetTitle("Column"); + g1overg2er_map->GetYaxis()->SetTitle("Row"); + g1overg2er_map->GetYaxis()->SetTitleOffset(0.7); + g1overg2er_map->GetZaxis()->SetRangeUser(0,2); + g1overg2er_map->Draw("colz"); + sprintf(savename,"plots/M%s/CurrentSource/module_g1overg2er_map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g1map_fit_er->GetXaxis()->SetTitle("Column"); + g1map_fit_er->GetYaxis()->SetTitle("Row"); + g1map_fit_er->GetYaxis()->SetTitleOffset(0.7); + g1map_fit_er->Draw("colz"); + g1map_fit_er->GetZaxis()->SetRangeUser(0,50); + /* + if (module_str == "032") { + g1map_fit_er->GetZaxis()->SetRangeUser(0,100); + } else if (module_str == "006" || module_str == "008" || module_str == "022" || module_str == "021") { + g1map_fit_er->GetZaxis()->SetRangeUser(0,2000); + } + */ + sprintf(savename,"plots/M%s/CurrentSource/module_g1map_er.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g2map_fit_er->GetXaxis()->SetTitle("Column"); + g2map_fit_er->GetYaxis()->SetTitle("Row"); + g2map_fit_er->GetYaxis()->SetTitleOffset(0.7); + g2map_fit_er->Draw("colz"); + g2map_fit_er->GetZaxis()->SetRangeUser(0,10); + sprintf(savename,"plots/M%s/CurrentSource/module_g2map_er.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + diffFromLinG1map->GetXaxis()->SetTitle("Column"); + diffFromLinG1map->GetYaxis()->SetTitle("Row"); + diffFromLinG1map->GetYaxis()->SetTitleOffset(0.7); + diffFromLinG1map->Draw("colz"); + diffFromLinG1map->GetZaxis()->SetRangeUser(-0.01,0.01); + sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG1map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + diffFromLinG2map->GetXaxis()->SetTitle("Column"); + diffFromLinG2map->GetYaxis()->SetTitle("Row"); + diffFromLinG2map->GetYaxis()->SetTitleOffset(0.7); + diffFromLinG2map->Draw("colz"); + diffFromLinG2map->GetZaxis()->SetRangeUser(-0.1,0.1); + sprintf(savename,"plots/M%s/CurrentSource/diffFromLinG2map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + sprintf(savename,"data/M%s/CS_ratio_M%s.root", module_str.c_str(), module_str.c_str()); + TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE"); + g1overg2_map->Write(); + g1overg2er_map->Write(); + saved_file2->Close(); + +} diff --git a/JFMC_DirectBeamScan.cpp b/JFMC_DirectBeamScan.cpp new file mode 100644 index 0000000..aff311d --- /dev/null +++ b/JFMC_DirectBeamScan.cpp @@ -0,0 +1,798 @@ +#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 "TGraph.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "TFile.h" + +int main(int argc, char* argv[]) { + jungfrauStyle(); + + if (argc != 5) { + cout << "Correct usage:" << endl; + cout << "arg 1: specify module number" << endl; + cout << "arg 2: specify data location" << endl; + cout << "arg 3: specify pede file" << endl; + cout << "arg 4: specify data file" << endl; + cout << "eg: ./JFMC_DirectBeamScan 006 /data_pool/Module_006_161116 100us_500Hz_QS_allgain_pede_1055 100us_500Hz_QS_beam_scan_a" << endl; + cout << " " << endl; + exit(1); + } + + string module_str = argv[1]; + string data_loc = argv[2]; + string pede_file = argv[3]; + string data_file = argv[4]; + + jungfrauFile *thisfile = new jungfrauFile(); + + jungfrauPixelMask *pixelMaskObject = new jungfrauPixelMask(); + static int pixel_mask [NCH]; + + jungfrauPedestal *pedestalObject = new jungfrauPedestal(); + pedestalObject->pedestalSetNFrames(1000); // using 1000 frames, rolling window + + TH2F* pedestalsG0 = new TH2F("pedestalsG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedestalsG1 = new TH2F("pedestalsG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedestalsG2 = new TH2F("pedestalsG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + TH2F* pedeRMSG0 = new TH2F("pedeRMSG0","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedeRMSG1 = new TH2F("pedeRMSG1","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* pedeRMSG2 = new TH2F("pedeRMSG2","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + char savename[128]; + + int createHistoFile = 1; + + // open pede file + /* KEEP shows where 1M data is + if (runAnalysis == 2) { + thisfile->open((char*)"/data_pool/test_module028_direct_beam_scan__01092016/200us_scan_a_pedestal_QS_%6.6d.dat", 0); + } else if (runAnalysis == 3) { + thisfile->open((char*)"/data_pool/1M_021_022_150916/200us_200Hz_QS_allgain_pede_scan_d_A_%6.6d.dat", 0); + } else if (runAnalysis == 4) { + thisfile->open((char*)"/data_pool/1M_021_022_150916/200us_200Hz_QS_allgain_pede_scan_d_B_%6.6d.dat", 0); + } else if (runAnalysis == 5) { + thisfile->open((char*)"/data_pool/Module_006_161116/100us_500Hz_QS_allgain_pede_1055_%6.6d.dat", 0); + } else if (runAnalysis == 6) { + thisfile->open((char*)"/data_pool/Module_008_181116/150us_500Hz_QS_allgain_pede_1332_%6.6d.dat", 0); + } + */ + sprintf(savename,"%s/%s_%%6.6d.dat", data_loc.c_str(), pede_file.c_str()); + thisfile->open((char*)savename, 0); + + // calculate pixel mask + pixelMaskObject->initialisePixelMask(pixel_mask); + if (module_str == "028") { + pixelMaskObject->maskChip(1, pixel_mask); + pixelMaskObject->maskChip(2, pixel_mask); + pixelMaskObject->maskChip(5, pixel_mask); + pixelMaskObject->maskChip(6, pixel_mask); + } else if (module_str == "006") { + pixelMaskObject->maskChip(6, pixel_mask); + pixelMaskObject->maskSupercolumn(4,1, pixel_mask); + pixelMaskObject->maskSupercolumn(4,2, pixel_mask); + } + + int nevents = 0; + while (thisfile->readNextFrame()) { + nevents++; + } + thisfile->rewind(); + cout << "read " << nevents << " events" << endl; + + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(0, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(1, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + for (int i = 0; i < 999; i++) { + thisfile->readNextFrame(); + pixelMaskObject->maskIfGainNot(3, thisfile->getFrameDataHandle(), (int*)(&pixel_mask)); + } + thisfile->rewind(); + + /* KEEP example of where 1M plots are + if (runAnalysis == 2) { + sprintf(savename,"plots/module28/pixelmask_%d.png", runAnalysis); + } else if (runAnalysis == 3 || runAnalysis == 4) { + sprintf(savename,"plots/1M/DirectBeam/pixelmask_%d.png", runAnalysis); + } else if (runAnalysis == 5) { + sprintf(savename,"plots/Module006Calib/DirectBeam/pixelmask_%d.png", runAnalysis); + } else if (runAnalysis == 6) { + sprintf(savename,"plots/Module008Calib/DirectBeam/pixelmask_%d.png", runAnalysis); + }*/ + sprintf(savename,"plots/M%s/DirectBeam/pixelmask.png", module_str.c_str()); + pixelMaskObject->plotPixelMask(pixel_mask,savename); + cout << "after chip mask, n masked pixels is " << pixelMaskObject->getNMasked(pixel_mask) << endl; + + // caluclate pedestals + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG0->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG0->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalClear(); + + for (int i = 0; i < 1000; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG1->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG1->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalClear(); + + for (int i = 0; i < 999; i++) { + thisfile->readNextFrame(); + pedestalObject->addFrameToPedestalCalculation(thisfile->getFrameDataHandle()); + } + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + pedestalsG2->Fill(i%NC,i/NC,pedestalObject->pedestalOfChannel(i)); + pedeRMSG2->Fill(i%NC,i/NC,pedestalObject->rmsOfChannel(i)); + } + } + pedestalObject->pedestalClear(); + + thisfile->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); + + pedestalsG0->GetXaxis()->SetTitle("Column"); + pedestalsG0->GetYaxis()->SetTitle("Row"); + pedestalsG0->GetYaxis()->SetTitleOffset(0.7); + pedestalsG0->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeG0.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedestalsG1->GetXaxis()->SetTitle("Column"); + pedestalsG1->GetYaxis()->SetTitle("Row"); + pedestalsG1->GetYaxis()->SetTitleOffset(0.7); + pedestalsG1->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeG1.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedestalsG2->GetXaxis()->SetTitle("Column"); + pedestalsG2->GetYaxis()->SetTitle("Row"); + pedestalsG2->GetYaxis()->SetTitleOffset(0.7); + pedestalsG2->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeG2.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG0->GetXaxis()->SetTitle("Column"); + pedeRMSG0->GetYaxis()->SetTitle("Row"); + pedeRMSG0->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG0->GetZaxis()->SetRangeUser(0,25); + pedeRMSG0->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG0.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG1->GetXaxis()->SetTitle("Column"); + pedeRMSG1->GetYaxis()->SetTitle("Row"); + pedeRMSG1->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG1->GetZaxis()->SetRangeUser(0,15); + pedeRMSG1->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG1.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + pedeRMSG2->GetXaxis()->SetTitle("Column"); + pedeRMSG2->GetYaxis()->SetTitle("Row"); + pedeRMSG2->GetYaxis()->SetTitleOffset(0.7); + pedeRMSG2->GetZaxis()->SetRangeUser(0,15); + pedeRMSG2->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/pedeRMSG2.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + int npoints = 40; // 2 to 80 in steps of 2 + + TH2I *adc2d_histo_1_g0 = new TH2I("adc2d_histo_1_g0","",100,0,17000,65536,(65536*0-0.5),(65536*1-0.5)); + TH2I *adc2d_histo_2_g0 = new TH2I("adc2d_histo_2_g0","",100,0,17000,65536,(65536*1-0.5),(65536*2-0.5)); + TH2I *adc2d_histo_3_g0 = new TH2I("adc2d_histo_3_g0","",100,0,17000,65536,(65536*2-0.5),(65536*3-0.5)); + TH2I *adc2d_histo_4_g0 = new TH2I("adc2d_histo_4_g0","",100,0,17000,65536,(65536*3-0.5),(65536*4-0.5)); + TH2I *adc2d_histo_5_g0 = new TH2I("adc2d_histo_5_g0","",100,0,17000,65536,(65536*4-0.5),(65536*5-0.5)); + TH2I *adc2d_histo_6_g0 = new TH2I("adc2d_histo_6_g0","",100,0,17000,65536,(65536*5-0.5),(65536*6-0.5)); + TH2I *adc2d_histo_7_g0 = new TH2I("adc2d_histo_7_g0","",100,0,17000,65536,(65536*6-0.5),(65536*7-0.5)); + TH2I *adc2d_histo_8_g0 = new TH2I("adc2d_histo_8_g0","",100,0,17000,65536,(65536*7-0.5),(65536*8-0.5)); + + TH2I *adc2d_histo_1_g1 = new TH2I("adc2d_histo_1_g1","",100,0,17000,65536,(65536*0-0.5),(65536*1-0.5)); + TH2I *adc2d_histo_2_g1 = new TH2I("adc2d_histo_2_g1","",100,0,17000,65536,(65536*1-0.5),(65536*2-0.5)); + TH2I *adc2d_histo_3_g1 = new TH2I("adc2d_histo_3_g1","",100,0,17000,65536,(65536*2-0.5),(65536*3-0.5)); + TH2I *adc2d_histo_4_g1 = new TH2I("adc2d_histo_4_g1","",100,0,17000,65536,(65536*3-0.5),(65536*4-0.5)); + TH2I *adc2d_histo_5_g1 = new TH2I("adc2d_histo_5_g1","",100,0,17000,65536,(65536*4-0.5),(65536*5-0.5)); + TH2I *adc2d_histo_6_g1 = new TH2I("adc2d_histo_6_g1","",100,0,17000,65536,(65536*5-0.5),(65536*6-0.5)); + TH2I *adc2d_histo_7_g1 = new TH2I("adc2d_histo_7_g1","",100,0,17000,65536,(65536*6-0.5),(65536*7-0.5)); + TH2I *adc2d_histo_8_g1 = new TH2I("adc2d_histo_8_g1","",100,0,17000,65536,(65536*7-0.5),(65536*8-0.5)); + + + // arrays to remember mean and error + // for every step and every pixel + static double adc_g0_means[NCH][40]; + static double adc_g1_means[NCH][40]; + static double adc_g0_meaners[NCH][40]; + static double adc_g1_meaners[NCH][40]; + + // 2D hist for every point, to save + // mean, error and number of fills + // can be loaded from a file + TH2F* avg_adc_g0_map [npoints]; + TH2F* avg_adc_g1_map [npoints]; + TH2F* avg_adcer_g0_map [npoints]; + TH2F* avg_adcer_g1_map [npoints]; + + // creating the histograms + if (createHistoFile) { + + for (int j = 0; j < npoints; j++) { + sprintf(savename,"avg_adc_g0_map_%d", j); + avg_adc_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + sprintf(savename,"avg_adc_g1_map_%d", j); + avg_adc_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + sprintf(savename,"avg_adcer_g0_map_%d", j); + avg_adcer_g0_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + sprintf(savename,"avg_adcer_g1_map_%d", j); + avg_adcer_g1_map[j] = new TH2F(savename,"",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + } + + + int nfiles = 4; + int nframes = 0; + for (int filei = 0; filei < nfiles; filei++) { + + /* KEEP links to data + if (runAnalysis == 2) { + thisfile->open((char*)"/data_pool/test_module028_direct_beam_scan__01092016/200us_scan_a_beam_QS_%6.6d.dat",filei); + } else if (runAnalysis == 3) { + thisfile->open((char*)"/data_pool/1M_021_022_150916/200us_200Hz_QS_beam_scan_d_A_%6.6d.dat",filei); + } else if (runAnalysis == 4) { + thisfile->open((char*)"/data_pool/1M_021_022_150916/200us_200Hz_QS_beam_scan_d_B_%6.6d.dat",filei); + } else if (runAnalysis == 5) { + thisfile->open((char*)"/data_pool/Module_006_161116/100us_500Hz_QS_beam_scan_a_%6.6d.dat",filei); + } else if (runAnalysis == 6) { + thisfile->open((char*)"/data_pool/Module_008_181116/150us_500Hz_QS_beam_scan_b_%6.6d.dat",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()) { + + if (module_str == "028") { + if (nframes%1000 < 100) { // shutter skip + nframes++; + continue; + } + } + + uint16_t* imagedptr = thisfile->getFrameDataHandle(); + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + + uint16_t gain = (imagedptr[i]&0xc000) >> 14; + int adc = imagedptr[i]&0x3fff; + + if (gain == 0) { + if (i < (65536*1)) { + adc2d_histo_1_g0->Fill(adc,i); + } else if (i < (65536*2)) { + adc2d_histo_2_g0->Fill(adc,i); + } else if (i < (65536*3)) { + adc2d_histo_3_g0->Fill(adc,i); + } else if (i < (65536*4)) { + adc2d_histo_4_g0->Fill(adc,i); + } else if (i < (65536*5)) { + adc2d_histo_5_g0->Fill(adc,i); + } else if (i < (65536*6)) { + adc2d_histo_6_g0->Fill(adc,i); + } else if (i < (65536*7)) { + adc2d_histo_7_g0->Fill(adc,i); + } else if (i < (65536*8)) { + adc2d_histo_8_g0->Fill(adc,i); + } + + } else if (gain == 1) { + if (i < (65536*1)) { + adc2d_histo_1_g1->Fill(adc,i); + } else if (i < (65536*2)) { + adc2d_histo_2_g1->Fill(adc,i); + } else if (i < (65536*3)) { + adc2d_histo_3_g1->Fill(adc,i); + } else if (i < (65536*4)) { + adc2d_histo_4_g1->Fill(adc,i); + } else if (i < (65536*5)) { + adc2d_histo_5_g1->Fill(adc,i); + } else if (i < (65536*6)) { + adc2d_histo_6_g1->Fill(adc,i); + } else if (i < (65536*7)) { + adc2d_histo_7_g1->Fill(adc,i); + } else if (i < (65536*8)) { + adc2d_histo_8_g1->Fill(adc,i); + } + } + } + } + nframes++; + + if (nframes%1000 == 0) { + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + + TH1D* projg0; + TH1D* projg1; + if (i < (65536*1)) { + projg0 = adc2d_histo_1_g0->ProjectionX("bin1",i-(65536*(1-1))+1,i-(65536*(1-1))+1); + projg1 = adc2d_histo_1_g1->ProjectionX("bin2",i-(65536*(1-1))+1,i-(65536*(1-1))+1); + } else if (i < (65536*2)) { + projg0 = adc2d_histo_2_g0->ProjectionX("bin3",i-(65536*(2-1))+1,i-(65536*(2-1))+1); + projg1 = adc2d_histo_2_g1->ProjectionX("bin4",i-(65536*(2-1))+1,i-(65536*(2-1))+1); + } else if (i < (65536*3)) { + projg0 = adc2d_histo_3_g0->ProjectionX("bin5",i-(65536*(3-1))+1,i-(65536*(3-1))+1); + projg1 = adc2d_histo_3_g1->ProjectionX("bin6",i-(65536*(3-1))+1,i-(65536*(3-1))+1); + } else if (i < (65536*4)) { + projg0 = adc2d_histo_4_g0->ProjectionX("bin7",i-(65536*(4-1))+1,i-(65536*(4-1))+1); + projg1 = adc2d_histo_4_g1->ProjectionX("bin8",i-(65536*(4-1))+1,i-(65536*(4-1))+1); + } else if (i < (65536*5)) { + projg0 = adc2d_histo_5_g0->ProjectionX("bin9",i-(65536*(5-1))+1,i-(65536*(5-1))+1); + projg1 = adc2d_histo_5_g1->ProjectionX("bin10",i-(65536*(5-1))+1,i-(65536*(5-1))+1); + } else if (i < (65536*6)) { + projg0 = adc2d_histo_6_g0->ProjectionX("bin11",i-(65536*(6-1))+1,i-(65536*(6-1))+1); + projg1 = adc2d_histo_6_g1->ProjectionX("bin12",i-(65536*(6-1))+1,i-(65536*(6-1))+1); + } else if (i < (65536*7)) { + projg0 = adc2d_histo_7_g0->ProjectionX("bin13",i-(65536*(7-1))+1,i-(65536*(7-1))+1); + projg1 = adc2d_histo_7_g1->ProjectionX("bin14",i-(65536*(7-1))+1,i-(65536*(7-1))+1); + } else if (i < (65536*8)) { + projg0 = adc2d_histo_8_g0->ProjectionX("bin15",i-(65536*(8-1))+1,i-(65536*(8-1))+1); + projg1 = adc2d_histo_8_g1->ProjectionX("bin16",i-(65536*(8-1))+1,i-(65536*(8-1))+1); + } + + projg0->Draw(); + adc_g0_means[i][nframes/1000 -1] = projg0->GetMean(); + adc_g0_meaners[i][nframes/1000 -1] = projg0->GetMeanError(); + projg1->Draw(); + adc_g1_means[i][nframes/1000 -1] = projg1->GetMean(); + adc_g1_meaners[i][nframes/1000 -1] = projg1->GetMeanError(); + + projg0->Delete(); // necessary + projg1->Delete(); // necessary + + } + } + adc2d_histo_1_g0->Reset(); + adc2d_histo_1_g1->Reset(); + adc2d_histo_2_g0->Reset(); + adc2d_histo_2_g1->Reset(); + adc2d_histo_3_g0->Reset(); + adc2d_histo_3_g1->Reset(); + adc2d_histo_4_g0->Reset(); + adc2d_histo_4_g1->Reset(); + adc2d_histo_5_g0->Reset(); + adc2d_histo_5_g1->Reset(); + adc2d_histo_6_g0->Reset(); + adc2d_histo_6_g1->Reset(); + adc2d_histo_7_g0->Reset(); + adc2d_histo_7_g1->Reset(); + adc2d_histo_8_g0->Reset(); + adc2d_histo_8_g1->Reset(); + } + + } + + thisfile->close(); + } // end of files + + for (int j = 0; j < npoints; j++) { + int current = (j+1)*2; + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + // fill 2D hists from arrays + + avg_adc_g0_map[j]->SetBinContent((i%NC)+1,(i/NC)+1,adc_g0_means[i][j]); + avg_adc_g1_map[j]->SetBinContent((i%NC)+1,(i/NC)+1,adc_g1_means[i][j]); + + avg_adcer_g0_map[j]->SetBinContent((i%NC)+1,(i/NC)+1,adc_g0_meaners[i][j]); + avg_adcer_g1_map[j]->SetBinContent((i%NC)+1,(i/NC)+1,adc_g1_meaners[i][j]); + + } + } + + avg_adc_g0_map[j]->GetXaxis()->SetTitle("Column"); + avg_adc_g0_map[j]->GetYaxis()->SetTitle("Row"); + avg_adc_g0_map[j]->GetYaxis()->SetTitleOffset(0.7); + avg_adc_g0_map[j]->Draw("colz"); + avg_adc_g0_map[j]->SetMinimum(0); + sprintf(savename,"plots/M%s/DirectBeam/avg_adc_g0_map_%d.png", module_str.c_str(), current); + mapcanvas->SaveAs((const char *)(savename)); + + avg_adc_g1_map[j]->GetXaxis()->SetTitle("Column"); + avg_adc_g1_map[j]->GetYaxis()->SetTitle("Row"); + avg_adc_g1_map[j]->GetYaxis()->SetTitleOffset(0.7); + avg_adc_g1_map[j]->Draw("colz"); + avg_adc_g1_map[j]->SetMinimum(0); + sprintf(savename,"plots/M%s/DirectBeam/avg_adc_g1_map_%d.png", module_str.c_str(), current); + mapcanvas->SaveAs((const char *)(savename)); + + avg_adcer_g0_map[j]->GetXaxis()->SetTitle("Column"); + avg_adcer_g0_map[j]->GetYaxis()->SetTitle("Row"); + avg_adcer_g0_map[j]->GetYaxis()->SetTitleOffset(0.7); + avg_adcer_g0_map[j]->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/avg_adcer_g0_map_%d.png", module_str.c_str(), current); + mapcanvas->SaveAs((const char *)(savename)); + + avg_adcer_g1_map[j]->GetXaxis()->SetTitle("Column"); + avg_adcer_g1_map[j]->GetYaxis()->SetTitle("Row"); + avg_adcer_g1_map[j]->GetYaxis()->SetTitleOffset(0.7); + avg_adcer_g1_map[j]->Draw("colz"); + sprintf(savename,"plots/M%s/DirectBeam/avg_adcer_g1_map_%d.png", module_str.c_str(), current); + mapcanvas->SaveAs((const char *)(savename)); + + } + + // save histograms + /* KEEP links to data + TFile* saved_file; + if (runAnalysis == 2) { + saved_file = new TFile("module28.root","RECREATE"); moved to data/ + } else if (runAnalysis == 3) { + saved_file = new TFile("plots/1M/DirectBeam/1M_DB_A.root","RECREATE"); + } else if (runAnalysis == 4) { + saved_file = new TFile("plots/1M/DirectBeam/1M_DB_B.root","RECREATE"); + } else if (runAnalysis == 5) { + saved_file = new TFile("plots/Module006Calib/DirectBeam/Module006_DB.root","RECREATE"); moved to data/ + } else if (runAnalysis == 6) { + saved_file = new TFile("plots/Module008Calib/DirectBeam/Module008_DB.root","RECREATE"); moved to data/ + } + */ + sprintf(savename,"data/M%s/DB_histos_M%s.root", module_str.c_str(), module_str.c_str()); + TFile* saved_file = new TFile((const char *)(savename),"RECREATE"); + + for (int i = 0; i < npoints; i++) { + avg_adc_g0_map[i]->Write(); + avg_adc_g1_map[i]->Write(); + avg_adcer_g0_map[i]->Write(); + avg_adcer_g1_map[i]->Write(); + } + saved_file->Close(); + } else { + // load histos + cout << "LOADING HISTOS" << endl; + + /* KEEP links to data + TFile* saved_file; + if (runAnalysis == 2) { + saved_file = new TFile("module28.root","READ"); + } else if (runAnalysis == 3) { + saved_file = new TFile("plots/1M/DirectBeam/1M_DB_A.root","READ"); + } else if (runAnalysis == 4) { + saved_file = new TFile("plots/1M/DirectBeam/1M_DB_B.root","READ"); + } else if (runAnalysis == 5) { + saved_file = new TFile("plots/Module006Calib/DirectBeam/Module006_DB.root","READ"); + } else if (runAnalysis == 6) { + saved_file = new TFile("plots/Module008Calib/DirectBeam/Module008_DB.root","READ"); + } + */ + sprintf(savename,"data/M%s/DB_histos_M%s.root", module_str.c_str(), module_str.c_str()); + TFile* saved_file = new TFile((const char *)(savename),"READ"); + + for (int i = 0; i < npoints; i++) { + sprintf(savename,"avg_adc_g0_map_%d", i); + avg_adc_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename)); + sprintf(savename,"avg_adc_g1_map_%d", i); + avg_adc_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename)); + + sprintf(savename,"avg_adcer_g0_map_%d", i); + avg_adcer_g0_map[i] = (TH2F*)saved_file->Get((const char *)(savename)); + sprintf(savename,"avg_adcer_g1_map_%d", i); + avg_adcer_g1_map[i] = (TH2F*)saved_file->Get((const char *)(savename)); + } + } + + TH1F* g0overg1hist = new TH1F("g0overg1hist","",100,-34,-24); + TH2F* g0map = new TH2F("g0map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* g1map = new TH2F("g1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* g0overg1map = new TH2F("g0overg1map","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* g0ermap = new TH2F("g0ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH2F* g1ermap = new TH2F("g1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + TH1F* g0overg1erhist = new TH1F("g0overg1erhist","",100,0,2); + TH2F* g0overg1ermap = new TH2F("g0overg1ermap","",NC,-0.5,NC-0.5,NR,-0.5,NR-0.5); + + for (int i = 0; i < NCH; i++) { + if (pixel_mask[i] == 1) { + + vector r0_adc; + vector r0_adcer; + vector r0_filter; + vector r0_filterer; + vector r1_adc; + vector r1_adcer; + vector r1_filter; + vector r1_filterer; + vector r2_adc; + vector r2_adcer; + vector r2_filter; + vector r2_filterer; + + r0_filter.push_back(0); + r0_filterer.push_back(0); + r0_adc.push_back(pedestalsG0->GetBinContent((i%NC)+1,(i/NC)+1)); + r0_adcer.push_back(pedeRMSG0->GetBinContent((i%NC)+1,(i/NC)+1)); + r1_filter.push_back(0); + r1_filterer.push_back(0); + r1_adc.push_back(pedestalsG1->GetBinContent((i%NC)+1,(i/NC)+1)); + r1_adcer.push_back(pedeRMSG1->GetBinContent((i%NC)+1,(i/NC)+1)); + r2_filter.push_back(0); + r2_filterer.push_back(0); + r2_adc.push_back(pedestalsG2->GetBinContent((i%NC)+1,(i/NC)+1)); + r2_adcer.push_back(pedeRMSG2->GetBinContent((i%NC)+1,(i/NC)+1)); + + for (int j = 0; j < npoints; j++) { + int current = (j+1)*2; + + double this_g0avg = avg_adc_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1); + double this_g1avg = avg_adc_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1); + double this_g0avger = avg_adcer_g0_map[j]->GetBinContent((i%NC)+1,(i/NC)+1); + double this_g1avger = avg_adcer_g1_map[j]->GetBinContent((i%NC)+1,(i/NC)+1); + + if (this_g0avg != 0) { + r0_filter.push_back(current); + r0_filterer.push_back(0); + r0_adc.push_back(this_g0avg); + r0_adcer.push_back(this_g0avger); + } + if (this_g1avg != 0) { + r1_filter.push_back(current); + r1_filterer.push_back(0); + r1_adc.push_back(this_g1avg); + r1_adcer.push_back(this_g1avger); + } + + } + + TGraphErrors *grap_g0 = new TGraphErrors(r0_adc.size(),&(r0_filter[0]),&(r0_adc[0]),&(r0_filterer[0]),&(r0_adcer[0])); + TGraphErrors *grap_g1 = new TGraphErrors(r1_adc.size(),&(r1_filter[0]),&(r1_adc[0]),&(r1_filterer[0]),&(r1_adcer[0])); + TGraphErrors *grap_g2 = new TGraphErrors(r2_adc.size(),&(r2_filter[0]),&(r2_adc[0]),&(r2_filterer[0]),&(r2_adcer[0])); + + grap_g0->SetMarkerStyle(20); + grap_g1->SetMarkerStyle(20); + grap_g2->SetMarkerStyle(20); + grap_g0->SetMarkerColor(kBlue); + grap_g0->SetLineColor(kBlue); + grap_g1->SetMarkerColor(kGreen+2); + grap_g1->SetLineColor(kGreen+2); + grap_g2->SetMarkerColor(kRed); + grap_g2->SetLineColor(kRed); + + TF1 *fit1; + TF1 *fit2; + TF1 *fit1_e; + TF1 *fit2_e; + + if (r0_adc.size() > 1) { + //if (runAnalysis == 2) { + if (module_str == "028") { + fit1 = new TF1("fit1","[0]+[1]*x",2,10); + //} else if (runAnalysis == 3 || runAnalysis == 4) { + } else if (module_str == "021" || module_str == "022") { + fit1 = new TF1("fit1","[0]+[1]*x",2,8); + //fit1 = new TF1("fit1","[0]+[1]*x",2,8); + //} else if (runAnalysis == 5) { + } else if (module_str == "006") { + fit1 = new TF1("fit1","[0]+[1]*x",4,14); + //} else if (runAnalysis == 6) { + } else if (module_str == "008") { + fit1 = new TF1("fit1","[0]+[1]*x",2,20); + } + fit1->SetParameter(0, 100.); + fit1->SetParameter(1, 1); + fit1->SetLineColor(kBlue); + grap_g0->Fit(fit1,"QR",""); + + g0map->Fill(i%NC,i/NC,fit1->GetParameter(1)); + g0ermap->Fill(i%NC,i/NC,fit1->GetParError(1)); + + fit1_e = new TF1("fit1_e","[0]+[1]*x",0,40); + fit1_e->SetParameters(fit1->GetParameter(0),fit1->GetParameter(1)); + fit1_e->SetLineColor(kBlue); + fit1_e->SetLineStyle(2); + } + if (r1_adc.size() > 1) { + //if (runAnalysis == 2 || runAnalysis == 3 || runAnalysis == 4) { + if (module_str == "028" || module_str == "021" || module_str == "022") { + fit2 = new TF1("fit2","[0]+[1]*x",40,80); + //} else if (runAnalysis == 5) { + } else if (module_str == "006") { + fit2 = new TF1("fit2","[0]+[1]*x",60,80); + //} else if (runAnalysis == 6) { + } else if (module_str == "008") { + fit2 = new TF1("fit2","[0]+[1]*x",60,80); + } + fit2->SetParameter(0, 1000.); + fit2->SetParameter(1, -0.1); + fit2->SetLineColor(kGreen+2); + grap_g1->Fit(fit2,"QR",""); + + g1map->Fill(i%NC,i/NC,fit2->GetParameter(1)); + g1ermap->Fill(i%NC,i/NC,fit2->GetParError(1)); + + fit2_e = new TF1("fit2_e","[0]+[1]*x",0,80); + fit2_e->SetParameters(fit2->GetParameter(0),fit2->GetParameter(1)); + fit2_e->SetLineColor(kGreen+2); + fit2_e->SetLineStyle(2); + } + if (r0_adc.size() > 1 && r1_adc.size() > 1) { + g0overg1map->Fill(i%NC,i/NC,fit1->GetParameter(1)/fit2->GetParameter(1)); + g0overg1hist->Fill(fit1->GetParameter(1)/fit2->GetParameter(1)); + g0overg1ermap->Fill(i%NC,i/NC,abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2))); + g0overg1erhist->Fill(abs(fit1->GetParameter(1)/fit2->GetParameter(1))*sqrt(pow((fit1->GetParError(1)/fit1->GetParameter(1)),2) + pow((fit2->GetParError(1)/fit2->GetParameter(1)),2))); + } + + grap_g0->GetXaxis()->SetTitle("Current [mA]"); + grap_g0->GetYaxis()->SetTitle("Average ADC [ADU]"); + grap_g0->GetXaxis()->SetLimits(0,50); + grap_g0->Draw("AP"); + if (r0_adc.size() > 1) { + fit1->Draw("same"); + fit1_e->Draw("same"); + } + if ((i > 51900 && i < 51910) || (i > 52150 && i < 52160)) { + sprintf(savename,"plots/M%s/DirectBeam/grap0_%d.png", module_str.c_str(), i); + mapcanvas->SaveAs((const char *)(savename)); + } + + grap_g1->GetXaxis()->SetTitle("Current [mA]"); + grap_g1->GetYaxis()->SetTitle("Average ADC [ADU]"); + grap_g1->GetXaxis()->SetLimits(0,85); + grap_g1->Draw("AP"); + if (r1_adc.size() > 1) { + fit2->Draw("same"); + fit2_e->Draw("same"); + } + if ((i > 51900 && i < 51910) || (i > 52150 && i < 52160)) { + sprintf(savename,"plots/M%s/DirectBeam/grap1_%d.png", module_str.c_str(), i); + mapcanvas->SaveAs((const char *)(savename)); + } + + grap_g1->SetMinimum(0); + grap_g1->SetMaximum(17000); + grap_g1->Draw("AP"); + if (r1_adc.size() > 1) { + fit2->Draw("same"); + fit2_e->Draw("same"); + } + if (r0_adc.size() > 0) { + grap_g0->Draw("P"); + if (r0_adc.size() > 1) { + fit1->Draw("same"); + fit1_e->Draw("same"); + } + } + if (r2_adc.size() > 0) { + grap_g2->Draw("P"); + } + + if ((i > 51900 && i < 51910) || (i > 52150 && i < 52160) || (i > 123880 && i < 123890)) { + sprintf(savename,"plots/M%s/DirectBeam/graps_%d.png", module_str.c_str(), i); + mapcanvas->SaveAs((const char *)(savename)); + } + + if (r0_adc.size() > 1) { + delete fit1; + delete fit1_e; + } + if (r1_adc.size() > 1) { + delete fit2; + delete fit2_e; + } + + } + } + + g1map->GetXaxis()->SetTitle("Column"); + g1map->GetYaxis()->SetTitle("Row"); + g1map->GetYaxis()->SetTitleOffset(0.7); + g1map->Draw("colz"); + g1map->GetZaxis()->SetRangeUser(-40,0); + sprintf(savename,"plots/M%s/DirectBeam/g1map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0map->GetXaxis()->SetTitle("Column"); + g0map->GetYaxis()->SetTitle("Row"); + g0map->GetYaxis()->SetTitleOffset(0.7); + g0map->Draw("colz"); + g0map->GetZaxis()->SetRangeUser(200,800); + sprintf(savename,"plots/M%s/DirectBeam/g0map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g1ermap->GetXaxis()->SetTitle("Column"); + g1ermap->GetYaxis()->SetTitle("Row"); + g1ermap->GetYaxis()->SetTitleOffset(0.7); + g1ermap->Draw("colz"); + g1ermap->GetZaxis()->SetRangeUser(0,0.7); + sprintf(savename,"plots/M%s/DirectBeam/g1ermap.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0ermap->GetXaxis()->SetTitle("Column"); + g0ermap->GetYaxis()->SetTitle("Row"); + g0ermap->GetYaxis()->SetTitleOffset(0.7); + g0ermap->Draw("colz"); + g0ermap->GetZaxis()->SetRangeUser(2,8); + sprintf(savename,"plots/M%s/DirectBeam/g0ermap.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0overg1map->GetXaxis()->SetTitle("Column"); + g0overg1map->GetYaxis()->SetTitle("Row"); + g0overg1map->GetYaxis()->SetTitleOffset(0.7); + g0overg1map->Draw("colz"); + g0overg1map->GetZaxis()->SetRangeUser(-35,-20); + sprintf(savename,"plots/M%s/DirectBeam/g0overg1map.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0overg1ermap->GetXaxis()->SetTitle("Column"); + g0overg1ermap->GetYaxis()->SetTitle("Row"); + g0overg1ermap->GetYaxis()->SetTitleOffset(0.7); + g0overg1ermap->Draw("colz"); + g0overg1ermap->GetZaxis()->SetRangeUser(0,1); + sprintf(savename,"plots/M%s/DirectBeam/g0overg1ermap.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0overg1hist->GetXaxis()->SetTitle("G0 / G1"); + g0overg1hist->Draw(); + sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + g0overg1erhist->GetXaxis()->SetTitle("G0 / G1 uncert"); + g0overg1erhist->Draw(); + sprintf(savename,"plots/M%s/DirectBeam/g0overg1erhist.png", module_str.c_str()); + mapcanvas->SaveAs((const char *)(savename)); + + /* KEEP links to data + if (runAnalysis == 3 || runAnalysis == 4) { + sprintf(savename,"plots/1M/DirectBeam/DB_ratio_%d.root",runAnalysis); + TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE"); + g0overg1map->Write(); + saved_file2->Close(); + } else if (runAnalysis == 5) { + sprintf(savename,"plots/Module006Calib/DirectBeam/DB_ratio_%d.root",runAnalysis); + TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE"); + g0overg1map->Write(); + g0overg1ermap->Write(); + saved_file2->Close(); + } else if (runAnalysis == 6) { + sprintf(savename,"plots/Module008Calib/DirectBeam/DB_ratio_%d.root",runAnalysis); + TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE"); + g0overg1map->Write(); + g0overg1ermap->Write(); + saved_file2->Close(); + } + */ + sprintf(savename,"data/M%s/DB_ratio_M%s.root", module_str.c_str(), module_str.c_str()); + TFile* saved_file2 = new TFile((const char *)(savename),"RECREATE"); + g0overg1map->Write(); + g0overg1ermap->Write(); + saved_file2->Close(); + +} diff --git a/makefile b/makefile new file mode 100644 index 0000000..4b2a44c --- /dev/null +++ b/makefile @@ -0,0 +1,15 @@ + +JFMC_CuFluoPeak: JFMC_CuFluoPeak.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 JFMC_CuFluoPeak.cpp -o JFMC_CuFluoPeak + +JFMC_CuFluoPeakFit: JFMC_CuFluoPeakFit.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 JFMC_CuFluoPeakFit.cpp -o JFMC_CuFluoPeakFit + +JFMC_DirectBeamScan: JFMC_DirectBeamScan.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 JFMC_DirectBeamScan.cpp -o JFMC_DirectBeamScan + +JFMC_CurrentSourceScan: JFMC_CurrentSourceScan.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 JFMC_CurrentSourceScan.cpp -o JFMC_CurrentSourceScan + +JFMC_CalibWriter: JFMC_CalibWriter.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 JFMC_CalibWriter.cpp -o JFMC_CalibWriter