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

Last change on this file since 14755 was 14751, checked in by Jens Buss, 12 years ago
add shift in X functions, improved shifting
  • Property svn:executable set to *
File size: 28.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// TCanvas test_canvas;
195// test_canvas.Divide(2,1);
196// test_canvas.cd(1);
197// inputHisto->Draw();
198 for (int i = 0; i < numIt; i++)
199 {
200// test_canvas.cd(2);
201 TH1* tempHisto = (TH1*)inputHisto->Clone("tempHisto");
202// "",
203// "",
204// inputHisto->GetNbinsX(),
205// inputHisto->GetXaxis()->GetBinLowEdge( inputHisto->GetXaxis()->GetFirst() ),
206// inputHisto->GetXaxis()->GetBinUpEdge( inputHisto->GetXaxis()->GetLast() )
207// );
208 tempHisto->Reset();
209 sample.BootstrapTH1(inputHisto, tempHisto);
210
211 Mean[i] = tempHisto->GetMean();
212 Median[i] = MedianOfH1 ( tempHisto );
213
214 Max[i] = tempHisto->GetBinCenter( tempHisto->GetMaximumBin() );
215
216 //improved determination of maximum
217// TF1 gaus("fgaus", "gaus", Max[i]-10, Max[i]+10);
218// tempHisto->Fit("fgaus", "WWQRN0");
219// Max[i] = gaus.GetParameter(1);
220
221
222// sample.SetSeed(sample.GetSeed() + 1);
223
224 /*
225 //Process gui events asynchronously during input
226 tempHisto->Draw();
227 test_canvas.Modified();
228 test_canvas.Update();
229 cout << endl;
230 TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
231 timer.TurnOn();
232 TString input = Getline("Type 'q' to exit, <return> to go on: ");
233 timer.TurnOff();
234 if (input=="q\n")
235 {
236 return;
237 }
238 */
239
240
241
242 delete tempHisto;
243 }
244
245 parameter[0] = TMath::Mean( numIt, Max);
246 parameterErr[0] = TMath::RMS( numIt, Max);
247 parameter[1] = TMath::Mean( numIt, Median);
248 parameterErr[1] = TMath::RMS( numIt, Median);
249 parameter[2] = TMath::Mean( numIt, Mean);
250 parameterErr[2] = TMath::RMS( numIt, Mean);
251
252 for (int i = 0; i < 3; i++)
253 {
254 if (&par[i] != NULL)
255 {
256 par[i] = parameter[i];
257 }
258 else cout << "par[" << i << "] array to small" << endl;
259 if (&par[i] != NULL)
260 {
261 parErr[i] = parameterErr[i];
262 }
263 else cout << "parErr[" << i << "] array to small" << endl;
264 }
265
266 return;
267}
268
269void
270PlotMedianOfSlice(
271 Pixel* CurrentPixel,
272 TString overlayType,
273 int verbosityLevel
274 )
275{
276 if (verbosityLevel > 2) cout << endl
277 << "...calculating pulse shape of slice's Median"
278 << endl;
279
280 TH2F** hInputHistoArray = NULL;
281 TH1F** hOutputHistoArray = NULL;
282// TH1* hTempHisto = NULL;
283// float median = 0;
284
285 if (overlayType.Contains("Edge"))
286 {
287 hInputHistoArray = CurrentPixel->hEdgeOverlay;
288 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
289
290 }
291 else if (overlayType.Contains("Maximum"))
292 {
293 hInputHistoArray = CurrentPixel->hMaxOverlay;
294 hOutputHistoArray = CurrentPixel->hPixelMedian;
295 }
296 else
297 {
298 cout << endl << "Unknown Overlay Method-->aborting" << endl;
299 return;
300 }
301
302 for (int pulse_order = 0 ;
303 pulse_order < CurrentPixel->mMaxPulseOrder ;
304 pulse_order ++)
305 {
306 if (verbosityLevel > 2) cout << "\t...calculation of "
307 << "pulse order " << pulse_order;
308
309 MedianOfTH2Slices(
310 hInputHistoArray[pulse_order],
311 hOutputHistoArray[pulse_order],
312 verbosityLevel
313 );
314
315 ErrorMedianOfTH2Slices(
316 hInputHistoArray[pulse_order],
317 hOutputHistoArray[pulse_order],
318 5, //numIterations
319 verbosityLevel
320 );
321
322// Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
323
324// for (Int_t slice=1;slice<=nbins;slice++) {
325
326// hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
327// median = MedianOfH1(hTempHisto);
328
329// if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
330
331// hOutputHistoArray[pulse_order]->SetBinContent(slice, median );
332//// delete h1;
333// }
334
335 if (verbosityLevel > 2) cout << "\t...done" << endl;
336 }
337}
338// end of PlotMedianEachSliceOfPulse
339//----------------------------------------------------------------------------
340
341void
342MedianOfTH2Slices(
343 TH2* inputHisto,
344 TH1* outputHisto,
345 int verbosityLevel
346 )
347{
348 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
349 float median = 0;
350
351 for (Int_t slice=1;slice<=nbins;slice++) {
352
353 median = MedianOfH1( inputHisto->ProjectionY("",slice,slice) );
354
355 if (verbosityLevel > 2) printf("Median of Slice %d, Median=%g\n",slice,median);
356
357 outputHisto->SetBinContent(slice, median );
358// delete h1;
359 }
360 return;
361}
362// end of MedianOfTH2Slices
363//----------------------------------------------------------------------------
364
365double
366ErrorMedianOfH1(
367 TH1* inputHisto,
368 int numIterations,
369 int verbosityLevel
370 )
371{
372 if (verbosityLevel > 5) cout << "...calculating errors of median" << endl;
373// float MedianOfSliceMean;
374 double medianRMS = 0;
375
376 int sample_min = inputHisto->GetXaxis()->GetFirst();
377 int sample_max = inputHisto->GetXaxis()->GetLast();
378 int sample_size = inputHisto->GetXaxis()->GetNbins();
379// int sample_size = sample_max - sample_min;
380// cout << "sample_min " << sample_min
381// << ", sample_max " << sample_max
382// << ", sample_size " << sample_size
383// << endl;
384 double median[numIterations];
385 vector<int> rndList;
386 Sample sample(sample_size);
387
388 for (int i = 0; i < numIterations; i++)
389 {
390 sample.SetSeed(sample.GetSeed() + 1);
391 sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
392 median[i] = MedianOfH1withRndSlices(
393 inputHisto,
394 &rndList
395 );
396 }
397
398// MedianOfSliceMean = TMath::Mean(numIterations, median);
399 medianRMS = RMS(numIterations, median);
400// if (verbosityLevel > 2) printf("Mean Median=%g\n",MedianOfSliceMean);
401 if (verbosityLevel > 5) printf("medianRMS=%g\n",medianRMS);
402
403 return medianRMS;
404}
405
406void
407ErrorMedianOfTH2Slices(
408 TH2* inputHisto,
409 TH1* outputHisto,
410 int numIterations,
411 int verbosityLevel
412 )
413{
414 if (verbosityLevel > 2) cout << "...calculating errors of median" << endl;
415 Int_t nbins = inputHisto->GetXaxis()->GetNbins();
416 cout << "nbins " << nbins << endl;
417 float MedianOfSliceMean[nbins];
418 float MedianOfSliceRMS[nbins];
419 float median[numIterations];
420 vector<int> rndList;
421// rndList->clear();
422
423 TH1 *htemp = NULL;
424
425 for (Int_t slice=1;slice<=nbins;slice++) {
426 htemp = inputHisto->ProjectionY("",slice,slice);
427
428
429 int sample_min = htemp->GetXaxis()->GetFirst();
430 int sample_max = htemp->GetXaxis()->GetLast();
431 int sample_size = htemp->GetXaxis()->GetNbins();
432
433// cout << "sample_min " << sample_min
434// << ", sample_max " << sample_max
435// << ", sample_size " << sample_size
436// << endl;
437
438 Sample sample(sample_size);
439 for (int i = 0; i < numIterations; i++)
440 {
441 sample.SetSeed(sample.GetSeed() + 1);
442 sample.BootstrapSample(&rndList, sample_min, sample_max, sample_size);
443 median[i] = MedianOfH1withRndSlices(
444 htemp,
445 &rndList
446 );
447
448 }
449 MedianOfSliceMean[slice] = TMath::Mean(numIterations, median);
450 MedianOfSliceRMS[slice] = TMath::RMS(numIterations, median);
451 outputHisto->SetBinError(slice, MedianOfSliceRMS[slice]);
452// outputHisto->SetBinError(slice, RMS(numIterations, median) );
453 if (verbosityLevel > 2) printf("Mean Median of Slice %d, Mean Median=%g\n",slice,MedianOfSliceMean[slice]);
454 if (verbosityLevel > 2) printf("RMS of Median of Slice %d, MedianRMS=%g\n",slice,MedianOfSliceRMS[slice]);
455// outputHisto->SetBinContent(slice, median );
456// delete h1;
457 }
458 return;
459}
460// end of ErrorMedianOfTH2Slices
461//----------------------------------------------------------------------------
462
463void
464PlotMeanOfSlice(
465 Pixel* CurrentPixel,
466 TString overlayType,
467 int verbosityLevel
468 )
469{
470 if (verbosityLevel > 2) cout << endl
471 << "...calculating pulse shape of slice's Mean"
472 << endl;
473
474 TH2F** hInputHistoArray = NULL;
475 TH1F** hOutputHistoArray = NULL;
476 TH1* hTempHisto = NULL;
477 float mean = 0;
478
479 if (overlayType.Contains("Edge"))
480 {
481 hInputHistoArray = CurrentPixel->hEdgeOverlay;
482 hOutputHistoArray = CurrentPixel->hPixelEdgeMean;
483
484 }
485 else if (overlayType.Contains("Maximum"))
486 {
487 hInputHistoArray = CurrentPixel->hMaxOverlay;
488 hOutputHistoArray = CurrentPixel->hPixelMean;
489 }
490 else
491 {
492 cout << endl << "Unknown Overlay Method-->aborting" << endl;
493 return;
494 }
495
496 for (int pulse_order = 0 ;
497 pulse_order < CurrentPixel->mMaxPulseOrder ;
498 pulse_order ++)
499 {
500 if (verbosityLevel > 2) cout << "\t...calculation of "
501 << "pulse order " << pulse_order;
502
503 Int_t nbins = hInputHistoArray[pulse_order]->GetXaxis()->GetNbins();
504
505 for (Int_t slice=1;slice<=nbins;slice++) {
506
507 hTempHisto = hInputHistoArray[pulse_order]->ProjectionY("",slice,slice);
508 mean = hTempHisto->GetMean();
509
510 if (verbosityLevel > 2) printf("Mean of Slice %d, Mean=%g\n",slice,mean);
511
512 hOutputHistoArray[pulse_order]->SetBinContent(slice, mean );
513// delete h1;
514 }
515
516 if (verbosityLevel > 2) cout << "\t...done" << endl;
517 }
518}
519// end of CalcMeanOfSlice
520//----------------------------------------------------------------------------
521
522void
523ExtractPulseTemplate(
524 Pixel* CurrentPixel,
525 TString overlayType,
526 int pulse_order,
527 int verbosityLevel
528 )
529{
530 if (verbosityLevel > 2) cout << endl
531 << "...calculating pulse shape"
532 << " of max. propability"
533 << " of "
534 << overlayType
535 << " Overlay"
536 << endl;
537 TH2F* hInputHisto = NULL;
538 TH1F* hOutputMaxHisto = NULL;
539 TH1F* hOutputMeanHisto = NULL;
540 TH1F* hOutputMedianHisto = NULL;
541 TH1* hTempHisto = NULL;
542 float max_prop = 0;
543 float median = 0;
544 float mean = 0;
545 float max_prop_err = 0;
546 float median_err = 0;
547 float mean_err = 0;
548
549 if (verbosityLevel > 3)
550 {
551 cout << "...setting pointers to histogramm arrays ";
552 cout << " for " << overlayType << "Overlay" << endl;
553 }
554 if (overlayType.Contains("Edge"))
555 {
556 hInputHisto = (CurrentPixel->hEdgeOverlay[pulse_order]);
557
558 //Maximum Propability of Slices
559 hOutputMaxHisto = (CurrentPixel->hPixelEdgeMax[pulse_order]);
560
561 //Mean of Slices
562 hOutputMeanHisto = (CurrentPixel->hPixelEdgeMean[pulse_order]);
563
564 //Median of Slices
565 hOutputMedianHisto = (CurrentPixel->hPixelEdgeMedian[pulse_order]);
566 }
567 else if (overlayType.Contains("Maximum"))
568 {
569 hInputHisto = (CurrentPixel->hMaxOverlay[pulse_order]);
570
571 //Maximum Propability of Slices
572 hOutputMaxHisto = (CurrentPixel->hPixelMax[pulse_order]);
573
574 //Mean of Slices
575 hOutputMeanHisto = (CurrentPixel->hPixelMean[pulse_order]);
576
577 //Median of Slices
578 hOutputMedianHisto = (CurrentPixel->hPixelMedian[pulse_order]);
579 }
580 else
581 {
582 cout << endl << overlayType << "Unknown Overlay Method-->aborting" << endl;
583 return;
584 }
585 if (verbosityLevel > 3)
586 {
587 cout << "...done " << endl;
588 }
589
590 if (verbosityLevel > 2)
591 {
592 cout << "\t...calculation of "
593 << "pulse order " << pulse_order << endl;
594 }
595
596// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
597
598 int slice_min = hInputHisto->GetXaxis()->GetFirst();
599 int slice_max = hInputHisto->GetXaxis()->GetLast();
600
601 for (Int_t slice=slice_min;slice<=slice_max;slice++)
602 {
603
604 hTempHisto = hInputHisto->ProjectionY("",slice,slice);
605
606 //containers for errors
607 double slMean[3];
608 double slError[3];
609
610 //calculate Errors with bootstrapping
611 CalculateErrorsWithBootstrapping(hTempHisto, 10, slMean, slError);
612
613 if (verbosityLevel > 3)
614 {
615 cout << "\t\t...calculating maxProb of slice " << slice << endl;
616 }
617
618 //Get maximum of slice's distribution
619 max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
620// max_prop = slMean[0];
621
622 //improve result by< fitting gaussian to slices distribution
623 TF1 gaus("fgaus", "gaus", max_prop-30, max_prop+30);
624 hTempHisto->Fit("fgaus", "QRN0");
625 max_prop = gaus.GetParameter(1);
626
627 //calculate error of max prop
628// max_prop_err = gaus.GetParameter(2); //RMS of gaus fit of slices distribution
629 max_prop_err = slError[0]; //error calculated with bootstrapping
630// max_prop_err = hTempHisto->GetRMS(); //RMS of slice
631
632 if (verbosityLevel > 3)
633 {
634 cout << "\t\t...calculating Median of slice " << slice << endl;
635 }
636 median = MedianOfH1(hTempHisto);
637// median = slMean[1];
638 median_err = slError[1]; //error from bootstraping
639
640 if (verbosityLevel > 3) cout << "\t\t...calculating Mean of slice " << slice << endl;
641
642
643 mean = hTempHisto->GetMean();
644// mean = slMean[2];
645 mean_err = hTempHisto->GetRMS(); //RMS of slice
646// mean_err = slError[2]; //error from bootstraping
647
648 if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
649 hOutputMaxHisto->SetBinContent(slice, max_prop );
650 hOutputMaxHisto->SetBinError( slice, max_prop_err);
651
652 if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
653 hOutputMeanHisto->SetBinContent(slice, mean );
654 hOutputMeanHisto->SetBinError( slice, mean_err);
655
656 if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
657 hOutputMedianHisto->SetBinContent( slice, median );
658 hOutputMedianHisto->SetBinError( slice, median_err);
659 delete hTempHisto;
660
661 }//Loop over slices
662
663// ErrorMedianOfTH2Slices(
664// hInputHisto,
665// hOutputMedianHisto,
666// 100, //numIterations
667// verbosityLevel
668// );
669// hOutputMaxHisto->GetXaxis()->SetLimits(
670// 0,
671// 300
672// );
673
674// hOutputMeanHisto->GetXaxis()->SetLimits(
675// 0,
676// 300
677// );
678
679// hOutputMedianHisto->GetXaxis()->SetLimits(
680// 0,
681// 300
682// );
683
684 ShiftStartOfHistoInXTo(
685 hOutputMaxHisto,
686 0
687 );
688
689 ShiftStartOfHistoInXTo(
690 hOutputMeanHisto,
691 0
692 );
693
694 ShiftStartOfHistoInXTo(
695 hOutputMedianHisto,
696 0
697 );
698
699}
700// end of PlotMaxPropabilityPulse
701//----------------------------------------------------------------------------
702
703bool
704WritePixelTemplateToCsv(
705 Pixel* CurrentPixel,
706 TString path,
707 TString overlayMethod,
708 int order,
709 int verbosityLevel
710 )
711{
712// TSystem this_system();
713// this_system.cd(path);
714// cout << this_system.pwd();
715// this_system.mkdir("CSV");
716 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
717
718 TH1F* Max_histo = NULL;
719 TH1F* Median_histo = NULL;
720 TH1F* Mean_histo = NULL;
721
722 if (overlayMethod.Contains("Maximum"))
723 {
724 Max_histo = CurrentPixel->hPixelMax[order];
725 Median_histo = CurrentPixel->hPixelMedian[order];
726 Mean_histo = CurrentPixel->hPixelMean[order];
727 }
728 else if (overlayMethod.Contains("Edge"))
729 {
730 Max_histo = CurrentPixel->hPixelMax[order];
731 Median_histo = CurrentPixel->hPixelMedian[order];
732 Mean_histo = CurrentPixel->hPixelMean[order];
733 }
734 else
735 {
736 cout << endl << "Unknown Overlay Method-->aborting" << endl;
737 return 1;
738 }
739
740
741 Int_t nbins = Max_histo->GetXaxis()->GetNbins();
742
743 if (verbosityLevel > 2)
744 {
745 cout << "writing point-set to csv file: " ;
746 cout << path << endl;
747 cout << "...opening file" << endl;
748 }
749 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
750 ofstream out;
751 out.open( path );
752 out << "### point-set of a single photon pulse template" << endl;
753 out << "### template determined with pulse overlay at: "
754 << overlayMethod << endl;
755 out << "### Slice's Amplitude determined by calculating the " << endl
756 << "### value of maximum propability of slice -> AmplitudeMax " << endl
757 << "### mean of slice -> AmplitudeMean " << endl
758 << "### median of slice -> AmplitudeMedian " << endl
759 << "### for each slice" << endl;
760 out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
761 << endl;
762
763 out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
764
765 for (int slice=1;slice<=nbins;slice++)
766 {
767 out << slice << "," ;
768 out << Max_histo->GetBinContent(slice) << ",";
769 out << Mean_histo->GetBinContent(slice) << ",";
770 out << Median_histo->GetBinContent(slice) << endl;
771 }
772 out.close();
773 if (verbosityLevel > 2) cout << "...csv file closed" << endl;
774 return 0;
775}
776
777void
778FitMaxPropabilityPuls(
779 TH1F* hMaximumTemp,
780 int verbosityLevel
781 )
782 {
783 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
784 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
785 hMaximumTemp->Fit("landau", "", "", -50, 250);
786 if (verbosityLevel > 2) cout << "...done" << endl;
787 }
788
789//void
790//FitFallingEdge(
791// TString name,
792// TH1F* histo,
793// double xMin,
794// double xMax,
795// double* parameters
796// )
797//{
798// TF1* polExpFit = new TF1(name, PolExp, xMin, xMax, 3 );
799// polExpFit->SetParNames("Pol0", "Slope", "Shift");
800// polExpFit->SetLineColor(kGreen);
801// histo->Fit(polExpFit, "+RWM");
802// polExpFit->GetParameters(parameters);
803//}
804
805//void
806//FitRisingEdge(
807// TString name,
808// TH1F* histo,
809// double xMin,
810// double xMax,
811// double* parameters
812// )
813//{
814// TF1* polExpFit = new TF1(name, NegPolExp, xMin, xMax, 3 );
815// polExpFit->SetParNames("Pol0", "Slope", "Shift");
816// polExpFit->SetLineColor(kRed);
817// histo->Fit(polExpFit, "+RWWM");
818// polExpFit->GetParameters(parameters);
819//}
820
821//double
822//NegPolExp(
823// double* x,
824// double* par
825// )
826//{
827// return par[0]+(-1)*TMath::Exp(par[1]+par[2]*x[0]);
828//}
829
830//double
831//PolExp(
832// double* x,
833// double* par
834// )
835//{
836// return
837//// par[0]+
838// TMath::Exp(par[1]+par[2]*x[0]);
839//}
840
841//double
842//ChargeDiode(
843// double time,
844// double chargeVoltage,
845// double impedance,
846// double capacity
847// )
848//{
849// return chargeVoltage*(1 - TMath::Exp(time/(impedance*capacity)));
850//}
851
852//double
853//UnChargeDiode(
854// double* time,
855// double* chargeVoltage,
856// double* timeConstant
857// )
858//{
859// return chargeVoltage[0]+TMath::Exp(chargeVoltage[1]+timeConstant[2]*time[0]);
860//// return chargeVoltage[0] * (TMath::Exp(time[0]/timeConstant[0]));
861//}
862
863//double
864//template_function(
865// double* input_x,
866// double* par)
867//{
868// double returnval = 0.0;
869
870// // I introduce a few names
871// // double shift = par[0];
872// double bsl = par[0];
873// double beginOfRisingEdge = par[1];
874// double p0RisingEdge = par[6];
875// double p1RisingEdge = par[7];
876// double p2RisingEdge = par[8];
877// double p3RisingEdge = par[9];
878// double endOfRisingEdge = par[2];
879//// double pOFallingEdge = par[3];
880//// double expPar1FallingEdge = par[4];
881//// double expPar1FallingEdge = par[5];
882// /*
883// bool couted_once = false;
884// if not couted_once
885// {
886// couted_once = true;
887// cout << "shift:" << shift << endl;
888// cout << "bsl:" << bsl << endl;
889// cout << "expars:" << endl;
890// cout << "\t factor:" << exppar[0] << endl;
891// cout << "\t tau:" << exppar[1] << endl;
892// cout << "\t t0:" << exppar[2] << endl;
893// cout << "pol3pars:" << endl;
894// cout << "p[0] + x p[1] + x^2 p[2] + x^3 p[3]" << endl;
895// cout << pol3par[0] << "\t" << pol3par[1] << "\t" << pol3par[2] << "\t" << pol3par[3] << endl;
896// cout << "ranges:" << endl;
897// cout << "begin of pol3: " << range[0] << endl;
898// cout << "begin of exp: " << range[1] << endl;
899// }
900// */
901// double x = input_x[0];
902
903// // the baseline is added everywhere.
904// returnval += bsl;
905
906// if ( (x > beginOfRisingEdge) && (x <= endOfRisingEdge) )
907// {
908// // from this point on the pol3 is added
909// returnval += p0RisingEdge;
910// returnval += p1RisingEdge * x;
911// returnval += p2RisingEdge * x*x;
912// returnval += p3RisingEdge * x*x*x;
913// }
914// else if ( x > endOfRisingEdge )
915// {
916// // from this point on the exp-func is added
917//// returnval += exppar[0] * TMath::Exp( exppar[1] * ( x - exppar[2] ) );
918// returnval += PolExp(input_x, par+3);
919// }
920
921// return returnval;
922//}
923
924
925void
926ShiftHistoInX(
927 TH1* histo,
928 float shift
929 )
930{
931 histo->GetXaxis()->SetLimits(
932 histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() ) + shift,
933 histo->GetXaxis()->GetBinUpEdge( histo->GetXaxis()->GetLast() ) + shift
934 );
935}
936
937void
938ShiftStartOfHistoInXTo(
939 TH1* histo,
940 float value
941 )
942{
943 float min = histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() );
944 float shift = value - min;
945
946 ShiftHistoInX( histo, shift);
947}
Note: See TracBrowser for help on using the repository browser.