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

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