Changeset 14707

Timestamp:

11/30/12 16:24:47 (12 years ago)

Author:

Jens Buss

Message:

implemented error calculation for median

Location:

fact/tools/rootmacros/PulseTemplates

Files:

• templateextractors.C (modified) (8 diffs)
• templateextractors.h (modified) (3 diffs)

fact/tools/rootmacros/PulseTemplates/templateextractors.C

@@ -86 +86 @@
         if (verbosityLevel > 2) cout << "\t...calculation of "
                                      << "pulse order " << pulse_order;
-        //  vector max_value_of to number of timeslices in Overlay Spectra
+        //  vector max_value_of to number of slices in Overlay Spectra
         CalcMaxPropabilityOfSlice(
                     hInputHistoArray[pulse_order],
@@ -135 +135 @@
 //----------------------------------------------------------------------------
 
+Double_t MedianOfH1withRndSlices (
+        TH1*            inputHisto,     //histogram from which median will be calculated
+        vector<int>*    position        //array with random slice numbers
+        )
+{
+   //compute the median for 1-d histogram h1
+   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
+   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) );
+      y[i] = inputHisto->GetBinContent(position->at(i) );
+   }
+   Double_t median = TMath::Median(nbins,x,y);
+   delete [] x;
+   delete [] y;
+   return median;
+}
+// end of PlotMedianEachSliceOfPulse
+//----------------------------------------------------------------------------
+
 void
 PlotMedianOfSlice(
@@ -148 +169 @@
     TH2F**  hInputHistoArray    = NULL;
     TH1F**  hOutputHistoArray   = NULL;
-    TH1*    hTempHisto          = NULL;
-    float   median              = 0;
+//    TH1*    hTempHisto          = NULL;
+//    float   median              = 0;
 
     if (overlayType.Contains("Edge"))
@@ -175 +196 @@
                                      << "pulse order " << pulse_order;
 
-        Int_t   nbins       = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
-
-        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
-
-            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
-            median      = MedianOfH1(hTempHisto);
-
-            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
-
-            hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, median );
+        MedianOfTH2Slices(
+                hInputHistoArray[pulse_order],
+                hOutputHistoArray[pulse_order],
+                verbosityLevel
+                );
+
+        ErrorMedianOfTH2Slices(
+                hInputHistoArray[pulse_order],
+                hOutputHistoArray[pulse_order],
+                5,  //numIterations
+                verbosityLevel
+                );
+
+//        Int_t   nbins       = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
+
+//        for (Int_t slice=1;slice<=nbins;slice++) {
+
+//            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
+//            median      = MedianOfH1(hTempHisto);
+
+//            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
+
+//            hOutputHistoArray[pulse_order]->SetBinContent(slice, median );
+////            delete h1;
+//        }
+
+        if (verbosityLevel > 2) cout << "\t...done" << endl;
+    }
+}
+// end of PlotMedianEachSliceOfPulse
+//----------------------------------------------------------------------------
+
+void
+MedianOfTH2Slices(
+        TH2*        inputHisto,
+        TH1*        outputHisto,
+        int         verbosityLevel
+        )
+{
+    Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
+    float   median              = 0;
+
+    for (Int_t slice=1;slice<=nbins;slice++) {
+
+        median      = MedianOfH1( inputHisto->ProjectionY("",slice,slice) );
+
+        if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
+
+        outputHisto->SetBinContent(slice, median );
 //            delete h1;
+    }
+    return;
+}
+// end of MedianOfTH2Slices
+//----------------------------------------------------------------------------
+
+void
+ErrorMedianOfTH2Slices(
+        TH2*        inputHisto,
+        TH1*        outputHisto,
+        int         numIterations,
+        int         verbosityLevel
+        )
+{
+    Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
+
+//    float MedianOfSliceMean[nbins];
+//    float MedianOfSliceRMS[nbins];
+
+    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;
+
+
+        Sample sample(sample_size);
+        for (int i = 0; i < numIterations; i++)
+        {
+            sample.BootstrapSample(rndList, sample_min, sample_max, sample_size);
+            median[i] = MedianOfH1withRndSlices(
+                        inputHisto->ProjectionY("",slice,slice),
+                        rndList
+                        );
+
+            if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
         }
-
-        if (verbosityLevel > 2) cout << "\t...done" << endl;
-    }
+//        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
+//        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
+//        outputHisto->SetBinError(slice, MedianOfSliceRMS);
+        outputHisto->SetBinError(slice, RMS(numIterations, median) );
+
+//        outputHisto->SetBinContent(slice, median );
+//            delete h1;
+    }
+    return;
 }
 // end of PlotMedianEachSliceOfPulse
