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

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