]> jspc29.x-matter.uni-frankfurt.de Git - radhard.git/commitdiff
Added landau fits for Sr90-source and noise calculation with error
authorStefan Strohauer <sstrohauer@jspc31.x-matter.uni-frankfurt.de>
Thu, 29 Aug 2013 11:03:03 +0000 (13:03 +0200)
committerStefan Strohauer <sstrohauer@jspc31.x-matter.uni-frankfurt.de>
Thu, 29 Aug 2013 11:03:03 +0000 (13:03 +0200)
estimation.

PlotGraph/PlotGraph.C
ProcessMeasurements/analyzeRun.C

index 39bcc283216f48c5ada3bac2d7199c3561398656..1f1c75ad2c2f73b89199636832a0a5e0d2bce448 100644 (file)
@@ -21,6 +21,7 @@
 #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)
@@ -73,11 +78,17 @@ void PlotGraph(Float_t selTemp, Int_t selChip, std::string selSource, Int_t selM
   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");
@@ -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<NUMPIXELS; i++)
              notSeedSum += pixel[i];
-           if (abs(notSeedSum) > 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; 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 
        }
@@ -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<TH1F*>::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; 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;
   }
   
   
index acfafd0dea2a689d9610b59846de6a257178231f..98e3bea25689944a92510011a1ae5c6a1b2ecc3b 100644 (file)
@@ -31,6 +31,7 @@
 #include "MAPS.h"
 
 #define NUMEVENTSFORNOISE 1000
+#define OFFSET 1
 
 void analyzeRun(const char *const inputDataPath, const char *const outputDataPath, Int_t runNo, std::string matrixAB, Bool_t useSeperateNoiseRun, Bool_t runIsNoise, Int_t noiseRunNo, Bool_t overwrite, std::string logfileName)
 {
@@ -79,7 +80,7 @@ void analyzeRun(const char *const inputDataPath, const char *const outputDataPat
     
     if (runIsNoise) {
       MAPS *mapsi = new MAPS(inputDataPath,outputDataPath,runNo,matrixNoRows,matrixNoCoulmns,matrixAB.c_str(),submatrix);
-      mapsi->getNoise(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;