Index: /fact/tools/rootmacros/PulseTemplates/templateextractors.C =================================================================== --- /fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 13652) +++ /fact/tools/rootmacros/PulseTemplates/templateextractors.C (revision 13652) @@ -0,0 +1,430 @@ +#include +#include + +#include "templateextractors.h" +#include +using namespace std; + + +void +CalcMaxPropabilityOfSlice( + TH2* inputHisto, + TH1* outputHisto, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating slieces value of max. propability" + << endl; + + float max_value_of_slice = 0; + TH1D* hTempHisto = NULL; + int nbins = inputHisto->GetXaxis()->GetNbins(); + + if (verbosityLevel > 2) cout << "...generating projections" << endl; + + for (Int_t TimeSlice = 1; TimeSlice <= nbins; TimeSlice++) + { + if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice; + + //put maximumvalue of every Y-projection of every timeslice into vector + hTempHisto = inputHisto->ProjectionY( "", TimeSlice,TimeSlice); + max_value_of_slice = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() ); + + if (verbosityLevel > 2) cout << " with max value " + << max_value_of_slice << endl; + + outputHisto->SetBinContent(TimeSlice, max_value_of_slice ); + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; +} +// end of CalcMaxPropabilityOfSlice +//---------------------------------------------------------------------------- + + +void +PlotMaxPropabilityPulse( + Pixel* CurrentPixel, + TString overlayType, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of max. propability" + << endl; + TH2F** hInputHistoArray = NULL; + TH1F** hOutputHistoArray = NULL; + + if (overlayType.Contains("Edge")) + { + hInputHistoArray = CurrentPixel->hEdgeOverlay; + hOutputHistoArray = CurrentPixel->hPixelEdgeMax; + + } + else if (overlayType.Contains("Maximum")) + { + hInputHistoArray = CurrentPixel->hMaxOverlay; + hOutputHistoArray = CurrentPixel->hPixelMax; + } + else + { + cout << endl << "Unknown Overlay Method-->aborting" << endl; + return; + } + + for ( int pulse_order = 0 ; + pulse_order < CurrentPixel->mMaxPulseOrder ; + pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...calculation of " + << "pulse order " << pulse_order; + // vector max_value_of to number of timeslices in Overlay Spectra + CalcMaxPropabilityOfSlice( + hInputHistoArray[pulse_order], + hOutputHistoArray[pulse_order], + verbosityLevel); + + CalcMaxPropabilityOfSlice( + hInputHistoArray[pulse_order], + hOutputHistoArray[pulse_order], + verbosityLevel); + } +} +// end of PlotMaxPropabilityPulse +//---------------------------------------------------------------------------- + +void +FitMaxPropabilityPulse( + TH1F* inputHisto, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << "...fit Landau in histograms" ; + if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ; + inputHisto->Fit("landau", "", "", -50, 250); + if (verbosityLevel > 2) cout << "...done" << endl; +} +// end of FitMaxPropabilityPuls +//---------------------------------------------------------------------------- + +Double_t MedianOfH1 ( + TH1* inputHisto + ) +{ + //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(i+1); + y[i] = inputHisto->GetBinContent(i+1); + } + Double_t median = TMath::Median(nbins,x,y); + delete [] x; + delete [] y; + return median; +} +// end of PlotMedianEachSliceOfPulse +//---------------------------------------------------------------------------- + +void +PlotMedianOfSlice( + Pixel* CurrentPixel, + TString overlayType, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of slice's Median" + << endl; + + TH2F** hInputHistoArray = NULL; + TH1F** hOutputHistoArray = NULL; + TH1* hTempHisto = NULL; + float median = 0; + + if (overlayType.Contains("Edge")) + { + hInputHistoArray = CurrentPixel->hEdgeOverlay; + hOutputHistoArray = CurrentPixel->hPixelEdgeMean; + + } + else if (overlayType.Contains("Maximum")) + { + hInputHistoArray = CurrentPixel->hMaxOverlay; + hOutputHistoArray = CurrentPixel->hPixelMedian; + } + else + { + cout << endl << "Unknown Overlay Method-->aborting" << endl; + return; + } + + for (int pulse_order = 0 ; + pulse_order < CurrentPixel->mMaxPulseOrder ; + pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...calculation of " + << "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 ); +// delete h1; + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + } +} +// end of PlotMedianEachSliceOfPulse +//---------------------------------------------------------------------------- + +void +PlotMeanOfSlice( + Pixel* CurrentPixel, + TString overlayType, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape of slice's Mean" + << endl; + + TH2F** hInputHistoArray = NULL; + TH1F** hOutputHistoArray = NULL; + TH1* hTempHisto = NULL; + float mean = 0; + + if (overlayType.Contains("Edge")) + { + hInputHistoArray = CurrentPixel->hEdgeOverlay; + hOutputHistoArray = CurrentPixel->hPixelEdgeMean; + + } + else if (overlayType.Contains("Maximum")) + { + hInputHistoArray = CurrentPixel->hMaxOverlay; + hOutputHistoArray = CurrentPixel->hPixelMean; + } + else + { + cout << endl << "Unknown Overlay Method-->aborting" << endl; + return; + } + + for (int pulse_order = 0 ; + pulse_order < CurrentPixel->mMaxPulseOrder ; + pulse_order ++) + { + if (verbosityLevel > 2) cout << "\t...calculation of " + << "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); + mean = hTempHisto->GetMean(); + + if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",TimeSlice,mean); + + hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, mean ); +// delete h1; + } + + if (verbosityLevel > 2) cout << "\t...done" << endl; + } +} +// end of CalcMeanOfSlice +//---------------------------------------------------------------------------- + +void +ExtractPulseTemplate( + Pixel* CurrentPixel, + TString overlayType, + int pulse_order, + int verbosityLevel + ) +{ + if (verbosityLevel > 2) cout << endl + << "...calculating pulse shape" + << " of max. propability" + << " of " + << overlayType + << " Overlay" + << endl; + TH2F* hInputHistoArray = NULL; + TH1F* hOutputMaxHistoArray = NULL; + TH1F* hOutputMeanHistoArray = NULL; + TH1F* hOutputMedianHistoArray = NULL; + TH1* hTempHisto = NULL; + float max_prop = 0; + float median = 0; + float mean = 0; + + if (verbosityLevel > 3) + { + cout << "...setting pointers to histogramm arrays "; + } + if (overlayType.Contains("Edge")) + { + hInputHistoArray = CurrentPixel->hEdgeOverlay[pulse_order]; + + //Maximum Propability of Slices + hOutputMaxHistoArray = CurrentPixel->hPixelEdgeMax[pulse_order]; + + //Mean of Slices + hOutputMeanHistoArray = CurrentPixel->hPixelEdgeMean[pulse_order]; + + //Median of Slices + hOutputMedianHistoArray = CurrentPixel->hPixelEdgeMedian[pulse_order]; + } + else if (overlayType.Contains("Maximum")) + { + hInputHistoArray = CurrentPixel->hMaxOverlay[pulse_order]; + + //Maximum Propability of Slices + hOutputMaxHistoArray = CurrentPixel->hPixelMax[pulse_order]; + + //Mean of Slices + hOutputMeanHistoArray = CurrentPixel->hPixelMean[pulse_order]; + + //Median of Slices + hOutputMedianHistoArray = CurrentPixel->hPixelMedian[pulse_order]; + } + else + { + cout << endl << "Unknown Overlay Method-->aborting" << endl; + return; + } + if (verbosityLevel > 3) + { + cout << "...done " << endl; + } + + if (verbosityLevel > 2) cout << "\t...calculation of " + << "pulse order " << pulse_order << endl; + + cout << "################ " << endl; + // vector max_value_of to number of timeslices in Overlay Spectra + cout << "### " << hInputHistoArray; + cout << "################ " << endl; + +// int nbins = hInputHistoArray->GetXaxis()->GetNbins(); +// cout << "###" << nbins << endl; + + +// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl; + + for (Int_t TimeSlice=1;TimeSlice<=300;TimeSlice++) + { + + hTempHisto = hInputHistoArray->ProjectionY("",TimeSlice,TimeSlice); + + if (verbosityLevel > 3) + { + cout << "\t\t...calculating maxProb of slice " << TimeSlice << endl; + } + max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() ); + + if (verbosityLevel > 3) + { + cout << "\t\t...calculating Median of slice " << TimeSlice << endl; + } + median = MedianOfH1(hTempHisto); + + if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << TimeSlice << endl; + mean = hTempHisto->GetMean(); + + if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << TimeSlice << ": " << max_prop << endl; + hOutputMaxHistoArray->SetBinContent(TimeSlice, max_prop ); + + if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << TimeSlice << ": " << mean << endl; + hOutputMeanHistoArray->SetBinContent(TimeSlice, mean ); + + if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << TimeSlice << ": " << median << endl; + hOutputMedianHistoArray->SetBinContent(TimeSlice, median ); +// delete h1; + + }//Loop over Timeslices +} +// end of PlotMaxPropabilityPulse +//---------------------------------------------------------------------------- + +bool +WritePixelTemplateToCsv( + Pixel* CurrentPixel, + TString path, + TString overlayMethod, + int order, + int verbosityLevel + ) +{ + path = CurrentPixel->CsvFileName( path, overlayMethod, order); + + TH1F** Max_histo_array = NULL; + TH1F** Median_histo_array = NULL; + TH1F** Mean_histo_array = NULL; + + if (overlayMethod.Contains("Maximum")) + { + Max_histo_array = CurrentPixel->hPixelMax; + Median_histo_array = CurrentPixel->hPixelMedian; + Mean_histo_array = CurrentPixel->hPixelMean; + } + else if (overlayMethod.Contains("Edge")) + { + Max_histo_array = CurrentPixel->hPixelMax; + Median_histo_array = CurrentPixel->hPixelMedian; + Mean_histo_array = CurrentPixel->hPixelMean; + } + else + { + cout << endl << "Unknown Overlay Method-->aborting" << endl; + return 1; + } + + + Int_t nbins = Max_histo_array[order]->GetXaxis()->GetNbins(); + + if (verbosityLevel > 0) + { + cout << "writing point-set to csv file: " ; + cout << path << endl; + cout << "...opening file" << endl; + } + if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl; + ofstream out; + out.open( path ); + out << "### point-set of a single photon pulse template" << endl; + out << "### template determined with pulse overlay at: " + << overlayMethod << endl; + out << "### Slice's Amplitude determined by calculating the " << endl + << "### value of maximum propability of slice -> AmplitudeMax " << endl + << "### mean of slice -> AmplitudeMean " << endl + << "### median of slice -> AmplitudeMedian " << endl + << "### for each slice" << endl; + out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl + << endl; + + out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl; + + for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++) + { + out << TimeSlice << "," ; + out << Max_histo_array[order]->GetBinContent(TimeSlice) << ","; + out << Mean_histo_array[order]->GetBinContent(TimeSlice) << ","; + out << Median_histo_array[order]->GetBinContent(TimeSlice) << endl; + } + out.close(); + if (verbosityLevel > 0) cout << "...file closed" << endl; + return 0; +}