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

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