detectors/JFCalibration
0
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
6ca856c5bb620000bca876a51f8ab13163361498
JFCalibration/JFMC_DirectBeamScan.cpp
redford_s 6ca856c5bb Style and removing old code.
2017-01-09 14:34:44 +01:00

717 lines
27 KiB
C++
Raw Blame History

#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 = 0;
// 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 == "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();
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++) {
// 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
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;
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<double> r0_adc;
vector<double> r0_adcer;
vector<double> r0_filter;
vector<double> r0_filterer;
vector<double> r1_adc;
vector<double> r1_adcer;
vector<double> r1_filter;
vector<double> r1_filterer;
vector<double> r2_adc;
vector<double> r2_adcer;
vector<double> r2_filter;
vector<double> 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));
TCanvas* c1 = new TCanvas("c1","");
g0overg1hist->GetXaxis()->SetTitle("G0 / G1");
g0overg1hist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/g0overg1hist.png", module_str.c_str());
c1->SaveAs((const char *)(savename));
g0overg1erhist->GetXaxis()->SetTitle("G0 / G1 uncert");
g0overg1erhist->Draw();
sprintf(savename,"plots/M%s/DirectBeam/g0overg1erhist.png", module_str.c_str());
c1->SaveAs((const char *)(savename));
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();
}
Reference in New Issue View Git Blame Copy Permalink