From: Dennis Doering Date: Tue, 10 Sep 2013 15:51:47 +0000 (+0200) Subject: PlotGraph edit for MIMOSA34 with bugs X-Git-Url: https://jspc29.x-matter.uni-frankfurt.de/git/?a=commitdiff_plain;h=2eb128b44bfb59d320aec5535879e1b9743a2b23;p=radhard.git PlotGraph edit for MIMOSA34 with bugs --- diff --git a/newCOMBI/PlotGraph.C b/newCOMBI/PlotGraph.C new file mode 100755 index 0000000..a6ca3ea --- /dev/null +++ b/newCOMBI/PlotGraph.C @@ -0,0 +1,428 @@ +#include +#include +#include +#include +//#include +//#include +//#include +#include + +#include "CSVRow.h" +//#ifdef __CINT__ +#include "CSVRow.C" // bad C-style... but root wants it this way... no correct separation of declaration and implementation possible. +//#endif + +#include "TFile.h" +#include "TTree.h" +#include "TNtuple.h" +#include "TBranch.h" +#include "TString.h" +#include "TCanvas.h" +#include "TH1F.h" +#include "TLegend.h" +#include "TF1.h" +#include "TMath.h" + + + +#define NUMPIXELS 25 +#define RIGHT_BOUNDARY 400 + +#define VETO_THRESHOLD 20 + + +#define DATAPATH "/d/salt/data_recovered/maps/DennisDoering/root/Analysedata/" +#define USE_SEPARATE_NOISE_RUN false +#define DRAW_SINGLE_HISTOGRAMS true +#define DRAW_SUMMED_HISTOGRAMS false +#define DRAW_VETO_HISTOGRAMS false +#define DRAW_PIXEL_DISTRIBUTION_HISTOGRAMS false +// examples for fit functions: gaus, landau +#define FIT_FUNC "landau" +#define DRAW_FITS true +#define DRAW_NOISE_HISTOGRAM true +#define PRINT_HEADER true + +void PlotGraph() +{ + +Int_t runNoStart=341247; +Int_t anzahlRuns=2; + Int_t runNo; + + + + std::vector histNtuple; + std::vector histNtupleSum; + 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"); + leg->SetFillStyle(0); + +for(Int_t runLauf=0;runLaufIsZombie()) { + // if we cannot open the file, print an error messageForm("hist%i",nHistNtuple) and return immediatly + printf("Error: cannot open %s\n", path.Data()); + return; + } + + // get a pointer to the tree + TNtuple* hitNtuple = (TNtuple*) f->Get("hit"); + TNtuple* hitsNtuple = (TNtuple*) f->Get("hits"); + + // one array element and one branch for each pixel + Float_t pixel[NUMPIXELS]; + TBranch* pixelBranch[NUMPIXELS]; + for (Int_t cnt=0; cntSetBranchAddress(Form("p%i",cnt+1), &pixel[cnt], &pixelBranch[cnt]); + } + UInt_t frame; + TBranch* frameBranch; + hitNtuple->SetBranchAddress("frame", &frame, &frameBranch); + UInt_t seedPixel; + TBranch* seedPixelBranch; + hitNtuple->SetBranchAddress("pixel", &seedPixel, &seedPixelBranch); + + //Float_t frameCounts; + //TBranch* frameCountsBranch; + //hitsNtuple->SetBranchAddress("counts", &frameCounts, &frameCountsBranch); + + // create histogram + histNtuple.push_back(new TH1F(Form("hist%i",nHistNtuple), "Histogram title", 200, 0, RIGHT_BOUNDARY)); + histNtupleSum.push_back(new TH1F(Form("histSum%i",nHistNtuple), "Histogram title", 200, 0, RIGHT_BOUNDARY)); + histNtupleVeto.push_back(new TH1F(Form("histVeto%i",nHistNtuple), "Histogram title", 200, 0, RIGHT_BOUNDARY)); + histPixelDistribution1.push_back(new TH1F(Form("histPixelDistribution1_%i",nHistNtuple), "Histogram title", 50000, 0, 50000)); + histPixelDistribution2.push_back(new TH1F(Form("histPixelDistribution2_%i",nHistNtuple), "Histogram title", 50000, 0, 50000)); + + // loop over all hits to fill histogram + Int_t nentries = hitNtuple->GetEntries(); +// Double_t meanNotSeedSum=0; +// Int_t nhits=0; + + for (Int_t cnt=0; cntGetEntry(cnt); + + // take only frames with a maximum of one hit +// hitsNtuple->GetEntry((Int_t) frame); +// if (frameCounts > 1) +// continue; + + + // histogram with the single pixel + histNtuple[nHistNtuple]->Fill(pixel[12]); + + // histogram with the summed pixels + Double_t pixelSum = 0; + for (Int_t i=0; iFill(pixelSum); + + // search for systematics in pixel distribution + if (pixelSum > 100 && pixelSum < 150) { + histPixelDistribution1[nHistNtuple]->Fill(frame); + } + else if (pixelSum > 180 && pixelSum < 300) { + histPixelDistribution2[nHistNtuple]->Fill(frame); + } + + + + // histogram with the single pixel and with "veto trigger" + // for the "histNtupleVeto" take only hits where only the seed pixel contains charge + Double_t notSeedSum = 0; + for (Int_t i=0; i<12; i++) + notSeedSum += pixel[i]; + for (Int_t i=13; i VETO_THRESHOLD) + continue; +// meanNotSeedSum+=notSeedSum; +// nhits++; + histNtupleVeto[nHistNtuple]->Fill(pixel[12]); // histogram with the single pixel + + } +// cout<<"Mittelwert"<SetLineColor(nHistNtuple+1); + histNtuple[nHistNtuple]->SetStats(kFALSE); + histNtupleSum[nHistNtuple]->SetLineColor(nHistNtuple+1); + histNtupleSum[nHistNtuple]->SetStats(kFALSE); + histNtupleSum[nHistNtuple]->SetLineStyle(2); // make summed graphs dashed + histNtupleVeto[nHistNtuple]->SetLineColor(nHistNtuple+1); + histNtupleVeto[nHistNtuple]->SetStats(kFALSE); + histNtupleVeto[nHistNtuple]->SetLineStyle(3); // make summed graphs dotted + histPixelDistribution1[nHistNtuple]->SetLineColor(1); // black line + histPixelDistribution2[nHistNtuple]->SetLineColor(2); // red line +// TString legendEntry = Form("#splitline{Run %i: T_{set}=%.1f, T_{sens}=%.1f, Chip=%i, }{Source=%s, Matrix=%i, Submatrix=%i, RadDose=%.2E}", runNo, temperature, tempSens, chip, TString(source).Data(), matrix6480, submatrix, radDose); + TString legendEntry = TString(Form("Run %i: ",runNo )); +// TString legendEntry = TString(Form("test") + (selTemp==-1 ? Form(", T_{set}=%.1f",temperature):"");// + Form(", T_{sens}=%.1f",tempSens) + (selChip==-1 ? Form(", Chip=%i",chip):"") + (selSource=="-1" ? Form(", Source=%s",source):"") + (selMatrix6480==-1 ? Form(", Matrix=%i",matrix6480):"") + (selSubmatrix==-1 ? Form(", Submatrix=%i",submatrix):"") + (selRadDose==-1 ? Form(", RadDose=%.2E",radDose):""); + leg->AddEntry(histNtuple[nHistNtuple], legendEntry); +// leg->AddEntry((TObject*)0, "", ""); + +// histNtuple[nHistNtuple]->SetLineColor(nHistNtuple); +// histNtuple[nHistNtuple]->Draw(); + +// cout << "Hits per Frame: " << nentries/frame << endl; + + + //------------------------------------------------- + // calculate median of noise (with error estimation) + + + // open the file + TFile* f_N; + TString path_N; + + path_N = TString(DATAPATH) + Form("/%i_0.root", runNo); + 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 + + + + + + + // determine maximum y-value + Float_t maxY = 0; + if (DRAW_SINGLE_HISTOGRAMS){ + for (std::vector::iterator it = histNtuple.begin(); it != histNtuple.end(); ++it){ + Int_t thismax = (*it)->GetMaximum(); + if (thismax > maxY){ + maxY = thismax; + // maxHistNum = std::distance(histNtuple.begin(), it); + } + } + } + if (DRAW_SUMMED_HISTOGRAMS){ + for (std::vector::iterator it = histNtupleSum.begin(); it != histNtupleSum.end(); ++it){ + Int_t thismax = (*it)->GetMaximum(); + if (thismax > maxY){ + maxY = thismax; + // maxHistNum = std::distance(histNtuple.begin(), it); + } + } + } + if (DRAW_VETO_HISTOGRAMS){ + for (std::vector::iterator it = histNtupleVeto.begin(); it != histNtupleVeto.end(); ++it){ + Int_t thismax = (*it)->GetMaximum(); + if (thismax > maxY){ + maxY = thismax; + // maxHistNum = std::distance(histNtuple.begin(), it); + } + } + } + + + TCanvas* c1 = new TCanvas(); + c1->SetTitle("Testtitel"); + 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){ + Int_t cnt = std::distance(histNtuple.begin(), it) ; + // (*it)->SetLineColor(cnt); + (*it)->SetMaximum(1.05*maxY); + (*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()); + 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); + + legendEntry = TString(Form("Run %i: %.1f %.1f %.1f",runNo,posMax,minFitMax,maxFitMax )); + + leg->AddEntry((TObject*) 0, legendEntry, ""); + leg->SetTextSize(0.02); + leg->Draw(); + + } + (*it)->GetXaxis()->UnZoom(); + // cout << "Maximum bin: " << (*it)->GetBinCenter((*it)->GetMaximumBin()) << endl; + // cout << "Maximum of gaussian fit: " << fitFunc->GetParameter(1) << "\n\n"; + fitMaxPosX.push_back(posMax); + // runDetails[cnt] = runDetails[cnt] + Form("%.2f", fitMaxPosX); + } + } + } + // plot all summed histograms into the same histogram + if (DRAW_SUMMED_HISTOGRAMS){ + for (std::vector::iterator it = histNtupleSum.begin(); it != histNtupleSum.end(); ++it){ + Int_t cnt = std::distance(histNtupleSum.begin(), it) ; + (*it)->SetMaximum(1.05*maxY); + if (!DRAW_SINGLE_HISTOGRAMS) + (*it)->Draw((cnt==0)?"":"same"); + else + (*it)->Draw("same"); + + if (DRAW_FITS) { + // (*it)->GetXaxis()->SetRange((*it)->GetXaxis()->FindBin(fitMaxPosX[cnt] + 25),(*it)->GetXaxis()->FindBin(RIGHT_BOUNDARY)); // look only for maxima with x greater than 25 more right than the peak of the seed pixel + (*it)->GetXaxis()->SetRange((*it)->GetXaxis()->FindBin(fitMaxPosX[cnt]),(*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-15, xValMax+15); + (*it)->GetXaxis()->UnZoom(); + fitFunc->DrawCopy("same"); + // cout << "Maximum bin: " << (*it)->GetBinCenter((*it)->GetMaximumBin()) << endl; + // cout << "Maximum of gaussian fit: " << fitFunc->GetParameter(1) << "\n\n"; + fitMaxPosXSum.push_back(fitFunc->GetParameter(1)); + } + } + } + // plot single pixel histograms with veto trigger in the same histogram + if (DRAW_VETO_HISTOGRAMS){ + for (std::vector::iterator it = histNtupleVeto.begin(); it != histNtupleVeto.end(); ++it){ + Int_t cnt = std::distance(histNtupleVeto.begin(), it) ; + (*it)->SetMaximum(1.05*maxY); + if (!DRAW_SINGLE_HISTOGRAMS && !DRAW_SUMMED_HISTOGRAMS) + (*it)->Draw((cnt==0)?"":"same"); + else + (*it)->Draw("same"); + + if (DRAW_FITS) { + // (*it)->GetXaxis()->SetRange((*it)->GetXaxis()->FindBin(fitMaxPosX[cnt] + 25),(*it)->GetXaxis()->FindBin(RIGHT_BOUNDARY)); // look only for maxima with x greater than 25 more right than the peak of the seed pixel + (*it)->GetXaxis()->SetRange((*it)->GetXaxis()->FindBin(fitMaxPosX[cnt]),(*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-15, xValMax+15); + (*it)->GetXaxis()->UnZoom(); + fitFunc->DrawCopy("same"); + // cout << "Maximum bin: " << (*it)->GetBinCenter((*it)->GetMaximumBin()) << endl; + // cout << "Maximum of gaussian fit: " << fitFunc->GetParameter(1) << "\n\n"; + fitMaxPosXVeto.push_back(fitFunc->GetParameter(1)); + } + } + } + // add legend + legendEntry = TString(Form("Run %i: ",runNo )); + + leg->AddEntry((TObject*) 0, legendEntry, ""); + leg->SetTextSize(0.02); + leg->Draw(); + + + // plot histogram for analysis of pixel distribution in another histogram + if (DRAW_PIXEL_DISTRIBUTION_HISTOGRAMS){ + TCanvas* c2 = new TCanvas(); + for (std::vector::iterator it = histPixelDistribution1.begin(); it != histPixelDistribution1.end(); ++it){ + Int_t cnt = std::distance(histPixelDistribution1.begin(), it) ; + (*it)->Draw((cnt==0)?"":"same"); + } + for (std::vector::iterator it = histPixelDistribution2.begin(); it != histPixelDistribution2.end(); ++it){ + (*it)->Draw("same"); + } + } + + } //for-Schleife runnumber + // print all information + TString runDetailsHeader = TString("runNo\t") + (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