source: fact/tools/rootmacros/PulseTemplates/templateextractors.cpp@ 13736

Last change on this file since 13736 was 13637, checked in by Jens Buss, 13 years ago
initial commit: functions for extracting information from TH2 and putting into TH1 to compute a Pulse Template
File size: 7.6 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include "templateextractors.h"
5#include <stdlib.h>
6using namespace std;
7
8
9void
10CalcMaxPropabilityOfSlice(
11 TH2* inputHisto,
12 TH1* outputHisto,
13 int verbosityLevel
14 )
15{
16 if (verbosityLevel > 2) cout << endl
17 << "...calculating slieces value of max. propability"
18 << endl;
19
20 float max_value_of_slice = 0;
21 TH1D* hTempHisto = NULL;
22 int nbins = inputHisto->GetXaxis()->GetNbins();
23
24 if (verbosityLevel > 2) cout << "...generating projections" << endl;
25
26 for (Int_t TimeSlice = 1; TimeSlice <= nbins; TimeSlice++)
27 {
28 if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
29
30 //put maximumvalue of every Y-projection of every timeslice into vector
31 hTempHisto = inputHisto->ProjectionY( "", TimeSlice,TimeSlice);
32 max_value_of_slice = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
33
34 if (verbosityLevel > 2) cout << " with max value "
35 << max_value_of_slice << endl;
36
37 outputHisto->SetBinContent(TimeSlice, max_value_of_slice );
38 }
39
40 if (verbosityLevel > 2) cout << "\t...done" << endl;
41}
42// end of CalcMaxPropabilityOfSlice
43//----------------------------------------------------------------------------
44
45
46void
47PlotMaxPropabilityPulse(
48 Pixel* CurrentPixel,
49 int verbosityLevel
50 )
51{
52 if (verbosityLevel > 2) cout << endl
53 << "...calculating pulse shape of max. propability"
54 << endl;
55 for ( int pulse_order = 0 ;
56 pulse_order < CurrentPixel->mMaxPulseOrder ;
57 pulse_order ++)
58 {
59 if (verbosityLevel > 2) cout << "\t...calculation of "
60 << "pulse order " << pulse_order;
61 // vector max_value_of to number of timeslices in Overlay Spectra
62 CalcMaxPropabilityOfSlice(
63 CurrentPixel->hMaxOverlay[pulse_order],
64 CurrentPixel->hPixelMax[pulse_order],
65 verbosityLevel);
66
67 CalcMaxPropabilityOfSlice(
68 CurrentPixel->hEdgeOverlay[pulse_order],
69 CurrentPixel->hPixelEdgeMax[pulse_order],
70 verbosityLevel);
71 }
72}
73// end of PlotMaxPropabilityPulse
74//----------------------------------------------------------------------------
75
76void
77FitMaxPropabilityPulse(
78 TH1F* inputHisto,
79 int verbosityLevel
80 )
81{
82 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
83 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
84 inputHisto->Fit("landau", "", "", -50, 250);
85 if (verbosityLevel > 2) cout << "...done" << endl;
86}
87// end of FitMaxPropabilityPuls
88//----------------------------------------------------------------------------
89
90Double_t MedianOfH1 (
91 TH1* inputHisto
92 )
93{
94 //compute the median for 1-d histogram h1
95 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
96 Double_t *x = new Double_t[nbins];
97 Double_t *y = new Double_t[nbins];
98 for (Int_t i=0;i<nbins;i++) {
99 x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
100 y[i] = inputHisto->GetBinContent(i+1);
101 }
102 Double_t median = TMath::Median(nbins,x,y);
103 delete [] x;
104 delete [] y;
105 return median;
106}
107// end of PlotMedianEachSliceOfPulse
108//----------------------------------------------------------------------------
109
110void
111CalcMedianOfSlice(
112 Pixel* CurrentPixel,
113 TString overlayType,
114 int verbosityLevel
115 )
116{
117 if (verbosityLevel > 2) cout << endl
118 << "...calculating pulse shape of slice's Median"
119 << endl;
120
121 TH2F** hInputHistoArray = NULL;
122 TH1F** hOutputHistoArray = NULL;
123 TH1* hTempHisto = NULL;
124 float median = 0;
125
126 if (overlayType.Contains("Edge"))
127 {
128 hInputHistoArray = CurrentPixel->hEdgeOverlay;
129 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
130
131 }
132 else if (overlayType.Contains("Max"))
133 {
134 hInputHistoArray = CurrentPixel->hMaxOverlay;
135 hOutputHistoArray = CurrentPixel->hPixelMedian;
136 }
137
138 for (int pulse_order = 0 ;
139 pulse_order < CurrentPixel->mMaxPulseOrder ;
140 pulse_order ++)
141 {
142 if (verbosityLevel > 2) cout << "\t...calculation of "
143 << "pulse order " << pulse_order;
144
145 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
146
147 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
148
149 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
150 median = MedianOfH1(hTempHisto);
151
152 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
153
154 hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, median );
155// delete h1;
156 }
157
158 if (verbosityLevel > 2) cout << "\t...done" << endl;
159 }
160}
161// end of PlotMedianEachSliceOfPulse
162//----------------------------------------------------------------------------
163
164bool
165WritePixelTemplateToCsv(
166 Pixel* CurrentPixel,
167 TString path,
168 TString overlayMethod,
169 int order,
170 int verbosityLevel
171 )
172{
173 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
174
175 TH1F** Max_histo_array = NULL;
176 TH1F** Median_histo_array = NULL;
177 TH1F** Mean_histo_array = NULL;
178
179 if (overlayMethod.Contains("Maximum"))
180 {
181 Max_histo_array = CurrentPixel->hPixelMax;
182 Median_histo_array = CurrentPixel->hPixelMedian;
183 Mean_histo_array = CurrentPixel->hPixelMean;
184 }
185 else if (overlayMethod.Contains("Edge"))
186 {
187 Max_histo_array = CurrentPixel->hPixelMax;
188 Median_histo_array = CurrentPixel->hPixelMedian;
189 Mean_histo_array = CurrentPixel->hPixelMean;
190 }
191 else
192 {
193 cout << endl << "Unknown Overlay Method-->aborting" << endl;
194 return 1;
195 }
196
197
198 Int_t nbins = Max_histo_array[order]->GetXaxis()->GetNbins();
199
200 if (verbosityLevel > 0)
201 {
202 cout << "writing point-set to csv file: " ;
203 cout << path << endl;
204 cout << "...opening file" << endl;
205 }
206 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
207 ofstream out;
208 out.open( path );
209 out << "### point-set of a single photon pulse template" << endl;
210 out << "### template determined with pulse overlay at: "
211 << overlayMethod << endl;
212 out << "### Slice's Amplitude determined by calculating the " << endl
213 << "### value of maximum propability of slice -> AmplitudeMax " << endl
214 << "### mean of slice -> AmplitudeMean " << endl
215 << "### median of slice -> AmplitudeMedian " << endl
216 << "### for each slice" << endl;
217 out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
218 << endl;
219
220 out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
221
222 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
223 {
224 out << TimeSlice << "," ;
225 out << Max_histo_array[order]->GetBinContent(TimeSlice) << ",";
226 out << Mean_histo_array[order]->GetBinContent(TimeSlice) << ",";
227 out << Median_histo_array[order]->GetBinContent(TimeSlice) << endl;
228 }
229 out.close();
230 if (verbosityLevel > 0) cout << "...file closed" << endl;
231 return 0;
232}
Note: See TracBrowser for help on using the repository browser.