Changeset 14742 for fact/tools/rootmacros

Timestamp:

12/12/12 18:29:42 (12 years ago)

Author:

Jens Buss

Message:

add ExtractTH1EnriesToVector,  CalculateErrorsWithBootstrapping;
modified ErrorMedianOfH1 to use BootstrapSample; add gaussian fit to
slice distribution to improve max_prob result; add error calculation of
slices' values with bootstrapping

File:

• fact/tools/rootmacros/PulseTemplates/templateextractors.C (modified) (10 diffs)

fact/tools/rootmacros/PulseTemplates/templateextractors.C

@@ -2 +2 @@
 #include <fstream>
 #include <stdlib.h>
+#include <iostream>
+
+#include <TCanvas.h>
+#include <TTimer.h>
+#include <Getline.h>
 
 #include "templateextractors.h"
 
 using namespace std;
+using namespace TMath;
 
 
@@ -144 +150 @@
    Double_t *x = new Double_t[nbins];
    Double_t *y = new Double_t[nbins];
+
    for (Int_t i=0;i<nbins;i++) {
       x[i] = inputHisto->GetXaxis()->GetBinCenter(position->at(i) );
@@ -156 +163 @@
 //----------------------------------------------------------------------------
 
+int ExtractTH1EnriesToVector (
+        TH1*            inputHisto,     //histogram from which median will be calculated
+        vector<int>*    entries        //array with random slice numbers
+        )
+{
+   //compute number of bins in imput histo
+   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
+   int quantity;
+   for (Int_t i=0;i<nbins;i++) {
+       quantity = inputHisto->GetBinContent(i);
+       for (int j = 0; j < quantity; j++){
+           entries->push_back(inputHisto->GetXaxis()->GetBinCenter(i));
+       }
+   }
+   return entries->size();
+}
+// end of PlotMedianEachSliceOfPulse
+//----------------------------------------------------------------------------
+
+void CalculateErrorsWithBootstrapping(TH1* inputHisto, int numIt, double* par, double* parErr)
+{
+    cout << "...calculating errors of " << inputHisto->GetName() << endl;
+    Sample  sample;
+    double  Mean[numIt];
+    double  Median[numIt];
+    double  Max[numIt];
+    double  parameter[3];
+    double  parameterErr[3];
+
+    for (int i = 0; i < numIt; i++)
+    {
+        TH1D    tempHisto(
+                    "tempHisto",
+                    "tempHisto",
+                    inputHisto->GetNbinsX(),
+                    inputHisto->GetXaxis()->GetFirst(),
+                    inputHisto->GetXaxis()->GetLast()
+                    );
+
+        sample.BootstrapTH1(inputHisto, &tempHisto);
+
+        Mean[i]     = tempHisto.GetMean();
+        Median[i]   = MedianOfH1 ( &tempHisto );
+
+        Max[i]      = tempHisto.GetBinCenter( tempHisto.GetMaximumBin() );
+
+        //improved determination of maximum
+//        TF1 gaus("fgaus", "gaus", Max[i]-10, Max[i]+10);
+//        tempHisto.Fit("fgaus", "WWQRN0");
+//        Max[i]      = gaus.GetParameter(1);
+
+
+        sample.SetSeed(sample.GetSeed() + 1);
+
+    }
+
+    parameter[0]    = TMath::Mean(  numIt, Max);
+    parameterErr[0] = TMath::RMS(   numIt, Max);
+    parameter[1]    = TMath::Mean(  numIt, Median);
+    parameterErr[1] = TMath::RMS(   numIt, Median);
+    parameter[2]    = TMath::Mean(  numIt, Mean);
+    parameterErr[2] = TMath::RMS(   numIt, Mean);
+
+    for (int i = 0; i < 3; i++)
+    {
+        if (&par[i] != NULL)
+        {
+            par[i] =  parameter[i];
+        }
+        else cout << "par[" << i << "] array to small" << endl;
+        if (&par[i] != NULL)
+        {
+            parErr[i] =  parameterErr[i];
+        }
+        else cout << "parErr[" << i << "] array to small" << endl;
+    }
+
+    return;
+}
+
 void
 PlotMedianOfSlice(
@@ -252 +339 @@
 //----------------------------------------------------------------------------
 
+double
+ErrorMedianOfH1(
+        TH1*        inputHisto,
+        int         numIterations,
+        int         verbosityLevel
+        )
+{
+    if (verbosityLevel > 5) cout << "...calculating errors of median" << endl;
+//    float MedianOfSliceMean;
+    double medianRMS = 0;
+
+    int sample_min  = inputHisto->GetXaxis()->GetFirst();
+    int sample_max  = inputHisto->GetXaxis()->GetLast();
+    int sample_size = inputHisto->GetXaxis()->GetNbins();
+//    int sample_size = sample_max - sample_min;
+//    cout << "sample_min " << sample_min
+//         << ", sample_max " << sample_max
+//         << ", sample_size " << sample_size
+//         << endl;
+    double   median[numIterations];
+    vector<int> rndList;
+    Sample sample(sample_size);
+
+    for (int i = 0; i < numIterations; i++)
+    {
+        sample.SetSeed(sample.GetSeed() + 1);
+        sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
+        median[i] = MedianOfH1withRndSlices(
+                    inputHisto,
+                    &rndList
+                    );
+    }
+
+//    MedianOfSliceMean    = TMath::Mean(numIterations, median);
+    medianRMS     = RMS(numIterations, median);
+//    if (verbosityLevel > 2) printf("Mean Median=%g\n",MedianOfSliceMean);
+    if (verbosityLevel > 5) printf("medianRMS=%g\n",medianRMS);
+
+    return medianRMS;
+}
+
 void
 ErrorMedianOfTH2Slices(
@@ -260 +388 @@
         )
 {
+    if (verbosityLevel > 2) cout << "...calculating errors of median" << endl;
     Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
-
-//    float MedianOfSliceMean[nbins];
-//    float MedianOfSliceRMS[nbins];
+    cout << "nbins " << nbins << endl;
+    float MedianOfSliceMean[nbins];
+    float MedianOfSliceRMS[nbins];
+    float median[numIterations];
+    vector<int> rndList;
+//    rndList->clear();
+
+    TH1 *htemp = NULL;
 
     for (Int_t slice=1;slice<=nbins;slice++) {
-
-        float   median[numIterations];
-        int sample_min  = inputHisto->GetXaxis()->GetFirst();
-        int sample_max  = inputHisto->GetXaxis()->GetLast();
-        int sample_size = sample_max - sample_min;
-        vector<int> rndList;
-
+        htemp   = inputHisto->ProjectionY("",slice,slice);
+
+
+        int sample_min  = htemp->GetXaxis()->GetFirst();
+        int sample_max  = htemp->GetXaxis()->GetLast();
+        int sample_size = htemp->GetXaxis()->GetNbins();
+
+//        cout << "sample_min " << sample_min
+//             << ", sample_max " << sample_max
+//             << ", sample_size " << sample_size
+//             << endl;
 
         Sample sample(sample_size);
         for (int i = 0; i < numIterations; i++)
         {
-            sample.BootstrapSample(rndList, sample_min, sample_max, sample_size);
+            sample.SetSeed(sample.GetSeed() + 1);
+            sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
             median[i] = MedianOfH1withRndSlices(
-                        inputHisto->ProjectionY("",slice,slice),
-                        rndList
+                        htemp,
+                        &rndList
                         );
 
-            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
         }
-//        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
-//        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
-//        outputHisto->SetBinError(slice, MedianOfSliceRMS);
-        outputHisto->SetBinError(slice, RMS(numIterations, median) );
-
+        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
+        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
+        outputHisto->SetBinError(slice, MedianOfSliceRMS[slice]);
+//        outputHisto->SetBinError(slice, RMS(numIterations, median) );
+        if (verbosityLevel > 2) printf("Mean Median of Slice %d, Mean Median=%g\n",slice,MedianOfSliceMean[slice]);
+        if (verbosityLevel > 2) printf("RMS of Median of Slice %d, MedianRMS=%g\n",slice,MedianOfSliceRMS[slice]);
 //        outputHisto->SetBinContent(slice, median );
 //            delete h1;
@@ -295 +434 @@
     return;
 }
-// end of PlotMedianEachSliceOfPulse
+// end of ErrorMedianOfTH2Slices
 //----------------------------------------------------------------------------
 
@@ -380 +519 @@
     float   median                  = 0;
     float   mean                    = 0;
+    float   max_prop_err            = 0;
+    float   median_err              = 0;
+    float   mean_err                = 0;
 
     if (verbosityLevel > 3)
@@ -430 +572 @@
 //    if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
 
-    for (Int_t slice=1;slice<=300;slice++)
+    int slice_min  = hInputHisto->GetXaxis()->GetFirst();
+    int slice_max  = hInputHisto->GetXaxis()->GetLast();
+
+    for (Int_t slice=slice_min;slice<=slice_max;slice++)
     {
 
         hTempHisto  = hInputHisto->ProjectionY("",slice,slice);
+
+        //containers for errors
+        double slMean[3];
+        double slError[3];
+
+        //calculate Errors with bootstrapping
+        CalculateErrorsWithBootstrapping(hTempHisto, 100, slMean, slError);
 
         if (verbosityLevel > 3)
@@ -439 +591 @@
             cout << "\t\t...calculating maxProb of slice " << slice << endl;
         }
+
+        //Get maximum of slice's distribution
         max_prop    = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
+
+        //improve result by< fitting gaussian to slices distribution
+        TF1 gaus("fgaus", "gaus", max_prop-30, max_prop+30);
+        hTempHisto->Fit("fgaus", "QRN0");
+        max_prop        = gaus.GetParameter(1);
+
+        //calculate error of max prop
+//        max_prop_err    = gaus.GetParameter(2); //RMS of gaus fit of slices distribution
+        max_prop_err    = slError[0];         //error calculated with bootstrapping
+//        max_prop_err    = hTempHisto->GetRMS();   //RMS of slice
 
         if (verbosityLevel > 3)
@@ -446 +610 @@
         }
         median      = MedianOfH1(hTempHisto);
-
-        if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << slice << endl;
+        median_err    = slError[1]; //error from bootstraping
+
+        if (verbosityLevel > 3) cout << "\t\t...calculating Mean of slice " << slice << endl;
+
+
         mean        = hTempHisto->GetMean();
+        mean_err    = hTempHisto->GetRMS(); //RMS of slice
+//        mean_err    = slError[2];         //error from bootstraping
 
         if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
         hOutputMaxHisto->SetBinContent(slice, max_prop );
+        hOutputMaxHisto->SetBinError(    slice,  max_prop_err);
 
         if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
         hOutputMeanHisto->SetBinContent(slice, mean );
+        hOutputMeanHisto->SetBinError(    slice,  mean_err);
 
         if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
-        hOutputMedianHisto->SetBinContent(slice, median );
+        hOutputMedianHisto->SetBinContent(  slice,  median );
+        hOutputMedianHisto->SetBinError(    slice,  median_err);
         delete hTempHisto;
 
     }//Loop over slices
+
+//    ErrorMedianOfTH2Slices(
+//            hInputHisto,
+//            hOutputMedianHisto,
+//            100,  //numIterations
+//            verbosityLevel
+//            );
 }
 // end of PlotMaxPropabilityPulse