@@ -236 +340 @@
         Int_t   nbins       = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
 
-        for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
-
-            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
+        for (Int_t slice=1;slice<=nbins;slice++) {
+
+            hTempHisto  = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
             mean      = hTempHisto->GetMean();
 
-            if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",TimeSlice,mean);
-
-            hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, mean );
+            if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",slice,mean);
+
+            hOutputHistoArray[pulse_order]->SetBinContent(slice, mean );
 //            delete h1;
         }
@@ -326 +430 @@
 //    if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
 
-    for (Int_t TimeSlice=1;TimeSlice<=300;TimeSlice++)
-    {
-
-        hTempHisto  = hInputHisto->ProjectionY("",TimeSlice,TimeSlice);
+    for (Int_t slice=1;slice<=300;slice++)
+    {
+
+        hTempHisto  = hInputHisto->ProjectionY("",slice,slice);
 
         if (verbosityLevel > 3)
         {
-            cout << "\t\t...calculating maxProb of slice " << TimeSlice << endl;
+            cout << "\t\t...calculating maxProb of slice " << slice << endl;
         }
         max_prop    = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
@@ -339 +443 @@
         if (verbosityLevel > 3)
         {
-            cout << "\t\t...calculating Median of slice " << TimeSlice << endl;
+            cout << "\t\t...calculating Median of slice " << slice << endl;
         }
         median      = MedianOfH1(hTempHisto);
 
-        if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << TimeSlice << endl;
+        if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << slice << endl;
         mean        = hTempHisto->GetMean();
 
-        if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << TimeSlice << ": " << max_prop << endl;
-        hOutputMaxHisto->SetBinContent(TimeSlice, max_prop );
-
-        if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << TimeSlice << ": " << mean << endl;
-        hOutputMeanHisto->SetBinContent(TimeSlice, mean );
-
-        if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << TimeSlice << ": " << median << endl;
-        hOutputMedianHisto->SetBinContent(TimeSlice, median );
+        if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
+        hOutputMaxHisto->SetBinContent(slice, max_prop );
+
+        if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
+        hOutputMeanHisto->SetBinContent(slice, mean );
+
+        if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
+        hOutputMedianHisto->SetBinContent(slice, median );
         delete hTempHisto;
 
-    }//Loop over Timeslices
+    }//Loop over slices
 }
 // end of PlotMaxPropabilityPulse
@@ -423 +527 @@
     out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
 
-    for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
-    {
-        out << TimeSlice << "," ;
-        out << Max_histo->GetBinContent(TimeSlice) << ",";
-        out << Mean_histo->GetBinContent(TimeSlice) << ",";
-        out << Median_histo->GetBinContent(TimeSlice) << endl;
+    for (int slice=1;slice<=nbins;slice++)
+    {
+        out << slice << "," ;
+        out << Max_histo->GetBinContent(slice) << ",";
+        out << Mean_histo->GetBinContent(slice) << ",";
+        out << Median_histo->GetBinContent(slice) << endl;
     }
     out.close();

fact/tools/rootmacros/PulseTemplates/templateextractors.h

@@ -13 +13 @@
 
 #include "pixel.h"
+#include "Sample.h"
 
 /** Assignment operator.
@@ -57 +58 @@
         );
 
+Double_t MedianOfH1withRndSlices (
+        TH1*            inputHisto,     //histogram from which median will be calculated
+        vector<int> *position        //array with random slice numbers
+        );
+
 void
 PlotMedianOfSlice(
@@ -69 +75 @@
         TString         overlayType,
         int             verbosityLevel
+        );
+
+void
+MedianOfTH2Slices(
+        TH2*        inputHisto,
+        TH1*        outputHisto,
+        int         verbosityLevel
+        );
+
+void
+ErrorMedianOfTH2Slices(
+        TH2*        inputHisto,
+        TH1*        outputHisto,
+        int         numIterations,
+        int         verbosityLevel
         );