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

Last change on this file since 14946 was 14946, checked in by Jens Buss, 12 years ago
added dynamic startvalue calculation
  • Property svn:executable set to *
File size: 28.9 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 numBootstrapIt,// = 10, //number of iterations for bootstrapping
528 int verbosityLevel
529 )
530{
531 if (verbosityLevel > 2) cout << endl
532 << "...calculating pulse shape"
533 << " of max. propability"
534 << " of "
535 << overlayType
536 << " Overlay"
537 << endl;
538 TH2F* hInputHisto = NULL;
539 TH1F* hOutputMaxHisto = NULL;
540 TH1F* hOutputMeanHisto = NULL;
541 TH1F* hOutputMedianHisto = NULL;
542 TH1D* hTempHisto = NULL;
543 double max_prop = 0;
544 double median = 0;
545 double mean = 0;
546 double max_prop_err = 0;
547 double median_err = 0;
548 double mean_err = 0;
549
550 if (verbosityLevel > 3)
551 {
552 cout << "...setting pointers to histogramm arrays ";
553 cout << " for " << overlayType << "Overlay" << endl;
554 }
555 if (overlayType.Contains("Edge"))
556 {
557 hInputHisto = (CurrentPixel->hEdgeOverlay[pulse_order]);
558
559 //Maximum Propability of Slices
560 hOutputMaxHisto = (CurrentPixel->hPixelEdgeMax[pulse_order]);
561
562 //Mean of Slices
563 hOutputMeanHisto = (CurrentPixel->hPixelEdgeMean[pulse_order]);
564
565 //Median of Slices
566 hOutputMedianHisto = (CurrentPixel->hPixelEdgeMedian[pulse_order]);
567 }
568 else if (overlayType.Contains("Maximum"))
569 {
570 hInputHisto = (CurrentPixel->hMaxOverlay[pulse_order]);
571
572 //Maximum Propability of Slices
573 hOutputMaxHisto = (CurrentPixel->hPixelMax[pulse_order]);
574
575 //Mean of Slices
576 hOutputMeanHisto = (CurrentPixel->hPixelMean[pulse_order]);
577
578 //Median of Slices
579 hOutputMedianHisto = (CurrentPixel->hPixelMedian[pulse_order]);
580 }
581 else
582 {
583 cout << endl << overlayType << "Unknown Overlay Method-->aborting" << endl;
584 return;
585 }
586 if (verbosityLevel > 3)
587 {
588 cout << "...done " << endl;
589 }
590
591 if (verbosityLevel > 2)
592 {
593 cout << "\t...calculation of "
594 << "pulse order " << pulse_order << endl;
595 }
596
597// if (verbosityLevel > 2) cout << "\t...# slices processed " << nbins << endl;
598
599 int slice_min = hInputHisto->GetXaxis()->GetFirst();
600 int slice_max = hInputHisto->GetXaxis()->GetLast();
601
602 for (Int_t slice=slice_min;slice<=slice_max;slice++)
603 {
604
605 hTempHisto = hInputHisto->ProjectionY("",slice,slice);
606
607 //containers for errors
608 double slMean[3];
609 double slError[3];
610
611 //calculate Errors with bootstrapping
612 CalculateErrorsWithBootstrapping(hTempHisto, numBootstrapIt, slMean, slError);
613
614 if (verbosityLevel > 3)
615 {
616 cout << "\t\t...calculating maxProb of slice " << slice << endl;
617 }
618
619 //Get maximum of slice's distribution
620// max_prop = hTempHisto->GetBinCenter( hTempHisto->GetMaximumBin() );
621 max_prop = slMean[0];
622
623 //improve result by< fitting gaussian to slices distribution
624// TF1 gaus("fgaus", "gaus", max_prop-30, max_prop+30);
625// hTempHisto->Fit("fgaus", "QRN0");
626// max_prop = gaus.GetParameter(1);
627
628 //calculate error of max prop
629// max_prop_err = gaus.GetParameter(2); //RMS of gaus fit of slices distribution
630 max_prop_err = slError[0]; //error calculated with bootstrapping
631// max_prop_err = hTempHisto->GetRMS(); //RMS of slice
632
633 if (verbosityLevel > 3)
634 {
635 cout << "\t\t...calculating Median of slice " << slice << endl;
636 }
637// median = MedianOfH1(hTempHisto);
638 median = slMean[1];
639 median_err = slError[1]; //error from bootstraping
640
641 if (verbosityLevel > 3) cout << "\t\t...calculating Mean of slice " << slice << endl;
642
643
644// mean = hTempHisto->GetMean();
645 mean = slMean[2];
646// mean_err = hTempHisto->GetRMS()/hTempHisto->GetEntries(); //RMS of slice
647 mean_err = slError[2]; //error from bootstraping
648
649 if (verbosityLevel > 4) cout << "\t\t\t\t MaxProb of Slice " << slice << ": " << max_prop << endl;
650 hOutputMaxHisto->SetBinContent(slice, max_prop );
651 hOutputMaxHisto->SetBinError( slice, max_prop_err);
652
653 if (verbosityLevel > 4) cout << "\t\t\t\t Mean of Slice " << slice << ": " << mean << endl;
654 hOutputMeanHisto->SetBinContent(slice, mean );
655 hOutputMeanHisto->SetBinError( slice, mean_err);
656
657 if (verbosityLevel > 4) cout << "\t\t\t\t Median of Slice " << slice << ": " << median << endl;
658 hOutputMedianHisto->SetBinContent( slice, median );
659 hOutputMedianHisto->SetBinError( slice, median_err);
660 delete hTempHisto;
661
662 }//Loop over slices
663
664// ErrorMedianOfTH2Slices(
665// hInputHisto,
666// hOutputMedianHisto,
667// 100, //numIterations
668// verbosityLevel
669// );
670// hOutputMaxHisto->GetXaxis()->SetLimits(
671// 0,
672// 300
673// );
674
675// hOutputMeanHisto->GetXaxis()->SetLimits(
676// 0,
677// 300
678// );
679
680// hOutputMedianHisto->GetXaxis()->SetLimits(
681// 0,
682// 300
683// );
684
685 ShiftStartOfHistoInXTo(
686 hOutputMaxHisto,
687 0
688 );
689
690 ShiftStartOfHistoInXTo(
691 hOutputMeanHisto,
692 0
693 );
694
695 ShiftStartOfHistoInXTo(
696 hOutputMedianHisto,
697 0
698 );
699
700 ShiftStartOfHistoInXTo(
701 hInputHisto,
702 0
703 );
704
705}
706// end of PlotMaxPropabilityPulse
707//----------------------------------------------------------------------------
708
709bool
710WritePixelTemplateToCsv(
711 Pixel* CurrentPixel,
712 TString path,
713 TString overlayMethod,
714 int order,
715 int verbosityLevel
716 )
717{
718// TSystem this_system();
719// this_system.cd(path);
720// cout << this_system.pwd();
721// this_system.mkdir("CSV");
722 path = CurrentPixel->CsvFileName( path, overlayMethod, order);
723
724 TH1F* Max_histo = NULL;
725 TH1F* Median_histo = NULL;
726 TH1F* Mean_histo = NULL;
727
728 if (overlayMethod.Contains("Maximum"))
729 {
730 Max_histo = CurrentPixel->hPixelMax[order];
731 Median_histo = CurrentPixel->hPixelMedian[order];
732 Mean_histo = CurrentPixel->hPixelMean[order];
733 }
734 else if (overlayMethod.Contains("Edge"))
735 {
736 Max_histo = CurrentPixel->hPixelMax[order];
737 Median_histo = CurrentPixel->hPixelMedian[order];
738 Mean_histo = CurrentPixel->hPixelMean[order];
739 }
740 else
741 {
742 cout << endl << "Unknown Overlay Method-->aborting" << endl;
743 return 1;
744 }
745
746
747 Int_t nbins = Max_histo->GetXaxis()->GetNbins();
748
749 if (verbosityLevel > 2)
750 {
751 cout << "writing point-set to csv file: " ;
752 cout << path << endl;
753 cout << "...opening file" << endl;
754 }
755 if (verbosityLevel > 2) cout << "...number of bins " << nbins << endl;
756 ofstream out;
757 out.open( path );
758 out << "### point-set of a single photon pulse template" << endl;
759 out << "### template determined with pulse overlay at: "
760 << overlayMethod << endl;
761 out << "### Slice's Amplitude determined by calculating the " << endl
762 << "### value of maximum propability of slice -> AmplitudeMax " << endl
763 << "### mean of slice -> AmplitudeMean " << endl
764 << "### median of slice -> AmplitudeMedian " << endl
765 << "### for each slice" << endl;
766 out << "### Pixel number (CHid): " << CurrentPixel->mChid << endl
767 << endl;
768
769 out << "time [slices],AmplitudeMax [mV],AmplitudeMean [mV],AmplitudeMedian [mV]" << endl;
770
771 for (int slice=1;slice<=nbins;slice++)
772 {
773 out << slice << "," ;
774 out << Max_histo->GetBinContent(slice) << ",";
775 out << Mean_histo->GetBinContent(slice) << ",";
776 out << Median_histo->GetBinContent(slice) << endl;
777 }
778 out.close();
779 if (verbosityLevel > 2) cout << "...csv file closed" << endl;
780 return 0;
781}
782
783void
784FitMaxPropabilityPuls(
785 TH1F* hMaximumTemp,
786 int verbosityLevel
787 )
788 {
789 if (verbosityLevel > 2) cout << "...fit Landau in histograms" ;
790 if (verbosityLevel > 2) cout << "\t...calculating Landaufit" ;
791 hMaximumTemp->Fit("landau", "", "", -50, 250);
792 if (verbosityLevel > 2) cout << "...done" << endl;
793 }
794
795//void
796//FitFallingEdge(
797// TString name,
798// TH1F* histo,
799// double xMin,
800// double xMax,
801// double* parameters
802// )
803//{
804// TF1* polExpFit = new TF1(name, PolExp, xMin, xMax, 3 );
805// polExpFit->SetParNames("Pol0", "Slope", "Shift");
806// polExpFit->SetLineColor(kGreen);
807// histo->Fit(polExpFit, "+RWM");
808// polExpFit->GetParameters(parameters);
809//}
810
811//void
812//FitRisingEdge(
813// TString name,
814// TH1F* histo,
815// double xMin,
816// double xMax,
817// double* parameters
818// )
819//{
820// TF1* polExpFit = new TF1(name, NegPolExp, xMin, xMax, 3 );
821// polExpFit->SetParNames("Pol0", "Slope", "Shift");
822// polExpFit->SetLineColor(kRed);
823// histo->Fit(polExpFit, "+RWWM");
824// polExpFit->GetParameters(parameters);
825//}
826
827//double
828//NegPolExp(
829// double* x,
830// double* par
831// )
832//{
833// return par[0]+(-1)*TMath::Exp(par[1]+par[2]*x[0]);
834//}
835
836//double
837//PolExp(
838// double* x,
839// double* par
840// )
841//{
842// return
843//// par[0]+
844// TMath::Exp(par[1]+par[2]*x[0]);
845//}
846
847//double
848//ChargeDiode(
849// double time,
850// double chargeVoltage,
851// double impedance,
852// double capacity
853// )
854//{
855// return chargeVoltage*(1 - TMath::Exp(time/(impedance*capacity)));
856//}
857
858//double
859//UnChargeDiode(
860// double* time,
861// double* chargeVoltage,
862// double* timeConstant
863// )
864//{
865// return chargeVoltage[0]+TMath::Exp(chargeVoltage[1]+timeConstant[2]*time[0]);
866//// return chargeVoltage[0] * (TMath::Exp(time[0]/timeConstant[0]));
867//}
868
869//double
870//template_function(
871// double* input_x,
872// double* par)
873//{
874// double returnval = 0.0;
875
876// // I introduce a few names
877// // double shift = par[0];
878// double bsl = par[0];
879// double beginOfRisingEdge = par[1];
880// double p0RisingEdge = par[6];
881// double p1RisingEdge = par[7];
882// double p2RisingEdge = par[8];
883// double p3RisingEdge = par[9];
884// double endOfRisingEdge = par[2];
885//// double pOFallingEdge = par[3];
886//// double expPar1FallingEdge = par[4];
887//// double expPar1FallingEdge = par[5];
888// /*
889// bool couted_once = false;
890// if not couted_once
891// {
892// couted_once = true;
893// cout << "shift:" << shift << endl;
894// cout << "bsl:" << bsl << endl;
895// cout << "expars:" << endl;
896// cout << "\t factor:" << exppar[0] << endl;
897// cout << "\t tau:" << exppar[1] << endl;
898// cout << "\t t0:" << exppar[2] << endl;
899// cout << "pol3pars:" << endl;
900// cout << "p[0] + x p[1] + x^2 p[2] + x^3 p[3]" << endl;
901// cout << pol3par[0] << "\t" << pol3par[1] << "\t" << pol3par[2] << "\t" << pol3par[3] << endl;
902// cout << "ranges:" << endl;
903// cout << "begin of pol3: " << range[0] << endl;
904// cout << "begin of exp: " << range[1] << endl;
905// }
906// */
907// double x = input_x[0];
908
909// // the baseline is added everywhere.
910// returnval += bsl;
911
912// if ( (x > beginOfRisingEdge) && (x <= endOfRisingEdge) )
913// {
914// // from this point on the pol3 is added
915// returnval += p0RisingEdge;
916// returnval += p1RisingEdge * x;
917// returnval += p2RisingEdge * x*x;
918// returnval += p3RisingEdge * x*x*x;
919// }
920// else if ( x > endOfRisingEdge )
921// {
922// // from this point on the exp-func is added
923//// returnval += exppar[0] * TMath::Exp( exppar[1] * ( x - exppar[2] ) );
924// returnval += PolExp(input_x, par+3);
925// }
926
927// return returnval;
928//}
929
930
931void
932ShiftHistoInX(
933 TH1* histo,
934 float shift
935 )
936{
937 histo->GetXaxis()->SetLimits(
938 histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() ) + shift,
939 histo->GetXaxis()->GetBinUpEdge( histo->GetXaxis()->GetLast() ) + shift
940 );
941}
942
943void
944ShiftStartOfHistoInXTo(
945 TH1* histo,
946 float value
947 )
948{
949 float min = histo->GetXaxis()->GetBinLowEdge( histo->GetXaxis()->GetFirst() );
950 float shift = value - min;
951
952 ShiftHistoInX( histo, shift);
953}
Note: See TracBrowser for help on using the repository browser.