From 6675cfcf42b449e76a1a7da31c417100c57f3b79 Mon Sep 17 00:00:00 2001 From: Stefan Strohauer Date: Thu, 29 Aug 2013 13:03:03 +0200 Subject: [PATCH] Added landau fits for Sr90-source and noise calculation with error estimation. --- PlotGraph/PlotGraph.C | 173 ++++++++++++++++++++++++++++--- ProcessMeasurements/analyzeRun.C | 5 +- 2 files changed, 164 insertions(+), 14 deletions(-) diff --git a/PlotGraph/PlotGraph.C b/PlotGraph/PlotGraph.C index 39bcc28..1f1c75a 100644 --- a/PlotGraph/PlotGraph.C +++ b/PlotGraph/PlotGraph.C @@ -21,6 +21,7 @@ #include "TH1F.h" #include "TLegend.h" #include "TF1.h" +#include "TMath.h" #define COLRUNNO 0 #define COLTEMP 2 @@ -38,11 +39,15 @@ #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) @@ -73,11 +78,17 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM std::vector histNtupleVeto; std::vector histPixelDistribution1; std::vector histPixelDistribution2; + std::vector histNoise; Int_t nHistNtuple = 0; std::vector runDetails; std::vector fitMaxPosX; std::vector fitMaxPosXSum; std::vector fitMaxPosXVeto; + std::vector fitLandauErrorLeft; + std::vector fitLandauErrorRight; + std::vector noiseMedian; + std::vector noiseLowest17percent; + std::vector 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"); @@ -203,7 +214,7 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM notSeedSum += pixel[i]; for (Int_t i=13; i VETO_THRESHOLD) + if (TMath::Abs(notSeedSum) > VETO_THRESHOLD) continue; // meanNotSeedSum+=notSeedSum; // nhits++; @@ -232,7 +243,113 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM // 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; cntGetEntry(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 } @@ -279,7 +396,7 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM 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::iterator it = histNtuple.begin(); it != histNtuple.end(); ++it){ @@ -289,14 +406,46 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM (*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); } } @@ -368,12 +517,12 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM // 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; igetNoise(NUMEVENTSFORNOISE,0,0); + mapsi->getNoise(NUMEVENTSFORNOISE,OFFSET,0); delete mapsi; logfile << runNo << "\t" << matrixAB << "\t" << submatrix << "\tprocessed as NOISE\n"; } @@ -87,7 +88,7 @@ void analyzeRun(const char *const inputDataPath, const char *const outputDataPat std::string noiseOutputDataPath = std::string(outputDataPath); if (!useSeperateNoiseRun){ MAPS *mapsiNoise = new MAPS(inputDataPath,std::string(outputDataPath)+"noise/",runNo,matrixNoRows,matrixNoCoulmns,matrixAB.c_str(),submatrix); - mapsiNoise->getNoise(NUMEVENTSFORNOISE,0,0); + mapsiNoise->getNoise(NUMEVENTSFORNOISE,OFFSET,0); delete mapsiNoise; logfile << runNo << "\t" << matrixAB << "\t" << submatrix << "\tprocessed as NOISE\n"; noiseRunNo = runNo; -- 2.43.0