source: fact/tools/rootmacros/PulseTemplates/templateextractors.C@ 13784

Last change on this file since 13784 was 13652, 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
  • Property svn:executable set to *
File size: 13.9 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 TString overlayType,
50 int verbosityLevel
51 )
52{
53 if (verbosityLevel > 2) cout << endl
54 << "...calculating pulse shape of max. propability"
55 << endl;
56 TH2F** hInputHistoArray = NULL;
57 TH1F** hOutputHistoArray = NULL;
58
59 if (overlayType.Contains("Edge"))
60 {
61 hInputHistoArray = CurrentPixel->hEdgeOverlay;
62 hOutputHistoArray = CurrentPixel->hPixelEdgeMax;
63
64 }
65 else if (overlayType.Contains("Maximum"))
66 {
67 hInputHistoArray = CurrentPixel->hMaxOverlay;
68 hOutputHistoArray = CurrentPixel->hPixelMax;
69 }
70 else
71 {
72 cout << endl << "Unknown Overlay Method-->aborting" << endl;
73 return;
74 }
75
76 for ( int pulse_order = 0 ;
77 pulse_order < CurrentPixel->mMaxPulseOrder ;
78 pulse_order ++)
79 {
80 if (verbosityLevel > 2) cout << "\t...calculation of "
81 << "pulse order " << pulse_order;
82 // vector max_value_of to number of timeslices in Overlay Spectra
83 CalcMaxPropabilityOfSlice(
84 hInputHistoArray[pulse_order],
85 hOutputHistoArray[pulse_order],
86 verbosityLevel);
87
88 CalcMaxPropabilityOfSlice(
89 hInputHistoArray[pulse_order],
90 hOutputHistoArray[pulse_order],
91 verbosityLevel);
92 }
93}
94// end of PlotMaxPropabilityPulse
95//----------------------------------------------------------------------------
96
97void
98FitMaxPropabilityPulse(
99 TH1F* inputHisto,
100 int verbosityLevel
101 )
102{
103 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
104 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
105 inputHisto->Fit("landau", "", "", -50, 250);
106 if (verbosityLevel > 2) cout << "...done" << endl;
107}
108// end of FitMaxPropabilityPuls
109//----------------------------------------------------------------------------
110
111Double_t MedianOfH1 (
112 TH1* inputHisto
113 )
114{
115 //compute the median for 1-d histogram h1
116 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
117 Double_t *x = new Double_t[nbins];
118 Double_t *y = new Double_t[nbins];
119 for (Int_t i=0;i<nbins;i++) {
120 x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
121 y[i] = inputHisto->GetBinContent(i+1);
122 }
123 Double_t median = TMath::Median(nbins,x,y);
124 delete [] x;
125 delete [] y;
126 return median;
127}
128// end of PlotMedianEachSliceOfPulse
129//----------------------------------------------------------------------------
130
131void
132PlotMedianOfSlice(
133 Pixel* CurrentPixel,
134 TString overlayType,
135 int verbosityLevel
136 )
137{
138 if (verbosityLevel > 2) cout << endl
139 << "...calculating pulse shape of slice's Median"
140 << endl;
141
142 TH2F** hInputHistoArray = NULL;
143 TH1F** hOutputHistoArray = NULL;
144 TH1* hTempHisto = NULL;
145 float median = 0;
146
147 if (overlayType.Contains("Edge"))
148 {
149 hInputHistoArray = CurrentPixel->hEdgeOverlay;
150 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
151
152 }
153 else if (overlayType.Contains("Maximum"))
154 {
155 hInputHistoArray = CurrentPixel->hMaxOverlay;
156 hOutputHistoArray = CurrentPixel->hPixelMedian;
157 }
158 else
159 {
160 cout << endl << "Unknown Overlay Method-->aborting" << endl;
161 return;
162 }
163
164 for (int pulse_order = 0 ;
165 pulse_order < CurrentPixel->mMaxPulseOrder ;
166 pulse_order ++)
167 {
168 if (verbosityLevel > 2) cout << "\t...calculation of "
169 << "pulse order " << pulse_order;
170
171 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
172
173 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
174
175 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
176 median = MedianOfH1(hTempHisto);
177
178 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
179
180 hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, median );
181// delete h1;
182 }
183
184 if (verbosityLevel > 2) cout << "\t...done" << endl;
185 }
186}
187// end of PlotMedianEachSliceOfPulse
188//----------------------------------------------------------------------------
189
190void
191PlotMeanOfSlice(
192 Pixel* CurrentPixel,
193 TString overlayType,
194 int verbosityLevel
195 )
196{
197 if (verbosityLevel > 2) cout << endl
198 << "...calculating pulse shape of slice's Mean"
199 << endl;
200
201 TH2F** hInputHistoArray = NULL;
202 TH1F** hOutputHistoArray = NULL;
203 TH1* hTempHisto = NULL;
204 float mean = 0;
205
206 if (overlayType.Contains("Edge"))
207 {
208 hInputHistoArray = CurrentPixel->hEdgeOverlay;
209 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
210
211 }
212 else if (overlayType.Contains("Maximum"))
213 {
214 hInputHistoArray = CurrentPixel->hMaxOverlay;
215 hOutputHistoArray = CurrentPixel->hPixelMean;
216 }
217 else
218 {
219 cout << endl << "Unknown Overlay Method-->aborting" << endl;
220 return;
221 }
222
223 for (int pulse_order = 0 ;
224 pulse_order < CurrentPixel->mMaxPulseOrder ;
225 pulse_order ++)
226 {
227 if (verbosityLevel > 2) cout << "\t...calculation of "
228 << "pulse order " << pulse_order;
229
230 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
231
232 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
233
234 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
235 mean = hTempHisto->GetMean();
236
237 if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",TimeSlice,mean);
238
239 hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, mean );
240// delete h1;
241 }
242
243 if (verbosityLevel > 2) cout << "\t...done" << endl;
244 }
245}
246// end of CalcMeanOfSlice
247//----------------------------------------------------------------------------
248
249void
250ExtractPulseTemplate(
251 Pixel* CurrentPixel,
252 TString overlayType,
253 int pulse_order,
254 int verbosityLevel
255 )
256{
257 if (verbosityLevel > 2) cout << endl
258 << "...calculating pulse shape"
259 << " of max. propability"
260 << " of "
261 << overlayType
262 << " Overlay"
263 << endl;
264 TH2F* hInputHistoArray = NULL;
265 TH1F* hOutputMaxHistoArray = NULL;
266 TH1F* hOutputMeanHistoArray = NULL;
267 TH1F* hOutputMedianHistoArray = NULL;
268 TH1* hTempHisto = NULL;
269 float max_prop = 0;
270 float median = 0;
271 float mean = 0;
272
273 if (verbosityLevel > 3)
274 {
275 cout << "...setting pointers to histogramm arrays ";
276 }
277 if (overlayType.Contains("Edge"))
278 {
279 hInputHistoArray = CurrentPixel->hEdgeOverlay[pulse_order];
280
281 //Maximum Propability of Slices
282 hOutputMaxHistoArray = CurrentPixel->hPixelEdgeMax[pulse_order];
283
284 //Mean of Slices
285 hOutputMeanHistoArray = CurrentPixel->hPixelEdgeMean[pulse_order];
286
287 //Median of Slices
288 hOutputMedianHistoArray = CurrentPixel->hPixelEdgeMedian[pulse_order];
289 }
290 else if (overlayType.Contains("Maximum"))
291 {
292 hInputHistoArray = CurrentPixel->hMaxOverlay[pulse_order];
293
294 //Maximum Propability of Slices
295 hOutputMaxHistoArray = CurrentPixel->hPixelMax[pulse_order];
296
297 //Mean of Slices
298 hOutputMeanHistoArray = CurrentPixel->hPixelMean[pulse_order];
299
300 //Median of Slices
301 hOutputMedianHistoArray = CurrentPixel->hPixelMedian[pulse_order];
302 }
303 else
304 {
305 cout << endl << "Unknown Overlay Method-->aborting" << endl;
306 return;
307 }
308 if (verbosityLevel > 3)
309 {
310 cout << "...done " << endl;
311 }
312
313 if (verbosityLevel > 2) cout << "\t...calculation of "
314 << "pulse order " << pulse_order << endl;
315
316 cout << "################ " << endl;
317 // vector max_value_of to number of timeslices in Overlay Spectra
318 cout << "### " << hInputHistoArray;
319 cout << "################ " << endl;
320
321// int nbins = hInputHistoArray->GetXaxis()->GetNbins();
322// cout << "###" << nbins << endl;
323
324
325// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
326
327 for (Int_t TimeSlice=1;TimeSlice<=300;TimeSlice++)
328 {
329
330 hTempHisto = hInputHistoArray->ProjectionY("",TimeSlice,TimeSlice);
331
332 if (verbosityLevel > 3)
333 {
334 cout << "\t\t...calculating maxProb of slice " << TimeSlice << endl;
335 }
336 max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
337
338 if (verbosityLevel > 3)
339 {
340 cout << "\t\t...calculating Median of slice " << TimeSlice << endl;
341 }
342 median = MedianOfH1(hTempHisto);
343
344 if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << TimeSlice << endl;
345 mean = hTempHisto->GetMean();
346
347 if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << TimeSlice << ": " << max_prop << endl;
348 hOutputMaxHistoArray->SetBinContent(TimeSlice, max_prop );
349
350 if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << TimeSlice << ": " << mean << endl;
351 hOutputMeanHistoArray->SetBinContent(TimeSlice, mean );
352
353 if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << TimeSlice << ": " << median << endl;
354 hOutputMedianHistoArray->SetBinContent(TimeSlice, median );
355// delete h1;
356
357 }//Loop over Timeslices
358}
359// end of PlotMaxPropabilityPulse
360//----------------------------------------------------------------------------
361
362bool
363WritePixelTemplateToCsv(
364 Pixel* CurrentPixel,
365 TString path,
366 TString overlayMethod,
367 int order,
368 int verbosityLevel
369 )
370{
371 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
372
373 TH1F** Max_histo_array = NULL;
374 TH1F** Median_histo_array = NULL;
375 TH1F** Mean_histo_array = NULL;
376
377 if (overlayMethod.Contains("Maximum"))
378 {
379 Max_histo_array = CurrentPixel->hPixelMax;
380 Median_histo_array = CurrentPixel->hPixelMedian;
381 Mean_histo_array = CurrentPixel->hPixelMean;
382 }
383 else if (overlayMethod.Contains("Edge"))
384 {
385 Max_histo_array = CurrentPixel->hPixelMax;
386 Median_histo_array = CurrentPixel->hPixelMedian;
387 Mean_histo_array = CurrentPixel->hPixelMean;
388 }
389 else
390 {
391 cout << endl << "Unknown Overlay Method-->aborting" << endl;
392 return 1;
393 }
394
395
396 Int_t nbins = Max_histo_array[order]->GetXaxis()->GetNbins();
397
398 if (verbosityLevel > 0)
399 {
400 cout << "writing point-set to csv file: " ;
401 cout << path << endl;
402 cout << "...opening file" << endl;
403 }
404 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
405 ofstream out;
406 out.open( path );
407 out << "### point-set of a single photon pulse template" << endl;
408 out << "### template determined with pulse overlay at: "
409 << overlayMethod << endl;
410 out << "### Slice's Amplitude determined by calculating the " << endl
411 << "### value of maximum propability of slice -> AmplitudeMax " << endl
412 << "### mean of slice -> AmplitudeMean " << endl
413 << "### median of slice -> AmplitudeMedian " << endl
414 << "### for each slice" << endl;
415 out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
416 << endl;
417
418 out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
419
420 for (int TimeSlice=1;TimeSlice<=nbins;TimeSlice++)
421 {
422 out << TimeSlice << "," ;
423 out << Max_histo_array[order]->GetBinContent(TimeSlice) << ",";
424 out << Mean_histo_array[order]->GetBinContent(TimeSlice) << ",";
425 out << Median_histo_array[order]->GetBinContent(TimeSlice) << endl;
426 }
427 out.close();
428 if (verbosityLevel > 0) cout << "...file closed" << endl;
429 return 0;
430}
Note: See TracBrowser for help on using the repository browser.