Changeset 14742


Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed

  • fact/tools/rootmacros/PulseTemplates/templateextractors.C

    r14707 r14742  
    22#include <fstream>
    33#include <stdlib.h>
     4#include <iostream>
     5
     6#include <TCanvas.h>
     7#include <TTimer.h>
     8#include <Getline.h>
    49
    510#include "templateextractors.h"
    611
    712using namespace std;
     13using namespace TMath;
    814
    915
    … …  
    144150   Double_t *x = new Double_t[nbins];
    145151   Double_t *y = new Double_t[nbins];
     152
    146153   for (Int_t i=0;i<nbins;i++) {
    147154      x[i] = inputHisto->GetXaxis()->GetBinCenter(position->at(i) );
    … …  
    156163//----------------------------------------------------------------------------
    157164
     165int ExtractTH1EnriesToVector (
     166        TH1*            inputHisto,     //histogram from which median will be calculated
     167        vector<int>*    entries        //array with random slice numbers
     168        )
     169{
     170   //compute number of bins in imput histo
     171   Int_t nbins = inputHisto->GetXaxis()->GetNbins();
     172   int quantity;
     173   for (Int_t i=0;i<nbins;i++) {
     174       quantity = inputHisto->GetBinContent(i);
     175       for (int j = 0; j < quantity; j++){
     176           entries->push_back(inputHisto->GetXaxis()->GetBinCenter(i));
     177       }
     178   }
     179   return entries->size();
     180}
     181// end of PlotMedianEachSliceOfPulse
     182//----------------------------------------------------------------------------
     183
     184void CalculateErrorsWithBootstrapping(TH1* inputHisto, int numIt, double* par, double* parErr)
     185{
     186    cout << "...calculating errors of " << inputHisto->GetName() << endl;
     187    Sample  sample;
     188    double  Mean[numIt];
     189    double  Median[numIt];
     190    double  Max[numIt];
     191    double  parameter[3];
     192    double  parameterErr[3];
     193
     194    for (int i = 0; i < numIt; i++)
     195    {
     196        TH1D    tempHisto(
     197                    "tempHisto",
     198                    "tempHisto",
     199                    inputHisto->GetNbinsX(),
     200                    inputHisto->GetXaxis()->GetFirst(),
     201                    inputHisto->GetXaxis()->GetLast()
     202                    );
     203
     204        sample.BootstrapTH1(inputHisto, &tempHisto);
     205
     206        Mean[i]     = tempHisto.GetMean();
     207        Median[i]   = MedianOfH1 ( &tempHisto );
     208
     209        Max[i]      = tempHisto.GetBinCenter( tempHisto.GetMaximumBin() );
     210
     211        //improved determination of maximum
     212//        TF1 gaus("fgaus", "gaus", Max[i]-10, Max[i]+10);
     213//        tempHisto.Fit("fgaus", "WWQRN0");
     214//        Max[i]      = gaus.GetParameter(1);
     215
     216
     217        sample.SetSeed(sample.GetSeed() + 1);
     218
     219    }
     220
     221    parameter[0]    = TMath::Mean(  numIt, Max);
     222    parameterErr[0] = TMath::RMS(   numIt, Max);
     223    parameter[1]    = TMath::Mean(  numIt, Median);
     224    parameterErr[1] = TMath::RMS(   numIt, Median);
     225    parameter[2]    = TMath::Mean(  numIt, Mean);
     226    parameterErr[2] = TMath::RMS(   numIt, Mean);
     227
     228    for (int i = 0; i < 3; i++)
     229    {
     230        if (&par[i] != NULL)
     231        {
     232            par[i] =  parameter[i];
     233        }
     234        else cout << "par[" << i << "] array to small" << endl;
     235        if (&par[i] != NULL)
     236        {
     237            parErr[i] =  parameterErr[i];
     238        }
     239        else cout << "parErr[" << i << "] array to small" << endl;
     240    }
     241
     242    return;
     243}
     244
    158245void
    159246PlotMedianOfSlice(
    … …  
    252339//----------------------------------------------------------------------------
    253340
     341double
     342ErrorMedianOfH1(
     343        TH1*        inputHisto,
     344        int         numIterations,
     345        int         verbosityLevel
     346        )
     347{
     348    if (verbosityLevel > 5) cout << "...calculating errors of median" << endl;
     349//    float MedianOfSliceMean;
     350    double medianRMS = 0;
     351
     352    int sample_min  = inputHisto->GetXaxis()->GetFirst();
     353    int sample_max  = inputHisto->GetXaxis()->GetLast();
     354    int sample_size = inputHisto->GetXaxis()->GetNbins();
     355//    int sample_size = sample_max - sample_min;
     356//    cout << "sample_min " << sample_min
     357//         << ", sample_max " << sample_max
     358//         << ", sample_size " << sample_size
     359//         << endl;
     360    double   median[numIterations];
     361    vector<int> rndList;
     362    Sample sample(sample_size);
     363
     364    for (int i = 0; i < numIterations; i++)
     365    {
     366        sample.SetSeed(sample.GetSeed() + 1);
     367        sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
     368        median[i] = MedianOfH1withRndSlices(
     369                    inputHisto,
     370                    &rndList
     371                    );
     372    }
     373
     374//    MedianOfSliceMean    = TMath::Mean(numIterations, median);
     375    medianRMS     = RMS(numIterations, median);
     376//    if (verbosityLevel > 2) printf("Mean Median=%g\n",MedianOfSliceMean);
     377    if (verbosityLevel > 5) printf("medianRMS=%g\n",medianRMS);
     378
     379    return medianRMS;
     380}
     381
    254382void
    255383ErrorMedianOfTH2Slices(
    … …  
    260388        )
    261389{
     390    if (verbosityLevel > 2) cout << "...calculating errors of median" << endl;
    262391    Int_t   nbins       = inputHisto->GetXaxis()->GetNbins();
    263 
    264 //    float MedianOfSliceMean[nbins];
    265 //    float MedianOfSliceRMS[nbins];
     392    cout << "nbins " << nbins << endl;
     393    float MedianOfSliceMean[nbins];
     394    float MedianOfSliceRMS[nbins];
     395    float median[numIterations];
     396    vector<int> rndList;
     397//    rndList->clear();
     398
     399    TH1 *htemp = NULL;
    266400
    267401    for (Int_t slice=1;slice<=nbins;slice++) {
    268 
    269         float   median[numIterations];
    270         int sample_min  = inputHisto->GetXaxis()->GetFirst();
    271         int sample_max  = inputHisto->GetXaxis()->GetLast();
    272         int sample_size = sample_max - sample_min;
    273         vector<int> rndList;
    274 
     402        htemp   = inputHisto->ProjectionY("",slice,slice);
     403
     404
     405        int sample_min  = htemp->GetXaxis()->GetFirst();
     406        int sample_max  = htemp->GetXaxis()->GetLast();
     407        int sample_size = htemp->GetXaxis()->GetNbins();
     408
     409//        cout << "sample_min " << sample_min
     410//             << ", sample_max " << sample_max
     411//             << ", sample_size " << sample_size
     412//             << endl;
    275413
    276414        Sample sample(sample_size);
    277415        for (int i = 0; i < numIterations; i++)
    278416        {
    279             sample.BootstrapSample(rndList, sample_min, sample_max, sample_size);
     417            sample.SetSeed(sample.GetSeed() + 1);
     418            sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
    280419            median[i] = MedianOfH1withRndSlices(
    281                         inputHisto->ProjectionY("",slice,slice),
    282                         rndList
     420                        htemp,
     421                        &rndList
    283422                        );
    284423
    285             if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
    286424        }
    287 //        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
    288 //        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
    289 //        outputHisto->SetBinError(slice, MedianOfSliceRMS);
    290         outputHisto->SetBinError(slice, RMS(numIterations, median) );
    291 
     425        MedianOfSliceMean[slice]    = TMath::Mean(numIterations, median);
     426        MedianOfSliceRMS[slice]     = TMath::RMS(numIterations, median);
     427        outputHisto->SetBinError(slice, MedianOfSliceRMS[slice]);
     428//        outputHisto->SetBinError(slice, RMS(numIterations, median) );
     429        if (verbosityLevel > 2) printf("Mean Median of Slice %d, Mean Median=%g\n",slice,MedianOfSliceMean[slice]);
     430        if (verbosityLevel > 2) printf("RMS of Median of Slice %d, MedianRMS=%g\n",slice,MedianOfSliceRMS[slice]);
    292431//        outputHisto->SetBinContent(slice, median );
    293432//            delete h1;
    … …  
    295434    return;
    296435}
    297 // end of PlotMedianEachSliceOfPulse
     436// end of ErrorMedianOfTH2Slices
    298437//----------------------------------------------------------------------------
    299438
    … …  
    380519    float   median                  = 0;
    381520    float   mean                    = 0;
     521    float   max_prop_err            = 0;
     522    float   median_err              = 0;
     523    float   mean_err                = 0;
    382524
    383525    if (verbosityLevel > 3)
    … …  
    430572//    if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
    431573
    432     for (Int_t slice=1;slice<=300;slice++)
     574    int slice_min  = hInputHisto->GetXaxis()->GetFirst();
     575    int slice_max  = hInputHisto->GetXaxis()->GetLast();
     576
     577    for (Int_t slice=slice_min;slice<=slice_max;slice++)
    433578    {
    434579
    435580        hTempHisto  = hInputHisto->ProjectionY("",slice,slice);
     581
     582        //containers for errors
     583        double slMean[3];
     584        double slError[3];
     585
     586        //calculate Errors with bootstrapping
     587        CalculateErrorsWithBootstrapping(hTempHisto, 100, slMean, slError);
    436588
    437589        if (verbosityLevel > 3)
    … …  
    439591            cout << "\t\t...calculating maxProb of slice " << slice << endl;
    440592        }
     593
     594        //Get maximum of slice's distribution
    441595        max_prop    = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
     596
     597        //improve result by< fitting gaussian to slices distribution
     598        TF1 gaus("fgaus", "gaus", max_prop-30, max_prop+30);
     599        hTempHisto->Fit("fgaus", "QRN0");
     600        max_prop        = gaus.GetParameter(1);
     601
     602        //calculate error of max prop
     603//        max_prop_err    = gaus.GetParameter(2); //RMS of gaus fit of slices distribution
     604        max_prop_err    = slError[0];         //error calculated with bootstrapping
     605//        max_prop_err    = hTempHisto->GetRMS();   //RMS of slice
    442606
    443607        if (verbosityLevel > 3)
    … …  
    446610        }
    447611        median      = MedianOfH1(hTempHisto);
    448 
    449         if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << slice << endl;
     612        median_err    = slError[1]; //error from bootstraping
     613
     614        if (verbosityLevel > 3) cout << "\t\t...calculating Mean of slice " << slice << endl;
     615
     616
    450617        mean        = hTempHisto->GetMean();
     618        mean_err    = hTempHisto->GetRMS(); //RMS of slice
     619//        mean_err    = slError[2];         //error from bootstraping
    451620
    452621        if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
    453622        hOutputMaxHisto->SetBinContent(slice, max_prop );
     623        hOutputMaxHisto->SetBinError(    slice,  max_prop_err);
    454624
    455625        if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
    456626        hOutputMeanHisto->SetBinContent(slice, mean );
     627        hOutputMeanHisto->SetBinError(    slice,  mean_err);
    457628
    458629        if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
    459         hOutputMedianHisto->SetBinContent(slice, median );
     630        hOutputMedianHisto->SetBinContent(  slice,  median );
     631        hOutputMedianHisto->SetBinError(    slice,  median_err);
    460632        delete hTempHisto;
    461633
    462634    }//Loop over slices
     635
     636//    ErrorMedianOfTH2Slices(
     637//            hInputHisto,
     638//            hOutputMedianHisto,
     639//            100,  //numIterations
     640//            verbosityLevel
     641//            );
    463642}
    464643// end of PlotMaxPropabilityPulse
Note: See TracChangeset for help on using the changeset viewer.