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

Last change on this file since 13851 was 13785, checked in by Jens Buss, 13 years ago
renamed temporary functions, set them to pointers
  • Property svn:executable set to *
File size: 13.7 KB
Line 
1#include <iostream>
2#include <fstream>
3
4#include "templateextractors.h"
5#include <stdlib.h>
6
7using namespace std;
8
9
10void
11CalcMaxPropabilityOfSlice(
12 TH2* inputHisto,
13 TH1* outputHisto,
14 int verbosityLevel
15 )
16{
17 if (verbosityLevel > 2) cout << endl
18 << "...calculating slieces value of max. propability"
19 << endl;
20
21 float max_value_of_slice = 0;
22 TH1D* hTempHisto = NULL;
23 int nbins = inputHisto->GetXaxis()->GetNbins();
24
25 if (verbosityLevel > 2) cout << "...generating projections" << endl;
26
27 for (Int_t TimeSlice = 1; TimeSlice <= nbins; TimeSlice++)
28 {
29 if (verbosityLevel > 2) cout << "...Timeslice: " << TimeSlice;
30
31 //put maximumvalue of every Y-projection of every timeslice into vector
32 hTempHisto = inputHisto->ProjectionY( "", TimeSlice,TimeSlice);
33 max_value_of_slice = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
34
35 if (verbosityLevel > 2) cout << " with max value "
36 << max_value_of_slice << endl;
37
38 outputHisto->SetBinContent(TimeSlice, max_value_of_slice );
39 }
40
41 if (verbosityLevel > 2) cout << "\t...done" << endl;
42}
43// end of CalcMaxPropabilityOfSlice
44//----------------------------------------------------------------------------
45
46
47void
48PlotMaxPropabilityPulse(
49 Pixel* CurrentPixel,
50 TString overlayType,
51 int verbosityLevel
52 )
53{
54 if (verbosityLevel > 2) cout << endl
55 << "...calculating pulse shape of max. propability"
56 << endl;
57 TH2F** hInputHistoArray = NULL;
58 TH1F** hOutputHistoArray = NULL;
59
60 if (overlayType.Contains("Edge"))
61 {
62 hInputHistoArray = CurrentPixel->hEdgeOverlay;
63 hOutputHistoArray = CurrentPixel->hPixelEdgeMax;
64
65 }
66 else if (overlayType.Contains("Maximum"))
67 {
68 hInputHistoArray = CurrentPixel->hMaxOverlay;
69 hOutputHistoArray = CurrentPixel->hPixelMax;
70 }
71 else
72 {
73 cout << endl << "Unknown Overlay Method-->aborting" << endl;
74 return;
75 }
76
77 for ( int pulse_order = 0 ;
78 pulse_order < CurrentPixel->mMaxPulseOrder ;
79 pulse_order ++)
80 {
81 if (verbosityLevel > 2) cout << "\t...calculation of "
82 << "pulse order " << pulse_order;
83 // vector max_value_of to number of timeslices in Overlay Spectra
84 CalcMaxPropabilityOfSlice(
85 hInputHistoArray[pulse_order],
86 hOutputHistoArray[pulse_order],
87 verbosityLevel);
88
89 CalcMaxPropabilityOfSlice(
90 hInputHistoArray[pulse_order],
91 hOutputHistoArray[pulse_order],
92 verbosityLevel);
93 }
94}
95// end of PlotMaxPropabilityPulse
96//----------------------------------------------------------------------------
97
98void
99FitMaxPropabilityPulse(
100 TH1F* inputHisto,
101 int verbosityLevel
102 )
103{
104 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
105 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
106 inputHisto->Fit("landau", "", "", -50, 250);
107 if (verbosityLevel > 2) cout << "...done" << endl;
108}
109// end of FitMaxPropabilityPuls
110//----------------------------------------------------------------------------
111
112Double_t MedianOfH1 (
113 TH1* inputHisto
114 )
115{
116 //compute the median for 1-d histogram h1
117 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
118 Double_t *x = new Double_t[nbins];
119 Double_t *y = new Double_t[nbins];
120 for (Int_t i=0;i<nbins;i++) {
121 x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
122 y[i] = inputHisto->GetBinContent(i+1);
123 }
124 Double_t median = TMath::Median(nbins,x,y);
125 delete [] x;
126 delete [] y;
127 return median;
128}
129// end of PlotMedianEachSliceOfPulse
130//----------------------------------------------------------------------------
131
132void
133PlotMedianOfSlice(
134 Pixel* CurrentPixel,
135 TString overlayType,
136 int verbosityLevel
137 )
138{
139 if (verbosityLevel > 2) cout << endl
140 << "...calculating pulse shape of slice's Median"
141 << endl;
142
143 TH2F** hInputHistoArray = NULL;
144 TH1F** hOutputHistoArray = NULL;
145 TH1* hTempHisto = NULL;
146 float median = 0;
147
148 if (overlayType.Contains("Edge"))
149 {
150 hInputHistoArray = CurrentPixel->hEdgeOverlay;
151 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
152
153 }
154 else if (overlayType.Contains("Maximum"))
155 {
156 hInputHistoArray = CurrentPixel->hMaxOverlay;
157 hOutputHistoArray = CurrentPixel->hPixelMedian;
158 }
159 else
160 {
161 cout << endl << "Unknown Overlay Method-->aborting" << endl;
162 return;
163 }
164
165 for (int pulse_order = 0 ;
166 pulse_order < CurrentPixel->mMaxPulseOrder ;
167 pulse_order ++)
168 {
169 if (verbosityLevel > 2) cout << "\t...calculation of "
170 << "pulse order " << pulse_order;
171
172 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
173
174 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
175
176 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
177 median = MedianOfH1(hTempHisto);
178
179 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",TimeSlice,median);
180
181 hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, median );
182// delete h1;
183 }
184
185 if (verbosityLevel > 2) cout << "\t...done" << endl;
186 }
187}
188// end of PlotMedianEachSliceOfPulse
189//----------------------------------------------------------------------------
190
191void
192PlotMeanOfSlice(
193 Pixel* CurrentPixel,
194 TString overlayType,
195 int verbosityLevel
196 )
197{
198 if (verbosityLevel > 2) cout << endl
199 << "...calculating pulse shape of slice's Mean"
200 << endl;
201
202 TH2F** hInputHistoArray = NULL;
203 TH1F** hOutputHistoArray = NULL;
204 TH1* hTempHisto = NULL;
205 float mean = 0;
206
207 if (overlayType.Contains("Edge"))
208 {
209 hInputHistoArray = CurrentPixel->hEdgeOverlay;
210 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
211
212 }
213 else if (overlayType.Contains("Maximum"))
214 {
215 hInputHistoArray = CurrentPixel->hMaxOverlay;
216 hOutputHistoArray = CurrentPixel->hPixelMean;
217 }
218 else
219 {
220 cout << endl << "Unknown Overlay Method-->aborting" << endl;
221 return;
222 }
223
224 for (int pulse_order = 0 ;
225 pulse_order < CurrentPixel->mMaxPulseOrder ;
226 pulse_order ++)
227 {
228 if (verbosityLevel > 2) cout << "\t...calculation of "
229 << "pulse order " << pulse_order;
230
231 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
232
233 for (Int_t TimeSlice=1;TimeSlice<=nbins;TimeSlice++) {
234
235 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",TimeSlice,TimeSlice);
236 mean = hTempHisto->GetMean();
237
238 if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",TimeSlice,mean);
239
240 hOutputHistoArray[pulse_order]->SetBinContent(TimeSlice, mean );
241// delete h1;
242 }
243
244 if (verbosityLevel > 2) cout << "\t...done" << endl;
245 }
246}
247// end of CalcMeanOfSlice
248//----------------------------------------------------------------------------
249
250void
251ExtractPulseTemplate(
252 Pixel* CurrentPixel,
253 TString overlayType,
254 int pulse_order,
255 int verbosityLevel
256 )
257{
258 if (verbosityLevel > 2) cout << endl
259 << "...calculating pulse shape"
260 << " of max. propability"
261 << " of "
262 << overlayType
263 << " Overlay"
264 << endl;
265 TH2F* hInputHisto = NULL;
266 TH1F* hOutputMaxHisto = NULL;
267 TH1F* hOutputMeanHisto = NULL;
268 TH1F* hOutputMedianHisto = NULL;
269 TH1* hTempHisto = NULL;
270 float max_prop = 0;
271 float median = 0;
272 float mean = 0;
273
274 if (verbosityLevel > 3)
275 {
276 cout << "...setting pointers to histogramm arrays ";
277 cout << " for " << overlayType << "Overlay" << endl;
278 }
279 if (overlayType.Contains("Edge"))
280 {
281 hInputHisto = (CurrentPixel->hEdgeOverlay[pulse_order]);
282
283 //Maximum Propability of Slices
284 hOutputMaxHisto = (CurrentPixel->hPixelEdgeMax[pulse_order]);
285
286 //Mean of Slices
287 hOutputMeanHisto = (CurrentPixel->hPixelEdgeMean[pulse_order]);
288
289 //Median of Slices
290 hOutputMedianHisto = (CurrentPixel->hPixelEdgeMedian[pulse_order]);
291 }
292 else if (overlayType.Contains("Maximum"))
293 {
294 hInputHisto = (CurrentPixel->hMaxOverlay[pulse_order]);
295
296 //Maximum Propability of Slices
297 hOutputMaxHisto = (CurrentPixel->hPixelMax[pulse_order]);
298
299 //Mean of Slices
300 hOutputMeanHisto = (CurrentPixel->hPixelMean[pulse_order]);
301
302 //Median of Slices
303 hOutputMedianHisto = (CurrentPixel->hPixelMedian[pulse_order]);
304 }
305 else
306 {
307 cout << endl << overlayType << "Unknown Overlay Method-->aborting" << endl;
308 return;
309 }
310 if (verbosityLevel > 3)
311 {
312 cout << "...done " << endl;
313 }
314
315 if (verbosityLevel > 2)
316 {
317 cout << "\t...calculation of "
318 << "pulse order " << pulse_order << endl;
319 }
320
321// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
322
323 for (Int_t TimeSlice=1;TimeSlice<=300;TimeSlice++)
324 {
325
326 hTempHisto = hInputHisto->ProjectionY("",TimeSlice,TimeSlice);
327
328 if (verbosityLevel > 3)
329 {
330 cout << "\t\t...calculating maxProb of slice " << TimeSlice << endl;
331 }
332 max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
333
334 if (verbosityLevel > 3)
335 {
336 cout << "\t\t...calculating Median of slice " << TimeSlice << endl;
337 }
338 median = MedianOfH1(hTempHisto);
339
340 if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << TimeSlice << endl;
341 mean = hTempHisto->GetMean();
342
343 if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << TimeSlice << ": " << max_prop << endl;
344 hOutputMaxHisto->SetBinContent(TimeSlice, max_prop );
345
346 if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << TimeSlice << ": " << mean << endl;
347 hOutputMeanHisto->SetBinContent(TimeSlice, mean );
348
349 if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << TimeSlice << ": " << median << endl;
350 hOutputMedianHisto->SetBinContent(TimeSlice, median );
351 delete hTempHisto;
352
353 }//Loop over Timeslices
354}
355// end of PlotMaxPropabilityPulse
356//----------------------------------------------------------------------------
357
358bool
359WritePixelTemplateToCsv(
360 Pixel* CurrentPixel,
361 TString path,
362 TString overlayMethod,
363 int order,
364 int verbosityLevel
365 )
366{
367// TSystem this_system();
368// this_system.cd(path);
369// cout << this_system.pwd();
370// this_system.mkdir("CSV");
371 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
372
373 TH1F* Max_histo = NULL;
374 TH1F* Median_histo = NULL;
375 TH1F* Mean_histo = NULL;
376
377 if (overlayMethod.Contains("Maximum"))
378 {
379 Max_histo = CurrentPixel->hPixelMax[order];
380 Median_histo = CurrentPixel->hPixelMedian[order];
381 Mean_histo = CurrentPixel->hPixelMean[order];
382 }
383 else if (overlayMethod.Contains("Edge"))
384 {
385 Max_histo = CurrentPixel->hPixelMax[order];
386 Median_histo = CurrentPixel->hPixelMedian[order];
387 Mean_histo = CurrentPixel->hPixelMean[order];
388 }
389 else
390 {
391 cout << endl << "Unknown Overlay Method-->aborting" << endl;
392 return 1;
393 }
394
395
396 Int_t nbins = Max_histo->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->GetBinContent(TimeSlice) << ",";
424 out << Mean_histo->GetBinContent(TimeSlice) << ",";
425 out << Median_histo->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.