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

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