Index: /fact/tools/rootmacros/PulseTemplates/templateextractors.C =================================================================== --- /fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 14741) +++ /fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 14742) @@ -2,8 +2,14 @@ #include #include +#include + +#include +#include +#include #include "templateextractors.h" using namespace std; +using namespace TMath; @@ -144,4 +150,5 @@ Double_t *x = new Double_t[nbins]; Double_t *y = new Double_t[nbins]; + for (Int_t i=0;iGetXaxis()->GetBinCenter(position->at(i) ); @@ -156,4 +163,84 @@ //---------------------------------------------------------------------------- +int ExtractTH1EnriesToVector ( + TH1* inputHisto, //histogram from which median will be calculated + vector* 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;iGetBinContent(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,4 +339,45 @@ //---------------------------------------------------------------------------- +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 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,34 +388,45 @@ ) { + 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 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 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,5 +434,5 @@ return; } -// end of PlotMedianEachSliceOfPulse +// end of ErrorMedianOfTH2Slices //---------------------------------------------------------------------------- @@ -380,4 +519,7 @@ float median = 0; float mean = 0; + float max_prop_err = 0; + float median_err = 0; + float mean_err = 0; if (verbosityLevel > 3) @@ -430,8 +572,18 @@ // 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,5 +591,17 @@ 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,19 +610,34 @@ } 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