Index: fact/tools/rootmacros/PulseTemplates/templateextractors.C =================================================================== --- fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 14706) +++ fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 14707) @@ -86,5 +86,5 @@ 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,4 +135,25 @@ //---------------------------------------------------------------------------- +Double_t MedianOfH1withRndSlices ( + TH1* inputHisto, //histogram from which median will be calculated + vector* 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;iGetXaxis()->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,6 +169,6 @@ TH2F** hInputHistoArray = NULL; TH1F** hOutputHistoArray = NULL; - TH1* hTempHisto = NULL; - float median = 0; +// TH1* hTempHisto = NULL; +// float median = 0; if (overlayType.Contains("Edge")) @@ -175,19 +196,102 @@ << "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 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,12 +340,12 @@ 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,12 +430,12 @@ // 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,22 +443,22 @@ 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,10 +527,10 @@ 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(); Index: fact/tools/rootmacros/PulseTemplates/templateextractors.h =================================================================== --- fact/tools/rootmacros/PulseTemplates/templateextractors.h (revision 14706) +++ fact/tools/rootmacros/PulseTemplates/templateextractors.h (revision 14707) @@ -13,4 +13,5 @@ #include "pixel.h" +#include "Sample.h" /** Assignment operator. @@ -57,4 +58,9 @@ ); +Double_t MedianOfH1withRndSlices ( + TH1* inputHisto, //histogram from which median will be calculated + vector *position //array with random slice numbers + ); + void PlotMedianOfSlice( @@ -69,4 +75,19 @@ TString overlayType, int verbosityLevel + ); + +void +MedianOfTH2Slices( + TH2* inputHisto, + TH1* outputHisto, + int verbosityLevel + ); + +void +ErrorMedianOfTH2Slices( + TH2* inputHisto, + TH1* outputHisto, + int numIterations, + int verbosityLevel );