--- /dev/null
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <string>
+//#include <sys/types.h>
+//#include <sys/stat.h>
+//#include <sys/time.h>
+#include <stdlib.h>
+
+#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<TH1F*> histNtuple;
+ std::vector<TH1F*> histNtupleSum;
+ 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");
+ leg->SetFillStyle(0);
+
+for(Int_t runLauf=0;runLauf<anzahlRuns;runLauf++)
+{
+Int_t runNo=runNoStart+runLauf;
+
+ cout<<runNo<<endl;
+
+ // open the file
+ TString path = TString(DATAPATH) + Form("%i_0.root", runNo);
+ TFile* f = TFile::Open(path);
+ if (f->IsZombie()) {
+ // 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; cnt<NUMPIXELS; cnt++) {
+ hitNtuple->SetBranchAddress(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; cnt<nentries; cnt++) {
+ hitNtuple->GetEntry(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; i<NUMPIXELS; i++){
+ pixelSum += pixel[i];
+ }
+ histNtupleSum[nHistNtuple]->Fill(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<NUMPIXELS; i++)
+ notSeedSum += pixel[i];
+ if (TMath::Abs(notSeedSum) > VETO_THRESHOLD)
+ continue;
+// meanNotSeedSum+=notSeedSum;
+// nhits++;
+ histNtupleVeto[nHistNtuple]->Fill(pixel[12]); // histogram with the single pixel
+
+ }
+// cout<<"Mittelwert"<<meanNotSeedSum/nhits<<endl;
+
+
+ histNtuple[nHistNtuple]->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; 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
+
+
+
+
+
+
+ // determine maximum y-value
+ Float_t maxY = 0;
+ if (DRAW_SINGLE_HISTOGRAMS){
+ for (std::vector<TH1F*>::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<TH1F*>::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<TH1F*>::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<TH1F*>::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<TH1F*>::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<TH1F*>::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<TH1F*>::iterator it = histPixelDistribution1.begin(); it != histPixelDistribution1.end(); ++it){
+ Int_t cnt = std::distance(histPixelDistribution1.begin(), it) ;
+ (*it)->Draw((cnt==0)?"":"same");
+ }
+ for (std::vector<TH1F*>::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<anzahlRuns; ++i){
+// cout << runDetails[i] << Form("\t%.2f", fitMaxPosX[i]) << Form("\t%.2f", fitMaxPosXSum[i]) << endl;
+ cout << runNoStart+i<<"\t"<<(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]) << endl;
+
+}
+
+}
+
+
jspc29 / radhard.git / commitdiff