#include "TH1F.h"
#include "TLegend.h"
#include "TF1.h"
+#include "TMath.h"
#define COLRUNNO 0
#define COLTEMP 2
#define LABORBUCH "RelevantRuns.csv"
#define DATAPATH "/local/sstrohauer/data/data_8.8.2013_19:21/"
+#define USE_SEPARATE_NOISE_RUN true
#define DRAW_SINGLE_HISTOGRAMS true
-#define DRAW_SUMMED_HISTOGRAMS true
-#define DRAW_VETO_HISTOGRAMS true
-#define DRAW_PIXEL_DISTRIBUTION_HISTOGRAMS false
-#define DRAW_FITS false
+#define DRAW_SUMMED_HISTOGRAMS false
+#define DRAW_VETO_HISTOGRAMS false
+#define DRAW_PIXEL_DISTRIBUTION_HISTOGRAMS true
+// examples for fit functions: gaus, landau
+#define FIT_FUNC "landau"
+#define DRAW_FITS true
+#define DRAW_NOISE_HISTOGRAM true
#define PRINT_HEADER false
void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selMatrix6480, Int_t selSubmatrix, Float_t selRadDose)
std::vector<TH1F*> histNtupleVeto;
std::vector<TH1F*> histPixelDistribution1;
std::vector<TH1F*> histPixelDistribution2;
+ std::vector<TH1F*> histNoise;
Int_t nHistNtuple = 0;
std::vector<std::string> runDetails;
std::vector<Float_t> fitMaxPosX;
std::vector<Float_t> fitMaxPosXSum;
std::vector<Float_t> fitMaxPosXVeto;
+ std::vector<Float_t> fitLandauErrorLeft;
+ std::vector<Float_t> fitLandauErrorRight;
+ std::vector<Float_t> noiseMedian;
+ std::vector<Float_t> noiseLowest17percent;
+ std::vector<Float_t> noiseHighest17percent;
TLegend* leg = new TLegend(0.5,0.5,0.89,0.89);//(0.6,0.7,0.89,0.89);
// leg->SetHeader();//"Legend Title");
notSeedSum += pixel[i];
for (Int_t i=13; i<NUMPIXELS; i++)
notSeedSum += pixel[i];
- if (abs(notSeedSum) > VETO_THRESHOLD)
+ if (TMath::Abs(notSeedSum) > VETO_THRESHOLD)
continue;
// meanNotSeedSum+=notSeedSum;
// nhits++;
// histNtuple[nHistNtuple]->SetLineColor(nHistNtuple);
// histNtuple[nHistNtuple]->Draw();
- cout << "Hits per Frame: " << nentries/frame << endl;
+// cout << "Hits per Frame: " << nentries/frame << endl;
+
+
+ //-------------------------------------------------
+ // calculate median of noise (with error estimation)
+ Int_t runNo_N; // run number
+ Float_t temperature_N; // desired temperature
+ Int_t chip_N; // chip number (eg. 1-6)
+ std::string source_N;
+ Int_t matrix6480_N; // which matrix?
+ Float_t radDose_N;
+ Bool_t runFinished_N; // measurement already done?
+ Int_t noiseRunNo = 0;
+
+ std::ifstream file_N(LABORBUCH);
+ CSVRow row_N;
+ file_N >> row_N; // read first row (header) of the laboratory book
+ while (file_N >> row_N) {
+ try {
+ if (row_N[COLRUNNO] == "") // no relevant data in this row
+ continue;
+ runNo_N = atoi(row_N[COLRUNNO].c_str()); // run number (conversion works like runNo = stoi(std::string ("300");)
+ temperature_N = atof(row_N[COLTEMP].c_str()); // temperature
+ chip_N = atoi(row_N[COLCHIP].c_str());
+ source_N = row_N[COLSOURCE];
+ matrix6480_N = atoi(row_N[COLMATRIX].c_str()); // which matrix?
+ radDose_N = atof(row_N[COLRADDOSE].c_str());
+ runFinished_N = row_N[COLDONE] == "y"; // measurement already done?
+
+ if (!runFinished_N)
+ continue;
+
+ if ( (temperature_N == selTemp || selTemp == -1)
+ && (chip_N == selChip || selChip == -1)
+ && (source_N == "none")
+ && (matrix6480_N == selMatrix6480 || selMatrix6480 == -1)
+ && (radDose_N == selRadDose || selRadDose == -1))
+ {
+ noiseRunNo = runNo_N;
+ }
+ }
+ catch(...){
+ cout << "Error while reading laboratory book (run number " << runNo << ")!\n";
+ }
+ }
+ if (noiseRunNo == 0)
+ cout << "No corresponding noise run found!\n";
+
+ // open the file
+ TFile* f_N;
+ TString path_N;
+ if (USE_SEPARATE_NOISE_RUN) {
+ path_N = TString(DATAPATH) + Form("%i/%i_0_%i.root", noiseRunNo, noiseRunNo, submatrix);
+ f_N = TFile::Open(path_N);
+ if (f_N->IsZombie()) {
+ // if we cannot open the file, print an error messageForm("hist%i",nHistNtuple) and return immediatly
+ printf("Error: cannot open %s\n", path_N.Data());
+ return;
+ }
+ }
+ else {
+ path_N = TString(DATAPATH) + "noise/" + Form("%i/%i_0_%i.root", noiseRunNo, noiseRunNo, submatrix);
+ f_N = TFile::Open(path_N);
+ if (f_N->IsZombie()) {
+ // if we cannot open the file, print an error messageForm("hist%i",nHistNtuple) and return immediatly
+ printf("Error: cannot open %s\n", path_N.Data());
+ return;
+ }
+ }
+
+ // get a pointer to the tree
+ TNtuple* noiseNtuple = (TNtuple*) f_N->Get("noise");
+
+ Float_t noise;
+ TBranch* noiseBranch;
+ noiseNtuple->SetBranchAddress("noise", &noise, &noiseBranch);
+ // create histogram
+ histNoise.push_back(new TH1F(Form("noise%i",nHistNtuple), "Noisy title", 10000, 0, 10));
+// histNoise.push_back(new TH1F(Form("noise%i",nHistNtuple), "Noisy title", 500, 0, 10));
+
+ Int_t nentries_N = noiseNtuple->GetEntries();
+ for (Int_t cnt=0; cnt<nentries_N; cnt++) {
+ noiseNtuple->GetEntry(cnt);
+ histNoise[nHistNtuple]->Fill(noise);
+ }
+
+ // get median and error estimation
+ Double_t const probabilities[] = {0.3415/2, 0.5, 1-0.3415/2}; // sigma/2 from gaus to the left and to the right //{0.17, 0.5, 1-0.17};
+ Double_t quantiles[3];
+ histNoise[nHistNtuple]->GetQuantiles( 3, quantiles, probabilities);
+ noiseLowest17percent.push_back((Float_t) quantiles[0]); // left error bar (lowest 17% noise values are below this quantile)
+ noiseMedian.push_back((Float_t) quantiles[1]); // median of the noise distribution
+ noiseHighest17percent.push_back((Float_t) quantiles[2]); // right error bar (highest 17% noise values are higher than this quantile)
+// cout << "q1: " << quantiles[1] - quantiles[0] << "\tq2: " << quantiles[1] << "\tq3: " << quantiles[2] - quantiles[1] << endl;
+
+ if (DRAW_NOISE_HISTOGRAM) {
+ new TCanvas();
+ histNoise[nHistNtuple]->Draw();
+ }
+ //-------------------------------------------------
+
+
+
+
+
+
+
nHistNtuple++; // increase
}
TCanvas* c1 = new TCanvas();
c1->SetTitle("Testtitel");
- TF1* fitFunc = new TF1("fitFunc","gaus",0,RIGHT_BOUNDARY);
+ TF1* fitFunc = new TF1("fitFunc",FIT_FUNC,0,RIGHT_BOUNDARY);
// plot every single pixel histogram in one canvas
if (DRAW_SINGLE_HISTOGRAMS){
for (std::vector<TH1F*>::iterator it = histNtuple.begin(); it != histNtuple.end(); ++it){
(*it)->Draw((cnt==0)?"":"same");
if (DRAW_FITS) {
+ Float_t posMax = 0;
(*it)->GetXaxis()->SetRange((*it)->GetXaxis()->FindBin(20),(*it)->GetXaxis()->FindBin(RIGHT_BOUNDARY)); // look only for maxima with x greater than 20
Int_t xValMax = (*it)->GetBinCenter((*it)->GetMaximumBin());
- (*it)->Fit(fitFunc, "N,Q,W", "", xValMax-10, xValMax+10);
+ if (TString(FIT_FUNC)=="gaus") {
+ (*it)->Fit(fitFunc, "N,Q,W", "", xValMax-10, xValMax+10);
+ posMax = fitFunc->GetParameter(1);
+ fitFunc->DrawCopy("same");
+ }
+ else if (TString(FIT_FUNC)=="landau") {
+ Float_t fitMax1 = 1000;
+ Float_t fitMax2 = 1000;
+ Float_t fitMax3 = 1000;
+ Float_t minFitMax = 1000;
+ Float_t maxFitMax = 1000;
+
+ (*it)->Fit(fitFunc, "N,Q,W", "", 20, RIGHT_BOUNDARY);
+ fitMax1 = fitFunc->GetParameter(1);
+ fitFunc->DrawCopy("same");
+ (*it)->Fit(fitFunc, "N,Q,W", "", 20, fitMax1);
+ fitMax2 = fitFunc->GetParameter(1);
+ fitFunc->SetLineColor(kBlue);
+ fitFunc->SetLineStyle(2); // dashed
+ fitFunc->DrawCopy("same");
+ (*it)->Fit(fitFunc, "N,Q,W", "", fitMax1, RIGHT_BOUNDARY);
+ fitMax3 = fitFunc->GetParameter(1);
+ fitFunc->SetLineColor(kGreen);
+ fitFunc->DrawCopy("same");
+ fitFunc->SetLineStyle(1); // normal for the following fits
+
+ // Sort the three fits and save error estimation
+ minFitMax = TMath::Min(TMath::Min(fitMax1,fitMax2),fitMax3);
+ maxFitMax = TMath::Max(TMath::Max(fitMax1,fitMax2),fitMax3);
+ posMax = fitMax1 + fitMax2 + fitMax3 - minFitMax - maxFitMax;
+ fitLandauErrorLeft.push_back(posMax - minFitMax);
+ fitLandauErrorRight.push_back(maxFitMax - posMax);
+ }
(*it)->GetXaxis()->UnZoom();
- fitFunc->DrawCopy("same");
// cout << "Maximum bin: " << (*it)->GetBinCenter((*it)->GetMaximumBin()) << endl;
// cout << "Maximum of gaussian fit: " << fitFunc->GetParameter(1) << "\n\n";
- fitMaxPosX.push_back(fitFunc->GetParameter(1));
+ fitMaxPosX.push_back(posMax);
// runDetails[cnt] = runDetails[cnt] + Form("%.2f", fitMaxPosX);
}
}
// print all information
- TString runDetailsHeader = TString("runNo\tT\tTSens\tChip\tSource\tMatrix\tRadDose\tSubmtrx") + (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS ? "\tColPeak":"") + (DRAW_SUMMED_HISTOGRAMS && DRAW_FITS ? "\tSumColPeak":"") + (DRAW_VETO_HISTOGRAMS && DRAW_FITS ? "\tCalibPeak":"") + "\n";
+ TString runDetailsHeader = TString("runNo\tT\tTSens\tChip\tSource\tMatrix\tRadDose\t\tSubmtrx") + (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS ? "\tColPeak":"") + (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS && TString(FIT_FUNC)=="landau" ? "\tErrL":"") + (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS && TString(FIT_FUNC)=="landau" ? "\tErrR":"") + (DRAW_SUMMED_HISTOGRAMS && DRAW_FITS ? "\tSumColPeak":"") + (DRAW_VETO_HISTOGRAMS && DRAW_FITS ? "\tCalPeak":"") + "\tNoise" + "\tNseLo17" + "\tNseHi17" + "\tS/N" + "\n";
if (PRINT_HEADER)
cout << endl << runDetailsHeader;
for (UInt_t i=0; i<runDetails.size(); ++i){
// cout << runDetails[i] << Form("\t%.2f", fitMaxPosX[i]) << Form("\t%.2f", fitMaxPosXSum[i]) << endl;
- cout << runDetails[i] << (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f", fitMaxPosX[i]):"") << (DRAW_SUMMED_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f", fitMaxPosXSum[i]):"") << (DRAW_VETO_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f", fitMaxPosXVeto[i]):"") << endl;
+ cout << runDetails[i] << (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f", fitMaxPosX[i]):"") << (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS && TString(FIT_FUNC)=="landau" ? Form("\t%.2f", fitLandauErrorLeft[i]):"") << (DRAW_SINGLE_HISTOGRAMS && DRAW_FITS && TString(FIT_FUNC)=="landau" ? Form("\t%.2f", fitLandauErrorRight[i]):"") << (DRAW_SUMMED_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f\t", fitMaxPosXSum[i]):"") << (DRAW_VETO_HISTOGRAMS && DRAW_FITS ? Form("\t%.2f", fitMaxPosXVeto[i]):"") << Form("\t%.2f", noiseMedian[i]) << Form("\t%.2f", noiseLowest17percent[i]) << Form("\t%.2f", noiseHighest17percent[i]) << Form("\t%.2f", fitMaxPosX[i]/noiseMedian[i]) /*S/N-ratio*/ << endl;
}
jspc29 / radhard.git / commitdiff