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

Last change on this file since 14708 was 14707, checked in by Jens Buss, 13 years ago
implemented error calculation for median
  • Property svn:executable set to *
File size: 20.8 KB
Line 
1#include <iostream>
2#include <fstream>
3#include <stdlib.h>
4
5#include "templateextractors.h"
6
7using namespace std;
8
9
10void
11CalcMaxPropabilityOfSlice(
12 TH2* inputHisto,
13 TH1* outputHisto,
14 int verbosityLevel
15 )
16{
17 if (verbosityLevel > 2)
18 cout << endl
19 << "...calculating slieces value of max. propability"
20 << endl;
21
22 float sliceMax = 0;
23 TH1D* hTemp = NULL;
24 int nbins = inputHisto->GetXaxis()->GetNbins();
25
26 if (verbosityLevel > 2)
27 cout << "...generating projections" << endl;
28
29 for (Int_t slice = 1; slice <= nbins; slice++)
30 {
31 if (verbosityLevel > 2)
32 cout << "...slice: " << slice;
33
34 //put maximumvalue of every Y-projection of every slice into vector
35 hTemp = inputHisto->ProjectionY( "", slice,slice);
36 sliceMax = hTemp->GetBinCenter( hTemp->GetMaximumBin() );
37
38 if (verbosityLevel > 2)
39 cout << " with max value "
40 << sliceMax << endl;
41
42 outputHisto->SetBinContent(slice, sliceMax );
43 }
44
45 if (verbosityLevel > 2)
46 cout << "\t...done" << endl;
47}
48// end of CalcMaxPropabilityOfSlice
49//----------------------------------------------------------------------------
50
51
52void
53PlotMaxPropabilityPulse(
54 Pixel* CurrentPixel,
55 TString overlayType,
56 int verbosityLevel
57 )
58{
59 if (verbosityLevel > 2) cout << endl
60 << "...calculating pulse shape of max. propability"
61 << endl;
62 TH2F** hInputHistoArray = NULL;
63 TH1F** hOutputHistoArray = NULL;
64
65 if (overlayType.Contains("Edge"))
66 {
67 hInputHistoArray = CurrentPixel->hEdgeOverlay;
68 hOutputHistoArray = CurrentPixel->hPixelEdgeMax;
69
70 }
71 else if (overlayType.Contains("Maximum"))
72 {
73 hInputHistoArray = CurrentPixel->hMaxOverlay;
74 hOutputHistoArray = CurrentPixel->hPixelMax;
75 }
76 else
77 {
78 cout << endl << "Unknown Overlay Method-->aborting" << endl;
79 return;
80 }
81
82 for ( int pulse_order = 0 ;
83 pulse_order < CurrentPixel->mMaxPulseOrder ;
84 pulse_order ++)
85 {
86 if (verbosityLevel > 2) cout << "\t...calculation of "
87 << "pulse order " << pulse_order;
88 // vector max_value_of to number of slices in Overlay Spectra
89 CalcMaxPropabilityOfSlice(
90 hInputHistoArray[pulse_order],
91 hOutputHistoArray[pulse_order],
92 verbosityLevel);
93
94 CalcMaxPropabilityOfSlice(
95 hInputHistoArray[pulse_order],
96 hOutputHistoArray[pulse_order],
97 verbosityLevel);
98 }
99}
100// end of PlotMaxPropabilityPulse
101//----------------------------------------------------------------------------
102
103void
104FitMaxPropabilityPulse(
105 TH1F* inputHisto,
106 int verbosityLevel
107 )
108{
109 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
110 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
111 inputHisto->Fit("landau", "", "", -50, 250);
112 if (verbosityLevel > 2) cout << "...done" << endl;
113}
114// end of FitMaxPropabilityPuls
115//----------------------------------------------------------------------------
116
117Double_t MedianOfH1 (
118 TH1* inputHisto
119 )
120{
121 //compute the median for 1-d histogram h1
122 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
123 Double_t *x = new Double_t[nbins];
124 Double_t *y = new Double_t[nbins];
125 for (Int_t i=0;i<nbins;i++) {
126 x[i] = inputHisto->GetXaxis()->GetBinCenter(i+1);
127 y[i] = inputHisto->GetBinContent(i+1);
128 }
129 Double_t median = TMath::Median(nbins,x,y);
130 delete [] x;
131 delete [] y;
132 return median;
133}
134// end of PlotMedianEachSliceOfPulse
135//----------------------------------------------------------------------------
136
137Double_t MedianOfH1withRndSlices (
138 TH1* inputHisto, //histogram from which median will be calculated
139 vector<int>* position //array with random slice numbers
140 )
141{
142 //compute the median for 1-d histogram h1
143 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
144 Double_t *x = new Double_t[nbins];
145 Double_t *y = new Double_t[nbins];
146 for (Int_t i=0;i<nbins;i++) {
147 x[i] = inputHisto->GetXaxis()->GetBinCenter(position->at(i) );
148 y[i] = inputHisto->GetBinContent(position->at(i) );
149 }
150 Double_t median = TMath::Median(nbins,x,y);
151 delete [] x;
152 delete [] y;
153 return median;
154}
155// end of PlotMedianEachSliceOfPulse
156//----------------------------------------------------------------------------
157
158void
159PlotMedianOfSlice(
160 Pixel* CurrentPixel,
161 TString overlayType,
162 int verbosityLevel
163 )
164{
165 if (verbosityLevel > 2) cout << endl
166 << "...calculating pulse shape of slice's Median"
167 << endl;
168
169 TH2F** hInputHistoArray = NULL;
170 TH1F** hOutputHistoArray = NULL;
171// TH1* hTempHisto = NULL;
172// float median = 0;
173
174 if (overlayType.Contains("Edge"))
175 {
176 hInputHistoArray = CurrentPixel->hEdgeOverlay;
177 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
178
179 }
180 else if (overlayType.Contains("Maximum"))
181 {
182 hInputHistoArray = CurrentPixel->hMaxOverlay;
183 hOutputHistoArray = CurrentPixel->hPixelMedian;
184 }
185 else
186 {
187 cout << endl << "Unknown Overlay Method-->aborting" << endl;
188 return;
189 }
190
191 for (int pulse_order = 0 ;
192 pulse_order < CurrentPixel->mMaxPulseOrder ;
193 pulse_order ++)
194 {
195 if (verbosityLevel > 2) cout << "\t...calculation of "
196 << "pulse order " << pulse_order;
197
198 MedianOfTH2Slices(
199 hInputHistoArray[pulse_order],
200 hOutputHistoArray[pulse_order],
201 verbosityLevel
202 );
203
204 ErrorMedianOfTH2Slices(
205 hInputHistoArray[pulse_order],
206 hOutputHistoArray[pulse_order],
207 5, //numIterations
208 verbosityLevel
209 );
210
211// Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
212
213// for (Int_t slice=1;slice<=nbins;slice++) {
214
215// hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
216// median = MedianOfH1(hTempHisto);
217
218// if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
219
220// hOutputHistoArray[pulse_order]->SetBinContent(slice, median );
221//// delete h1;
222// }
223
224 if (verbosityLevel > 2) cout << "\t...done" << endl;
225 }
226}
227// end of PlotMedianEachSliceOfPulse
228//----------------------------------------------------------------------------
229
230void
231MedianOfTH2Slices(
232 TH2* inputHisto,
233 TH1* outputHisto,
234 int verbosityLevel
235 )
236{
237 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
238 float median = 0;
239
240 for (Int_t slice=1;slice<=nbins;slice++) {
241
242 median = MedianOfH1( inputHisto->ProjectionY("",slice,slice) );
243
244 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
245
246 outputHisto->SetBinContent(slice, median );
247// delete h1;
248 }
249 return;
250}
251// end of MedianOfTH2Slices
252//----------------------------------------------------------------------------
253
254void
255ErrorMedianOfTH2Slices(
256 TH2* inputHisto,
257 TH1* outputHisto,
258 int numIterations,
259 int verbosityLevel
260 )
261{
262 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
263
264// float MedianOfSliceMean[nbins];
265// float MedianOfSliceRMS[nbins];
266
267 for (Int_t slice=1;slice<=nbins;slice++) {
268
269 float median[numIterations];
270 int sample_min = inputHisto->GetXaxis()->GetFirst();
271 int sample_max = inputHisto->GetXaxis()->GetLast();
272 int sample_size = sample_max - sample_min;
273 vector<int> rndList;
274
275
276 Sample sample(sample_size);
277 for (int i = 0; i < numIterations; i++)
278 {
279 sample.BootstrapSample(rndList, sample_min, sample_max, sample_size);
280 median[i] = MedianOfH1withRndSlices(
281 inputHisto->ProjectionY("",slice,slice),
282 rndList
283 );
284
285 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
286 }
287// MedianOfSliceMean[slice] = TMath::Mean(numIterations, median);
288// MedianOfSliceRMS[slice] = TMath::RMS(numIterations, median);
289// outputHisto->SetBinError(slice, MedianOfSliceRMS);
290 outputHisto->SetBinError(slice, RMS(numIterations, median) );
291
292// outputHisto->SetBinContent(slice, median );
293// delete h1;
294 }
295 return;
296}
297// end of PlotMedianEachSliceOfPulse
298//----------------------------------------------------------------------------
299
300void
301PlotMeanOfSlice(
302 Pixel* CurrentPixel,
303 TString overlayType,
304 int verbosityLevel
305 )
306{
307 if (verbosityLevel > 2) cout << endl
308 << "...calculating pulse shape of slice's Mean"
309 << endl;
310
311 TH2F** hInputHistoArray = NULL;
312 TH1F** hOutputHistoArray = NULL;
313 TH1* hTempHisto = NULL;
314 float mean = 0;
315
316 if (overlayType.Contains("Edge"))
317 {
318 hInputHistoArray = CurrentPixel->hEdgeOverlay;
319 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
320
321 }
322 else if (overlayType.Contains("Maximum"))
323 {
324 hInputHistoArray = CurrentPixel->hMaxOverlay;
325 hOutputHistoArray = CurrentPixel->hPixelMean;
326 }
327 else
328 {
329 cout << endl << "Unknown Overlay Method-->aborting" << endl;
330 return;
331 }
332
333 for (int pulse_order = 0 ;
334 pulse_order < CurrentPixel->mMaxPulseOrder ;
335 pulse_order ++)
336 {
337 if (verbosityLevel > 2) cout << "\t...calculation of "
338 << "pulse order " << pulse_order;
339
340 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
341
342 for (Int_t slice=1;slice<=nbins;slice++) {
343
344 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
345 mean = hTempHisto->GetMean();
346
347 if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",slice,mean);
348
349 hOutputHistoArray[pulse_order]->SetBinContent(slice, mean );
350// delete h1;
351 }
352
353 if (verbosityLevel > 2) cout << "\t...done" << endl;
354 }
355}
356// end of CalcMeanOfSlice
357//----------------------------------------------------------------------------
358
359void
360ExtractPulseTemplate(
361 Pixel* CurrentPixel,
362 TString overlayType,
363 int pulse_order,
364 int verbosityLevel
365 )
366{
367 if (verbosityLevel > 2) cout << endl
368 << "...calculating pulse shape"
369 << " of max. propability"
370 << " of "
371 << overlayType
372 << " Overlay"
373 << endl;
374 TH2F* hInputHisto = NULL;
375 TH1F* hOutputMaxHisto = NULL;
376 TH1F* hOutputMeanHisto = NULL;
377 TH1F* hOutputMedianHisto = NULL;
378 TH1* hTempHisto = NULL;
379 float max_prop = 0;
380 float median = 0;
381 float mean = 0;
382
383 if (verbosityLevel > 3)
384 {
385 cout << "...setting pointers to histogramm arrays ";
386 cout << " for " << overlayType << "Overlay" << endl;
387 }
388 if (overlayType.Contains("Edge"))
389 {
390 hInputHisto = (CurrentPixel->hEdgeOverlay[pulse_order]);
391
392 //Maximum Propability of Slices
393 hOutputMaxHisto = (CurrentPixel->hPixelEdgeMax[pulse_order]);
394
395 //Mean of Slices
396 hOutputMeanHisto = (CurrentPixel->hPixelEdgeMean[pulse_order]);
397
398 //Median of Slices
399 hOutputMedianHisto = (CurrentPixel->hPixelEdgeMedian[pulse_order]);
400 }
401 else if (overlayType.Contains("Maximum"))
402 {
403 hInputHisto = (CurrentPixel->hMaxOverlay[pulse_order]);
404
405 //Maximum Propability of Slices
406 hOutputMaxHisto = (CurrentPixel->hPixelMax[pulse_order]);
407
408 //Mean of Slices
409 hOutputMeanHisto = (CurrentPixel->hPixelMean[pulse_order]);
410
411 //Median of Slices
412 hOutputMedianHisto = (CurrentPixel->hPixelMedian[pulse_order]);
413 }
414 else
415 {
416 cout << endl << overlayType << "Unknown Overlay Method-->aborting" << endl;
417 return;
418 }
419 if (verbosityLevel > 3)
420 {
421 cout << "...done " << endl;
422 }
423
424 if (verbosityLevel > 2)
425 {
426 cout << "\t...calculation of "
427 << "pulse order " << pulse_order << endl;
428 }
429
430// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
431
432 for (Int_t slice=1;slice<=300;slice++)
433 {
434
435 hTempHisto = hInputHisto->ProjectionY("",slice,slice);
436
437 if (verbosityLevel > 3)
438 {
439 cout << "\t\t...calculating maxProb of slice " << slice << endl;
440 }
441 max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
442
443 if (verbosityLevel > 3)
444 {
445 cout << "\t\t...calculating Median of slice " << slice << endl;
446 }
447 median = MedianOfH1(hTempHisto);
448
449 if (verbosityLevel > 4) cout << "\t\t...calculating Mean of slice " << slice << endl;
450 mean = hTempHisto->GetMean();
451
452 if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
453 hOutputMaxHisto->SetBinContent(slice, max_prop );
454
455 if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
456 hOutputMeanHisto->SetBinContent(slice, mean );
457
458 if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
459 hOutputMedianHisto->SetBinContent(slice, median );
460 delete hTempHisto;
461
462 }//Loop over slices
463}
464// end of PlotMaxPropabilityPulse
465//----------------------------------------------------------------------------
466
467bool
468WritePixelTemplateToCsv(
469 Pixel* CurrentPixel,
470 TString path,
471 TString overlayMethod,
472 int order,
473 int verbosityLevel
474 )
475{
476// TSystem this_system();
477// this_system.cd(path);
478// cout << this_system.pwd();
479// this_system.mkdir("CSV");
480 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
481
482 TH1F* Max_histo = NULL;
483 TH1F* Median_histo = NULL;
484 TH1F* Mean_histo = NULL;
485
486 if (overlayMethod.Contains("Maximum"))
487 {
488 Max_histo = CurrentPixel->hPixelMax[order];
489 Median_histo = CurrentPixel->hPixelMedian[order];
490 Mean_histo = CurrentPixel->hPixelMean[order];
491 }
492 else if (overlayMethod.Contains("Edge"))
493 {
494 Max_histo = CurrentPixel->hPixelMax[order];
495 Median_histo = CurrentPixel->hPixelMedian[order];
496 Mean_histo = CurrentPixel->hPixelMean[order];
497 }
498 else
499 {
500 cout << endl << "Unknown Overlay Method-->aborting" << endl;
501 return 1;
502 }
503
504
505 Int_t nbins = Max_histo->GetXaxis()->GetNbins();
506
507 if (verbosityLevel > 2)
508 {
509 cout << "writing point-set to csv file: " ;
510 cout << path << endl;
511 cout << "...opening file" << endl;
512 }
513 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
514 ofstream out;
515 out.open( path );
516 out << "### point-set of a single photon pulse template" << endl;
517 out << "### template determined with pulse overlay at: "
518 << overlayMethod << endl;
519 out << "### Slice's Amplitude determined by calculating the " << endl
520 << "### value of maximum propability of slice -> AmplitudeMax " << endl
521 << "### mean of slice -> AmplitudeMean " << endl
522 << "### median of slice -> AmplitudeMedian " << endl
523 << "### for each slice" << endl;
524 out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
525 << endl;
526
527 out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
528
529 for (int slice=1;slice<=nbins;slice++)
530 {
531 out << slice << "," ;
532 out << Max_histo->GetBinContent(slice) << ",";
533 out << Mean_histo->GetBinContent(slice) << ",";
534 out << Median_histo->GetBinContent(slice) << endl;
535 }
536 out.close();
537 if (verbosityLevel > 2) cout << "...csv file closed" << endl;
538 return 0;
539}
540
541void
542FitMaxPropabilityPuls(
543 TH1F* hMaximumTemp,
544 int verbosityLevel
545 )
546 {
547 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
548 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
549 hMaximumTemp->Fit("landau", "", "", -50, 250);
550 if (verbosityLevel > 2) cout << "...done" << endl;
551 }
552
553//void
554//FitFallingEdge(
555// TString name,
556// TH1F* histo,
557// double xMin,
558// double xMax,
559// double* parameters
560// )
561//{
562// TF1* polExpFit = new TF1(name, PolExp, xMin, xMax, 3 );
563// polExpFit->SetParNames("Pol0", "Slope", "Shift");
564// polExpFit->SetLineColor(kGreen);
565// histo->Fit(polExpFit, "+RWM");
566// polExpFit->GetParameters(parameters);
567//}
568
569//void
570//FitRisingEdge(
571// TString name,
572// TH1F* histo,
573// double xMin,
574// double xMax,
575// double* parameters
576// )
577//{
578// TF1* polExpFit = new TF1(name, NegPolExp, xMin, xMax, 3 );
579// polExpFit->SetParNames("Pol0", "Slope", "Shift");
580// polExpFit->SetLineColor(kRed);
581// histo->Fit(polExpFit, "+RWWM");
582// polExpFit->GetParameters(parameters);
583//}
584
585//double
586//NegPolExp(
587// double* x,
588// double* par
589// )
590//{
591// return par[0]+(-1)*TMath::Exp(par[1]+par[2]*x[0]);
592//}
593
594//double
595//PolExp(
596// double* x,
597// double* par
598// )
599//{
600// return
601//// par[0]+
602// TMath::Exp(par[1]+par[2]*x[0]);
603//}
604
605//double
606//ChargeDiode(
607// double time,
608// double chargeVoltage,
609// double impedance,
610// double capacity
611// )
612//{
613// return chargeVoltage*(1 - TMath::Exp(time/(impedance*capacity)));
614//}
615
616//double
617//UnChargeDiode(
618// double* time,
619// double* chargeVoltage,
620// double* timeConstant
621// )
622//{
623// return chargeVoltage[0]+TMath::Exp(chargeVoltage[1]+timeConstant[2]*time[0]);
624//// return chargeVoltage[0] * (TMath::Exp(time[0]/timeConstant[0]));
625//}
626
627//double
628//template_function(
629// double* input_x,
630// double* par)
631//{
632// double returnval = 0.0;
633
634// // I introduce a few names
635// // double shift = par[0];
636// double bsl = par[0];
637// double beginOfRisingEdge = par[1];
638// double p0RisingEdge = par[6];
639// double p1RisingEdge = par[7];
640// double p2RisingEdge = par[8];
641// double p3RisingEdge = par[9];
642// double endOfRisingEdge = par[2];
643//// double pOFallingEdge = par[3];
644//// double expPar1FallingEdge = par[4];
645//// double expPar1FallingEdge = par[5];
646// /*
647// bool couted_once = false;
648// if not couted_once
649// {
650// couted_once = true;
651// cout << "shift:" << shift << endl;
652// cout << "bsl:" << bsl << endl;
653// cout << "expars:" << endl;
654// cout << "\t factor:" << exppar[0] << endl;
655// cout << "\t tau:" << exppar[1] << endl;
656// cout << "\t t0:" << exppar[2] << endl;
657// cout << "pol3pars:" << endl;
658// cout << "p[0] + x p[1] + x^2 p[2] + x^3 p[3]" << endl;
659// cout << pol3par[0] << "\t" << pol3par[1] << "\t" << pol3par[2] << "\t" << pol3par[3] << endl;
660// cout << "ranges:" << endl;
661// cout << "begin of pol3: " << range[0] << endl;
662// cout << "begin of exp: " << range[1] << endl;
663// }
664// */
665// double x = input_x[0];
666
667// // the baseline is added everywhere.
668// returnval += bsl;
669
670// if ( (x > beginOfRisingEdge) && (x <= endOfRisingEdge) )
671// {
672// // from this point on the pol3 is added
673// returnval += p0RisingEdge;
674// returnval += p1RisingEdge * x;
675// returnval += p2RisingEdge * x*x;
676// returnval += p3RisingEdge * x*x*x;
677// }
678// else if ( x > endOfRisingEdge )
679// {
680// // from this point on the exp-func is added
681//// returnval += exppar[0] * TMath::Exp( exppar[1] * ( x - exppar[2] ) );
682// returnval += PolExp(input_x, par+3);
683// }
684
685// return returnval;
686//}
Note: See TracBrowser for help on using the repository browser.