GammaWare  Head Version for release 0.9
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Peak1D.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * Copyright (C) 2004 by Olivier Stezowski, B. Rosse & Christian Finck *
3  * stezow(AT)ipnl.in2p3.fr, cfinck(AT)ires.in2p3.fr *
4  * *
5  * This program is free software; you can redistribute it and/or modif *
6  * it under the terms of the GNU General Public License as published by *
7  * the Free Software Foundation; either version 2 of the License, or *
8  * (at your option) any later version. *
9  * *
10  * This program is distributed in the hope that it will be useful, *
11  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
13  * GNU General Public License for more details. *
14  * *
15  * You should have received a copy of the GNU General Public License *
16  * along with this program; if not, write to the *
17  * Free Software Foundation, Inc., *
18  * 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
19  ***************************************************************************/
20 
23 #ifndef ROOT_TF1
24 #include "TF1.h"
25 #endif
26 
27 #ifndef ROOT_TH1F
28 #include "TH1.h"
29 #endif
30 
31 #ifndef ROOT_TList
32 #include "TList.h"
33 #endif
34 
35 #ifndef ROOT_TMarker
36 #include "TMarker.h"
37 #endif
38 
39 #ifndef ROOT_TMath
40 #include "TMath.h"
41 #endif
42 
43 #ifndef ROOT_TPad
44 #include "TPad.h"
45 #endif
46 
47 #ifndef ROOT_TPolyMarker
48 #include "TMarker.h"
49 #endif
50 
51 #ifndef ROOT_TROOT
52 #include "TROOT.h"
53 #endif
54 
55 #include <RVersion.h>
56 
57 #ifndef Gw_PadManager
58 #include "PadManager.h"
59 #endif
60 
61 #ifndef Gw_Peak1D
62 #include "Peak1D.h"
63 #endif
64 
65 #include "PeakCreator.h"
66 
67 #include <iostream>
68 #include <iomanip>
69 
70 // to work with FWHM instead of sigma
71 #define SigmaToFWHM 2.35482
72 
73 using namespace Gw;
74 using namespace std;
75 
76 // COULD NOT BE CHANGED OTHERWISE BREAK ROOT STREAMING
77 const Int_t Peak1D::fgkPointsPeak = 7;
78 const Int_t Peak1D::fgkPointsBkg = 4;
79 const Int_t Peak1D::fgkPeakIdx = 3;
80 
81 // Could be modified
82 Style_t Peak1D::fgFillStylePeak = 3004;
83 Color_t Peak1D::fgLineColorBkg = 39;
84 Style_t Peak1D::fgFillStyleBkg = 3001;
85 Color_t Peak1D::fgMarkerColor = 36;
86 Color_t Peak1D::fgFuncColor = 26;
87 
88 
89 //__________________________________________________________
90 Peak1D::Peak1D() :
91  BasePeak("Peak","Gw::Peak1D"),
92  fPosition(0.0),
93  fFWHM(0.0),
94  fBkgLeft1(0.0),
95  fBkgLeft2(0.0),
96  fBkgRight1(0.0),
97  fBkgRight2(0.0),
98  fBkgLevelLeft1(0.0),
99  fBkgLevelLeft2(0.0),
100  fBkgLevelRight1(0.0),
101  fBkgLevelRight2(0.0),
102  fBkgIntegral(0.0),
103  fPeakIntegral(0.0),
104  fPeakIntegralError(0.0),
105  fSubPeakIntegral(0.0),
106  fBkgFlag(false),
107  fBinWidth(1.),
108  fPolyLinePeak(),
109  fPolyLineBkgL(),
110  fPolyLineBkgR(),
111  fCPolyline(0x0),
112  fCopyPolyline(),
113  fMarkerList(),
114  fDrawAs(kPeak),
115  fSigFunc(0x0),
116  fBkgFunc(0x0),
117  fPeakFunc(0x0),
118  fFitCombined(0x0),
119  fIsAGammaSourceRay(false),
120  fGammaSourceRay(0x0),
121  fLog("Peak1D"),
122  fEfficiencyGraph(0x0),
123  fEfficiencyFunc(0x0),
124  fRefArea(1),
125  fRefAreaError(0)
126 {
127  //Default Constructor
128  fDimension = 1;
129 
130  //reset polyline
131  Double_t x[fgkPointsPeak] = {0,0,0,0,0,0,0};
132  Double_t y[fgkPointsPeak] = {0,0,0,0,0,0,0};
133  fPolyLinePeak.SetPolyLine(fgkPointsPeak, x, y);
134  fPolyLinePeak.SetLineWidth(2);
135 
136  Double_t xb[fgkPointsBkg] = {0,0,0,0};
137  Double_t yb[fgkPointsBkg] = {0,0,0,0};
138  fPolyLineBkgL.SetPolyLine(fgkPointsBkg, xb, yb);
139  fPolyLineBkgL.SetLineWidth(2);
140 
141  fPolyLineBkgR.SetPolyLine(fgkPointsBkg, xb, yb);
142  fPolyLineBkgR.SetLineWidth(2);
143 
144  // set default line color
145  SetLineColorPeak(gGateColor.k0);
146  SetLineColorBkg(fgLineColorBkg);
147 
148  //set default fill style
149  SetFillStylePeak(fgFillStylePeak);
150  SetFillStyleBkg(fgFillStyleBkg);
151 
152  // set default fill color
153  SetFillColorPeak(gGateColor.k0);
154  SetFillColorBkg(fgLineColorBkg);
155 
156  // Add markers
157  InitPeakMarkers();
158 }
159 
160 //__________________________________________________________
161 Peak1D::Peak1D(const char* name, const char* title) :
162  BasePeak(name, title),
163  fPosition(0.0),
164  fFWHM(0.0),
165  fBkgLeft1(0.0),
166  fBkgLeft2(0.0),
167  fBkgRight1(0.0),
168  fBkgRight2(0.0),
169  fBkgLevelLeft1(0.0),
170  fBkgLevelLeft2(0.0),
171  fBkgLevelRight1(0.0),
172  fBkgLevelRight2(0.0),
173  fBkgIntegral(0.0),
174  fPeakIntegral(0.0),
175  fPeakIntegralError(0.0),
176  fSubPeakIntegral(0.0),
177  fBkgFlag(false),
178  fBinWidth(1.),
179  fPolyLinePeak(),
180  fPolyLineBkgL(),
181  fPolyLineBkgR(),
182  fCPolyline(0x0),
183  fCopyPolyline(),
184  fMarkerList(),
185  fDrawAs(kPeak),
186  fSigFunc(0x0),
187  fBkgFunc(0x0),
188  fPeakFunc(0x0),
189  fFitCombined(0x0),
190  fIsAGammaSourceRay(false),
191  fGammaSourceRay(0x0),
192  fLog("Peak1D"),
193  fEfficiencyGraph(0x0),
194  fEfficiencyFunc(0x0),
195  fRefArea(1),
196  fRefAreaError(0)
197 {
198  //Default Constructor
199  fDimension = 1;
200 
201  //reset polyline
202  Double_t x[fgkPointsPeak] = {0,0,0,0,0,0,0};
203  Double_t y[fgkPointsPeak] = {0,0,0,0,0,0,0};
204  fPolyLinePeak.SetPolyLine(fgkPointsPeak, x, y);
205  fPolyLinePeak.SetLineWidth(2);
206 
207  Double_t xb[fgkPointsBkg] = {0,0,0,0};
208  Double_t yb[fgkPointsBkg] = {0,0,0,0};
209  fPolyLineBkgL.SetPolyLine(fgkPointsBkg, xb, yb);
210  fPolyLineBkgL.SetLineWidth(2);
211 
212  fPolyLineBkgR.SetPolyLine(fgkPointsBkg, xb, yb);
213  fPolyLineBkgR.SetLineWidth(2);
214 
215  // set default line color
216  SetLineColorPeak(gGateColor.k0);
217  SetLineColorBkg(fgLineColorBkg);
218 
219  //set default fill style
220  SetFillStylePeak(fgFillStylePeak);
221  SetFillStyleBkg(fgFillStyleBkg);
222 
223  // set default fill color
224  SetFillColorPeak(gGateColor.k0);
225  SetFillColorBkg(fgLineColorBkg);
226 
227  // Add markers
228  InitPeakMarkers();
229 }
230 
231 //__________________________________________________________
232 Peak1D::Peak1D(TPolyLine* polyline) :
233  BasePeak("Peak","Gw::Peak1D"),
234  fPosition(0.0),
235  fFWHM(0.0),
236  fBkgLeft1(0.0),
237  fBkgLeft2(0.0),
238  fBkgRight1(0.0),
239  fBkgRight2(0.0),
240  fBkgLevelLeft1(0.0),
241  fBkgLevelLeft2(0.0),
242  fBkgLevelRight1(0.0),
243  fBkgLevelRight2(0.0),
244  fBkgIntegral(0.0),
245  fPeakIntegral(0.0),
246  fPeakIntegralError(0.0),
247  fSubPeakIntegral(0.0),
248  fBkgFlag(false),
249  fBinWidth(1.),
250  fPolyLinePeak(),
251  fPolyLineBkgL(),
252  fPolyLineBkgR(),
253  fCPolyline(0x0),
254  fCopyPolyline(),
255  fMarkerList(new TList()),
256  fDrawAs(kPeak),
257  fSigFunc(0x0),
258  fBkgFunc(0x0),
259  fPeakFunc(0x0),
260  fFitCombined(0x0),
261  fIsAGammaSourceRay(false),
262  fGammaSourceRay(0x0),
263  fLog("Peak1D"),
264  fEfficiencyGraph(0x0),
265  fEfficiencyFunc(0x0),
266  fRefArea(1),
267  fRefAreaError(0)
268 {
269  //Default Constructor
270  fDimension = 1;
271 
272  //set peak polyline
273  SetPeakPoints(polyline);
274  UpdatePeak();
275 
276  Double_t xb[fgkPointsBkg] = {0,0,0,0};
277  Double_t yb[fgkPointsBkg] = {0,0,0,0};
278  fPolyLineBkgL.SetPolyLine(fgkPointsBkg, xb, yb);
279  fPolyLineBkgL.SetLineWidth(2);
280 
281  fPolyLineBkgR.SetPolyLine(fgkPointsBkg, xb, yb);
282  fPolyLineBkgR.SetLineWidth(2);
283 
284  // set default line color
285  SetLineColorPeak(gGateColor.k0);
286  SetLineColorBkg(fgLineColorBkg);
287 
288  //set default fill style
289  SetFillStylePeak(fgFillStylePeak);
290  SetFillStyleBkg(fgFillStyleBkg);
291 
292  // set default fill color
293  SetFillColorPeak(gGateColor.k0);
294  SetFillColorBkg(fgLineColorBkg);
295 
296  // Add markers
297  InitPeakMarkers();
298 }
299 
300 
301 //__________________________________________________________
302 Peak1D::Peak1D(const Peak1D& p) : BasePeak(p),
303  fPosition(p.fPosition),
304  fFWHM(p.fFWHM),
305  fBkgLeft1(p.fBkgLeft1),
306  fBkgLeft2(p.fBkgLeft2),
307  fBkgRight1(p.fBkgRight1),
308  fBkgRight2(p.fBkgRight2),
309  fBkgLevelLeft1(p.fBkgLevelLeft1),
310  fBkgLevelLeft2(p.fBkgLevelLeft2),
311  fBkgLevelRight1(p.fBkgLevelRight1),
312  fBkgLevelRight2(p.fBkgLevelRight2),
313  fBkgIntegral(p.fBkgIntegral),
314  fPeakIntegral(p.fPeakIntegral),
315  fPeakIntegralError(p.fPeakIntegralError),
316  fSubPeakIntegral(p.fSubPeakIntegral),
317  fBkgFlag(p.fBkgFlag),
318  fBinWidth(p.fBkgFlag),
319  fDrawAs(p.fDrawAs),
320  fIsAGammaSourceRay(p.fIsAGammaSourceRay),
321  fLog("Peak1D"),
322  fEfficiencyGraph(p.fEfficiencyGraph),
323  fEfficiencyFunc(p.fEfficiencyFunc),
324  fRefArea(p.fRefArea),
325  fRefAreaError(p.fRefArea)
326 
327 {
328  // copy constructor
329  fPolyLinePeak.Copy((TPolyLine&)p.fPolyLinePeak);
330  fPolyLineBkgL.Copy((TPolyLine&)p.fPolyLineBkgL);
331  fPolyLineBkgR.Copy((TPolyLine&)p.fPolyLineBkgR);
332  fMarkerList.Copy((TList&)p.fMarkerList);
333 
334  if ( fSigFunc )
335  delete fSigFunc;
336  fSigFunc = new TF1(*p.fSigFunc);
337  if ( fBkgFunc )
338  delete fBkgFunc;
339  fBkgFunc = new TF1(*p.fBkgFunc);
340  if ( fPeakFunc )
341  delete fPeakFunc;
342  fPeakFunc = new TF1(*p.fPeakFunc);
343  if ( fFitCombined )
344  delete fFitCombined;
345  fFitCombined = new TF1(*p.fFitCombined);
346 
347  // fGammaSourceRay = new Gw::GammaSourceRay(*p.fGammaSourceRay);
348 
349  InitPeakMarkers();
350 }
351 
352 //__________________________________________________________
353 void Peak1D::Copy(TObject &p) const
354 {
355  BasePeak::Copy(p);
356 
357  ((Peak1D&)p).fPosition = fPosition;
358  ((Peak1D&)p).fFWHM = fFWHM;
359  ((Peak1D&)p).fBkgLeft1 = fBkgLeft1;
360  ((Peak1D&)p).fBkgLeft2 = fBkgLeft2;
361  ((Peak1D&)p).fBkgRight1 = fBkgRight1;
362  ((Peak1D&)p).fBkgRight2 = fBkgRight2;
363  ((Peak1D&)p).fBkgLevelLeft1 = fBkgLevelLeft1;
364  ((Peak1D&)p).fBkgLevelLeft2 = fBkgLevelLeft2;
365  ((Peak1D&)p).fBkgLevelRight1 = fBkgLevelRight1;
366  ((Peak1D&)p).fBkgLevelRight2 = fBkgLevelRight2;
367  ((Peak1D&)p).fBkgIntegral = fBkgIntegral;
368  ((Peak1D&)p).fPeakIntegral = fPeakIntegral;
369  ((Peak1D&)p).fPeakIntegralError = fPeakIntegralError;
370  ((Peak1D&)p).fSubPeakIntegral = fSubPeakIntegral;
371  ((Peak1D&)p).fBkgFlag = fBkgFlag;
372  ((Peak1D&)p).fBinWidth = fBinWidth;
373  ((Peak1D&)p).fDrawAs = fDrawAs;
374  ((Peak1D&)p).fIsAGammaSourceRay = fIsAGammaSourceRay;
375  ((Peak1D&)p).fEfficiencyGraph = fEfficiencyGraph;
376  ((Peak1D&)p).fEfficiencyFunc = fEfficiencyFunc;
377  ((Peak1D&)p).fRefArea = fRefArea;
378  ((Peak1D&)p).fRefAreaError = fRefAreaError;
379 
380  fPolyLinePeak.Copy( ((Peak1D &)p).fPolyLinePeak);
381  fPolyLineBkgL.Copy( ((Peak1D &)p).fPolyLineBkgL);
382  fPolyLineBkgR.Copy( ((Peak1D &)p).fPolyLineBkgR);
383 
384  // (T(List&)p.fMarkerList).Copy(fMarkerList);
385 
386  if ( ((Peak1D&)p).fSigFunc )
387  delete ((Peak1D&)p).fSigFunc;
388  ((Peak1D&)p).fSigFunc = new TF1(*fSigFunc);
389  if ( ((Peak1D&)p).fBkgFunc )
390  delete ((Peak1D&)p).fBkgFunc;
391  ((Peak1D&)p).fBkgFunc = new TF1(*fBkgFunc);
392  if ( ((Peak1D&)p).fPeakFunc )
393  delete ((Peak1D&)p).fPeakFunc;
394  ((Peak1D&)p).fPeakFunc = new TF1(*fPeakFunc);
395  if ( ((Peak1D&)p).fFitCombined )
396  delete ((Peak1D&)p).fFitCombined;
397  ((Peak1D&)p).fFitCombined = new TF1(*fFitCombined);
398 
399  // if ( ((Peak1D&)p).fGammaSourceRay )
400  // delete ((Peak1D&)p).fGammaSourceRay;
401  // ((Peak1D&)p).fGammaSourceRay = new Gw::GammaSourceRay(*fGammaSourceRay);
402 
403  ((Peak1D&)p).InitPeakMarkers();
404 
405 }
406 
407 /*
408 //__________________________________________________________
409 Peak1D& Peak1D::operator=(const Peak1D& p)
410 {
411  // assignment operator
412  if (this == &p)
413  return *this;
414  BasePeak::operator=(p);
415 
416  fPolyLinePeak = p.fPolyLinePeak;
417  fPolyLineBkgL = p.fPolyLineBkgL;
418  fPolyLineBkgR = p.fPolyLineBkgR;
419  fMarkerList = p.fMarkerList;
420  fSigFunc = p.fSigFunc;
421  fBkgFunc = p.fBkgFunc;
422  fPeakFunc = p.fPeakFunc;
423  fFitCombined = p.fFitCombined;
424  fPosition = p.fPosition;
425  fFWHM = p.fFWHM;
426  fBkgLeft1 = p.fBkgLeft1;
427  fBkgLeft2 = p.fBkgLeft2;
428  fBkgRight1 = p.fBkgRight1;
429  fBkgRight2 = p.fBkgRight2;
430  fBkgLevelLeft1 = p.fBkgLevelLeft1;
431  fBkgLevelLeft2 = p.fBkgLevelLeft2;
432  fBkgLevelRight1 = p.fBkgLevelRight1;
433  fBkgLevelRight2 = p.fBkgLevelRight2;
434  fBkgFlag = p.fBkgFlag;
435  fDrawAs = p.fDrawAs;
436 
437  return *this;
438 }
439  */
440 
441 //__________________________________________________________
442 Peak1D::~Peak1D()
443 {
444  // default destructor
445  fMarkerList.Delete();
446 
447  if(fSigFunc) delete fSigFunc;
448  if(fBkgFunc) delete fBkgFunc;
449  if(fPeakFunc) delete fPeakFunc;
450  if(fFitCombined) delete fFitCombined;
451 }
452 
453 
454 //__________________________________________________________
455 Double_t Peak1D::Bkg(Double_t* x, Double_t* par)
456 {
457  // reject point outside the backgorund limits
458  if ( x[0] < fBkgLeft1 || x[0] > fBkgRight2 ) {
459  TF1::RejectPoint();
460 
461  return 0;
462  }
463  if ( x[0] > fBkgLeft2 && x[0] < fBkgRight1 ) {
464  TF1::RejectPoint();
465 
466  return 0;
467  }
468  return fBkgFunc->EvalPar(x, par);
469 }
470 
471 //__________________________________________________________
472 Double_t Peak1D::SignalAndBkg(Double_t* x, Double_t* par)
473 {
474  // Total function
475  Double_t f = fSigFunc->EvalPar(x, par);
476 
477  if ( fBkgFunc )
478  f += fBkgFunc->EvalPar(x, par + fSigFunc->GetNpar());
479 
480  return f;
481 }
482 
483 TF1 *Peak1D::SignalAndBkgFunction()
484 {
485  TString FuncName = fSigFunc->GetName();
486 
487  if(FuncName.Contains("DTGaus") && FuncName.Contains("Step"))
488  return fPeakFunc;
489 
490  // if already set
491  if (fPeakFunc)
492  delete fPeakFunc;
493 
494  // now set the global function
495  Double_t *allpar = new Double_t[fSigFunc->GetNpar()+fBkgFunc->GetNpar()];
496 
497  ::memcpy(allpar, fSigFunc->GetParameters(), fSigFunc->GetNpar()*sizeof(Double_t));
498  ::memcpy(allpar + fSigFunc->GetNpar(), fBkgFunc->GetParameters(), fBkgFunc->GetNpar()*sizeof(Double_t));
499 
500  TString tmp = Form("SigAndBkg_%p",this);
501  fPeakFunc = new TF1(tmp.Data(), this, &Peak1D::SignalAndBkg, fBkgLeft1, fBkgRight2, fSigFunc->GetNpar()+fBkgFunc->GetNpar(), "Peak1D", "SignalAndBkg");
502  fPeakFunc->SetParameters(allpar);
503  delete allpar;
504 
505  fPeakFunc->SetBit(TObject::kCannotPick);
506  gROOT->GetListOfFunctions()->Remove(fPeakFunc);
507 
508  fPeakFunc->SetLineColor(fgFuncColor);
509  return fPeakFunc;
510 }
511 
512 //__________________________________________________________
513 void Peak1D::WithBkg(Bool_t with_bg)
514 {
515  fBkgFlag = with_bg;
516 }
517 
518 //__________________________________________________________
519 void Peak1D::SetPeak(Double_t position, Double_t height, Double_t FWHM, Double_t intensity)
520 {
521  // constructor for members
522  fPosition = position;
523  fFWHM = FWHM;
524  fHeight = height;
525  fIntensity = intensity;
526 
527  SetModified(fPolyLinePeak);
528 }
529 
530 void Peak1D::SetPeakFromFit()
531 {
532  fHeight = fSigFunc->GetParameter(2);
533  fPosition = fSigFunc->GetParameter(3);
534  fFWHM = fSigFunc->GetX(fHeight/2.,fPosition,fPosition+5*fFWHM)-fSigFunc->GetX(fHeight/2.,fPosition-5*fFWHM,fPosition);
535 #if ROOT_VERSION_CODE > ROOT_VERSION(5,34,36)
536  fIntensity = fSigFunc->Integral(fBkgLeft1,fBkgRight2,1e-6)/fBinWidth;
537 #else
538  fIntensity = fSigFunc->Integral(fBkgLeft1,fBkgRight2)/fBinWidth;
539 #endif
540 
541  // modified
542  SetPeakPoints();
543 }
544 
545 //__________________________________________________________
546 void Peak1D::SetPosition(Double_t position, Option_t* /*axis*/)
547 {
548  // set position
549  fPosition = position;
550 
551  SetModified(fPolyLinePeak);
552 }
553 
554 //__________________________________________________________
555 void Peak1D::SetFWHM(Double_t FWHM, Option_t* /*axis*/)
556 {
557  if ( FWHM <= 0.0 )
558  return;
559 
560  // set FWHM
561  fFWHM = FWHM;
562 
563  SetModified(fPolyLinePeak);
564 }
565 
566 //__________________________________________________________
567 void Peak1D::SetHeight(Double_t height)
568 {
569  if ( height <= 0.0 )
570  return;
571 
572  // set height
573  BasePeak::SetHeight(height);
574 
575  SetModified(fPolyLinePeak);
576 }
577 
578 //__________________________________________________________
579 void Peak1D::SetIntensity(Double_t intensity)
580 {
581  if ( intensity <= 0.0 )
582  return;
583 
584  // set height
585  BasePeak::SetIntensity(intensity);
586 
587  SetModified(fPolyLinePeak);
588 }
589 
590 //__________________________________________________________
591 void Peak1D::SetBackground( Double_t bgLeft1, Double_t bgLeft2, Double_t bgRight1, Double_t bgRight2,
592  Double_t bgLevelLeft1, Double_t bgLevelLeft2, Double_t bgLevelRight1,Double_t bgLevelRight2)
593 {
594  // set bkg limits left & right from the peak
595  fBkgLeft1 = bgLeft1;
596  fBkgLeft2 = bgLeft2;
597  fBkgRight1 = bgRight1;
598  fBkgRight2 = bgRight2;
599 
600  fBkgLevelLeft1 = bgLevelLeft1;
601  fBkgLevelLeft2 = bgLevelLeft2;
602  fBkgLevelRight1 = bgLevelRight1;
603  fBkgLevelRight2 = bgLevelRight2;
604 
605  SetModified(fPolyLineBkgL);
606  SetModified(fPolyLineBkgR);
607 }
608 
609 void Peak1D::SetBackground( Double_t bgLeft1, Double_t bgLeft2, Double_t bgRight1, Double_t bgRight2, const TH1 *h )
610 {
611  SetBackground(bgLeft1,
612  bgLeft2,
613  bgRight1,
614  bgRight2,
615  h->GetBinContent(h->GetXaxis()->FindFixBin(bgLeft1)),
616  h->GetBinContent(h->GetXaxis()->FindFixBin(bgLeft2)),
617  h->GetBinContent(h->GetXaxis()->FindFixBin(bgRight1)),
618  h->GetBinContent(h->GetXaxis()->FindFixBin(bgRight2)));
619 }
620 
621 //__________________________________________________________
622 void Peak1D::ShiftPolyline(TPolyLine &poly, Double_t x, Double_t y)
623 {
624  for (Int_t i = 0; i < poly.GetN(); i++) {
625  poly.GetX()[i] += x;
626  poly.GetY()[i] += y;
627  }
628 }
629 
630 //__________________________________________________________
631 void Peak1D::SetAllPointsinPolyline(TPolyLine &poly, Double_t x, Double_t y)
632 {
633  for (Int_t i = 0; i < poly.GetN(); i++) {
634  poly.GetX()[i] = x;
635  poly.GetY()[i] = y;
636  }
637 }
638 
639 //__________________________________________________________
640 void Peak1D::SetMarkerColor(Int_t color)
641 {
642  // set marker color
643  TMarker* m = 0x0;
644  for (Int_t i = 0; i < fMarkerList.GetEntries(); ++i) {
645  m = static_cast<TMarker*> (fMarkerList.At(i));
646  m->SetMarkerColor(color);
647  }
648 }
649 
650 //__________________________________________________________
651 TF1 *Peak1D::SetSignalFunction(const TF1* func)
652 {
653  // set user signal function
654 
655  fLog.GetProcessMethod() = "SetSignalFunction(const TF1*)";
656 
657  if (func != 0x0) { //user defined function
658 
659  Bool_t ok = true;
660 
661  // check it is an expected kind of function
662  TString name = func->GetParName(0);
663  name.ToLower();
664  if (name != "height") {
665  fLog << error << "Wrong name or bad index for parameter 0 (height expected)" << nline;
666  ok = false;
667  }
668  name = func->GetParName(1);
669  name.ToLower();
670  if (name != "position") {
671  fLog << error << "Wrong name or bad index for parameter 1 (position expected)" << nline;
672  ok = false;
673  }
674  name = func->GetParName(2);
675  name.ToLower();
676  if (name != "fwhm") {
677  fLog << error << "Wrong name or bad index for parameter 2 (FWHM expected)" << nline;
678  ok = false;
679  }
680  if ( ok == false ) {
681  fLog << dolog;
682 
683  return fSigFunc;
684  }
685 
686  // try to build the new function
687  TF1 *tmp = new TF1(*func);
688 
689  name = Form("%s_Sig_%p",func->GetName(),this);
690  tmp->SetName(name.Data());
691 
692  // check if function is correct
693 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,4)
694  if ( !tmp->IsValid() ) {
695 #else
696  if ( tmp->Compile() ) {
697 #endif
698  fLog << error << "Bad formula for function" << dolog;
699  delete tmp;
700 
701  return fSigFunc;
702  }
703  if (fSigFunc) // delete previous fit
704  delete fSigFunc;
705 
706  fSigFunc = tmp;
707  gROOT->GetListOfFunctions()->Remove(fSigFunc);
708 
709  } else fLog << error << "No function defined" << nline;
710 
711  fLog << dolog;
712 
713  return fSigFunc;
714 }
715 //__________________________________________________________
716 TF1 *Peak1D::SetSignalFunction(const char* nameFunc, TH1 *h)
717 {
718  // set pre-defined signal function
719  Double_t xmin = 0, xmax = 0;
720 
721  // check first it has not yet been set with the same type
722  TString name = Form("%s_Sig_%p",nameFunc,this);
723 
724  // check if it is needed to creat a new function
725  Bool_t do_new_func = true;
726  if ( fSigFunc ) {
727  if ( name == fSigFunc->GetName() ) {
728  do_new_func = false;
729  }
730  else {
731  delete fSigFunc; fSigFunc = 0x0;
732  }
733  }
734  if (do_new_func ) {
735  if (name.Contains("gaus")) { // gaus function
736 
737  xmin = fPosition - 3*fFWHM; xmax = fPosition + 3*fFWHM;
738 
739  fSigFunc = new TF1(name.Data(),this,&Peak1D::Gaus,xmin,xmax,5,"Peak1D","Gaus");
740  fSigFunc->SetParName(0, "NumberOfPeaks");
741  fSigFunc->SetParName(1, "UsingBkg");
742  fSigFunc->SetParName(2, "Height");
743  fSigFunc->SetParName(3, "Position");
744  fSigFunc->SetParName(4, "FWHM");
745 
746  fSigFunc->SetBit(TObject::kCannotPick);
747  gROOT->GetListOfFunctions()->Remove(fSigFunc);
748  }
749  else if(name.Contains("land"))
750  {
751  xmin = fPosition - 3*fFWHM; xmax = fPosition + 10*fFWHM;
752 
753  fSigFunc = new TF1(name.Data(),this,&Peak1D::Landau,xmin,xmax,5,"Peak1D","Landau");
754 
755  fSigFunc->SetBit(TObject::kCannotPick);
756  gROOT->GetListOfFunctions()->Remove(fSigFunc);
757  }
758  else if(name.Contains("DSCB"))
759  {
760  xmin = fPosition - 5*fFWHM; xmax = fPosition + 5*fFWHM;
761 
762  fSigFunc = new TF1(name.Data(),this,&Peak1D::DoubleSidedCrystalBallFct,xmin,xmax,9,"Peak1D","DoubleSidedCrystalBallFct");
763 
764  fSigFunc->SetParName(5, "a1");
765  fSigFunc->SetParName(6, "p1");
766  fSigFunc->SetParName(7, "a2");
767  fSigFunc->SetParName(8, "p2");
768 
769  fSigFunc->SetBit(TObject::kCannotPick);
770  gROOT->GetListOfFunctions()->Remove(fSigFunc);
771  }
772  else if(name.Contains("DTGaus"))
773  {
774  xmin = fPosition - 5*fFWHM; xmax = fPosition + 5*fFWHM;
775 
776  if(name.Contains("Step"))
777  {
778  fSigFunc = new TF1(name.Data(),this,&Peak1D::DoubleTailedStepedGaussian,xmin,xmax,8,"Peak1D","DoubleTailedStepedGaussian");
779  fSigFunc->SetParName(7, "AmplitudeStep");
780  }
781  else
782  {
783  fSigFunc = new TF1(name.Data(),this,&Peak1D::DoubleTailedGaussian,xmin,xmax,7,"Peak1D","DoubleTailedGaussian");
784  }
785 
786  fSigFunc->SetParName(5, "LeftTail");
787  fSigFunc->SetParName(6, "RightTail");
788 
789  fSigFunc->SetBit(TObject::kCannotPick);
790  gROOT->GetListOfFunctions()->Remove(fSigFunc);
791  }
792  else
793  {
794  fLog.GetProcessMethod() = "SetSignalFunction(const char* )";
795  fLog << error << "Function not defined" << dolog;
796 
797  return fSigFunc;
798  }
799  }
800 
801  fSigFunc->SetParName(0, "NumberOfPeaks");
802  fSigFunc->SetParName(1, "UsingBkg");
803  fSigFunc->SetParName(2, "Height");
804  fSigFunc->SetParName(3, "Position");
805  fSigFunc->SetParName(4, "FWHM");
806 
807  fSigFunc->FixParameter(0, 1); // 1 peak
808  fSigFunc->FixParameter(1, 0); // No background
809  fSigFunc->SetParameter(2, fHeight);
810  fSigFunc->SetParameter(3, fPosition);
811  fSigFunc->SetParLimits(3, xmin, xmax);
812  fSigFunc->SetParameter(4, fFWHM);
813  fSigFunc->SetParLimits(4, fFWHM/2.,fFWHM*4);
814 
815  fSigFunc->SetRange(xmin,xmax);
816 
817  if ( fSigFunc && (name.Contains("DSCB") ) )
818  {
819  fSigFunc->SetParameter(5, 2.);
820  fSigFunc->SetParLimits(5, 0., 5.);
821  fSigFunc->SetParameter(6, 5.);
822  fSigFunc->SetParLimits(6, 0., 25.);
823  fSigFunc->SetParameter(7, 2.);
824  fSigFunc->SetParLimits(7, 0., 5.);
825  fSigFunc->SetParameter(8, 5.);
826  fSigFunc->SetParLimits(8, 0., 25.);
827  }
828  else if ( fSigFunc && (name.Contains("DTGaus") ) )
829  {
830  fSigFunc->SetParameter(5, -2);
831  fSigFunc->SetParLimits(5,-50,-0.1);
832  fSigFunc->FixParameter(6, 100);
833  // fSigFunc->SetParameter(6, 2);
834  // fSigFunc->SetParLimits(6,0.1,50.);
835 
836  if(name.Contains("Step"))
837  {
838  fSigFunc->SetParName(1, "Bgd_Amp");
839  fSigFunc->ReleaseParameter(1);
840  if(h!=0x0)
841  fSigFunc->SetParameter(1, h->GetBinContent(TMath::Min(h->FindBin(fPosition-4*fFWHM),h->FindBin(fPosition+4*fFWHM))));
842  else
843  fSigFunc->SetParameter(1, 0);
844 
845  fSigFunc->SetParameter(7, 0.);
846  fSigFunc->SetParLimits(7, 0., 1.);
847  }
848  }
849 
850  return fSigFunc;
851 }
852 
853 //__________________________________________________________
854 TF1 *Peak1D::SetBkgFunction(const TF1* func)
855 {
856  // set user background function
857  // the user has to provide the initial parameters
858  // fit is done within the background limits
859 
860  fLog.GetProcessMethod() = "SetBkgFunction(const TF1* )";
861 
862  if (func != 0x0) { //user defined function
863 
864  // try to build the new function
865  TF1 *tmp = (TF1*)func->Clone();
866  tmp->SetName(Form("%s_Bkg_%p",func->GetName(),this));
867 
868  // check if function is correct
869 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,4)
870  if ( !tmp->IsValid() ) {
871 #else
872  if ( tmp->Compile() ) {
873 #endif
874  fLog << error << "Bad formula for function" << dolog;
875  delete tmp;
876 
877  return fBkgFunc;
878  }
879 
880  if (fBkgFunc) // delete previous fit
881  delete fBkgFunc;
882 
883  fBkgFunc = tmp;
884  fBkgFunc->SetBit(TObject::kCannotPick);
885  gROOT->GetListOfFunctions()->Remove(fBkgFunc);
886 
887  } else fLog << error << "No function defined" << nline;
888 
889  fLog << dolog;
890 
891  return fBkgFunc;
892 }
893 
894 //__________________________________________________________
895 TF1 *Peak1D::SetBkgFunction(const char* nameFunc)
896 {
897  fLog.GetProcessMethod() = "SetBkgFunction(const char*)";
898 
899  TString namef = Form("%s_Bkg_%p",nameFunc,this);
900 
901  // check first if the function has not yet been allocated with the right type
902  // or if it could be modified
903  Bool_t do_new_func = true;
904  if ( fBkgFunc ) {
905  TString tmp = Form("%s_Bkg_%p",fBkgFunc->GetName(),this);
906 
907  if ( tmp == namef ) {
908  do_new_func = false;
909  }
910  else { delete fBkgFunc; fBkgFunc = 0x0; }
911  }
912  // no need to change it !
913  if ( do_new_func == false )
914  return fBkgFunc;
915 
916  // list of pre defined functions
917  if ( namef.Contains("cst") || namef.Contains("pol0") || ((TString)nameFunc) == "" || namef.Contains("-") ) {
918  fBkgFunc = new TF1(namef.Data(), "pol0", fBkgLeft1, fBkgRight2);
919  fBkgFunc->SetParameter(0,0.0); // so that is can be used in case no background substraction is required
920  }
921 
922  if ( namef.Contains("lin") || namef.Contains("pol1") )
923  fBkgFunc = new TF1(namef.Data(), "pol1", fBkgLeft1, fBkgRight2);
924 
925  if ( namef.Contains("quad") || namef.Contains("pol2") )
926  fBkgFunc = new TF1(namef.Data(), "pol2", fBkgLeft1, fBkgRight2);
927  // to do ... step function
928 
929  // try directly to compile the given expression
930  if ( fBkgFunc == 0x0 ) {
931  fBkgFunc = new TF1(namef.Data(), nameFunc, fBkgLeft1, fBkgRight2);
932  // check if function is correct
933 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,3,4)
934  if ( !fBkgFunc->IsValid() ) {
935 #else
936  if ( fBkgFunc->Compile() ) {
937 #endif
938  fLog << error << "Bad formula for function" << dolog;
939  delete fBkgFunc; fBkgFunc = 0x0;
940  }
941  }
942  if ( fBkgFunc == 0x0 ) {
943  fBkgFunc = new TF1("pol0_Bkg", "pol0", fBkgLeft1, fBkgRight2);
944  fBkgFunc->SetParameter(0,0.0);
945  fLog << warning
946  << "Cannot define " << namef.Data() << " , pol0 is set instead" << dolog;
947  }
948  fBkgFunc->SetBit(TObject::kCannotPick);
949  gROOT->GetListOfFunctions()->Remove(fBkgFunc);
950 
951  return fBkgFunc;
952 }
953 
954 //__________________________________________________________
955 Double_t Peak1D::GetPosition(Option_t* /*axis*/) const
956 {
957  // return position for given axis
958  return fPosition;
959 }
960 
961 //__________________________________________________________
962 Double_t Peak1D::GetFWHM(Option_t* /*axis*/) const
963 {
964  // return FWHM for given axis
965  return fFWHM;
966 }
967 
968 //__________________________________________________________
969 TMarker* Peak1D::GetMarker(Int_t markerId)
970 {
971  // get marker
972  if ( markerId >= fMarkerList.GetEntries() ) {
973  fLog.GetProcessMethod() = "GetMarker(Int_t )";
974  fLog << error << "Index out of range" << dolog;
975  }
976 
977  return ( static_cast<TMarker*> (fMarkerList.At(markerId)) );
978 }
979 
980 //__________________________________________________________
981 void Peak1D::FitMenu(const char* nameFunc, Option_t* optFit, Option_t* optBkg)
982 {
983  fLog.SetProcessMethod("FitMenu(const char* nameFunc, Option_t* optFit, Option_t* optBkg)");
984 
985  // check the current pad owns the peak
986  if ( TVirtualPad::Pad() == 0x0 )
987  return;
988  // check the current pad owns the peak
989  if ( TVirtualPad::Pad()->FindObject(this) == 0x0 )
990  return;
991 
992  TH1 *h = PadManager::GetHisto(TVirtualPad::Pad());
993  if (h == 0x0) {
994  fLog << error << "Histogram not defined" << dolog;
995  return;
996  }
997 
998  // set signal function
999  SetSignalFunction(nameFunc,h);
1000 
1001  // look for other peaks in the pad in case fit combined is required
1002  std::vector<Gw::Peak1D *> allpeaks;
1003  PadManager::Collect<Gw::Peak1D>(allpeaks,TVirtualPad::Pad());
1004 
1005  TList peaks; peaks.SetOwner(false);
1006  //
1007  for (size_t i = 0; i < allpeaks.size(); i++) {
1008  if ( allpeaks[i] == this )
1009  continue;
1010  if ( IsPointInBkg(allpeaks[i]->GetPosition()) < 1 ) {
1011  allpeaks[i]->SetSignalFunction(nameFunc);
1012  peaks.Add(allpeaks[i]);
1013  }
1014  }
1015  if ( peaks.GetSize() ) {
1016  fLog << " COMBINE called because several peaks are between this peaks limits : # "<< peaks.GetSize() << nline;
1017  FitCombined(h, &peaks, optFit, optBkg);
1018  }
1019  else
1020  Fit(h, optFit, optBkg);
1021  //
1022  TIter next(&peaks); TObject *obj;
1023  while ( (obj = next()) ) {
1024  peaks.Remove(obj);
1025  }
1026  // send to log
1027  fLog << dolog;
1028 }
1029 
1030 //__________________________________________________________
1031 void Peak1D::Area(Option_t * opt)
1032 {
1033  TString Opt(opt); TH1 *h = 0x0; Short_t choice = 0;
1034 
1035  fLog.SetProcessMethod("Area"); // fLog << warning << "Not yet implemented" << dolog;
1036 
1037  if ( Opt.Contains("h") ) {
1038  if ( ( h = PadManager::GetHisto(TVirtualPad::Pad()) ) ) // only if there is an histogram in the current pad
1039  choice = 1;
1040  }
1041 
1042  if ( Opt.Contains("f") ) {
1043  if (fSigFunc != 0x0 && fBkgFunc != 0x0)
1044  choice = 2;
1045  }
1046 
1047  if ( Opt.Contains("g") ) { // get area when peak is drawn as a gate
1048  if ( ( h = PadManager::GetHisto(TVirtualPad::Pad()) ) ) // only if there is an histogram in the current pad
1049  choice = 3;
1050  }
1051  Double_t position = 0.0, sig = 0.0, sig_bg = 0.0, bg1 = 0.0, bg2 = 0.0, bg_sig = 0.0, sig_err=0.;
1052 
1053  switch ( choice ) {
1054  case 1:
1055  { //
1056  sig_bg = h->Integral(h->GetXaxis()->FindFixBin(fBkgLeft2),h->GetXaxis()->FindFixBin(fBkgRight1));
1057  //
1058  bg1 = h->Integral(h->GetXaxis()->FindFixBin(fBkgLeft1),h->GetXaxis()->FindFixBin(fBkgLeft2));
1059  bg2 = h->Integral(h->GetXaxis()->FindFixBin(fBkgRight1),h->GetXaxis()->FindFixBin(fBkgRight2));
1060 
1061  // loop on bins into bg limits to compute sig by substracting an extrapolation of bg
1062  Double_t slope = (fBkgLevelLeft2-fBkgLevelLeft1) / (fBkgLeft2-fBkgLeft1);
1063  Double_t offset = fBkgLevelLeft1 - slope*fBkgLeft1;
1064  //
1065  Double_t sum = 0.0;
1066  position = 0.0;
1067  sig = 0.0;
1068  for (Int_t i = h->GetXaxis()->FindFixBin(fBkgLeft2); i < h->GetXaxis()->FindFixBin(fBkgRight1); i++) {
1069 
1070  Double_t bin_cont = h->GetBinContent(i);
1071  Double_t bin_cent = h->GetBinCenter(i);
1072 
1073  sig += ( bin_cont - bin_cent*slope - offset );
1074  bg_sig += bin_cent*slope + offset;
1075  position += bin_cont*bin_cent;
1076  sum += bin_cont;
1077 
1078  }
1079  if ( sum )
1080  position /= sum;
1081  else {
1082  position = 0.0;
1083  }
1084  }
1085  break;
1086 
1087  case 2:
1088  {
1089 #if ROOT_VERSION_CODE > ROOT_VERSION(5,34,36)
1090  sig = fSigFunc->Integral(fPosition-4*fFWHM, fPosition+4*fFWHM,1e-6)/fBinWidth;
1091 #else
1092  sig = fSigFunc->Integral(fPosition-4*fFWHM, fPosition+4*fFWHM)/fBinWidth;
1093 #endif
1094 
1095  sig_err = 2*sqrt(sig);
1096 
1097  // if(Opt.Contains("e")) sig_err = 2*fSigFunc->IntegralError(fPosition-4*fFWHM, fPosition+4*fFWHM)/fBinWidth;
1098  // else
1099  // {
1100  // Float_t e1 = fSigFunc->GetParameter(4)*fSigFunc->GetParError(2);
1101  // Float_t e2 = fSigFunc->GetParameter(2)*fSigFunc->GetParError(4);
1102 
1103  // sig_err = TMath::Max(2*sqrt(2.*TMath::Pi())*sqrt(e1*e1+e2*e2),2*sqrt(sig));
1104  // }
1105  sig_bg = fBkgFunc->Integral(fPosition-fFWHM/2., fPosition+fFWHM/2.)/fBinWidth + sig;
1106  bg1 = fBkgFunc->Integral(fBkgLeft1, fBkgLeft2)/fBinWidth;
1107  bg2 = fBkgFunc->Integral(fBkgRight1, fBkgRight2)/fBinWidth;
1108  position = fPosition;
1109  }
1110  break;
1111 
1112  case 3:
1113  {
1114  //
1115  sig_bg = h->Integral(h->GetXaxis()->FindFixBin(fPosition-fFWHM/2.),h->GetXaxis()->FindFixBin(fPosition+fFWHM/2.));
1116  //
1117  bg1 = h->Integral(h->GetXaxis()->FindFixBin(fBkgLeft1),h->GetXaxis()->FindFixBin(fBkgLeft2));
1118  bg2 = h->Integral(h->GetXaxis()->FindFixBin(fBkgRight1),h->GetXaxis()->FindFixBin(fBkgRight2));
1119 
1120  // loop on bins into bg limits to compute sig by substracting an extrapolation of bg
1121  Double_t slope = (fBkgLevelLeft2-fBkgLevelLeft1) / (fBkgLeft2-fBkgLeft1);
1122  Double_t offset = fBkgLevelLeft1 - slope*fBkgLeft1;
1123  //
1124  Double_t sum = 0.0;
1125  position = 0.0;
1126  sig = 0.0;
1127  for (Int_t i = h->GetXaxis()->FindFixBin(fPosition-fFWHM/2.); i < h->GetXaxis()->FindFixBin(fPosition+fFWHM/2.); i++) {
1128 
1129  Double_t bin_cont = h->GetBinContent(i);
1130  Double_t bin_cent = h->GetBinCenter(i);
1131 
1132  sig += ( bin_cont - bin_cent*slope - offset );
1133  bg_sig += bin_cent*slope + offset;
1134  position += bin_cont*bin_cent;
1135  sum += bin_cont;
1136  }
1137  if ( sum )
1138  position /= sum;
1139  else {
1140  position = 0.0;
1141  }
1142  }
1143  break;
1144 
1145  default:
1146  break;
1147  }
1148  if ( Opt.Contains("v") )
1149  {
1150  fLog << setw(7)<<"Mean" <<setw(20)<<"Area_SIG"<<setw(20)<<"Area_SIG+BG" <<setw(20)<<"Area_BG" <<setw(20)<<"Area_BG_LEFT"<<setw(20)<<"Area_BG_RIGHT" << nline;
1151  fLog <<setprecision(6)<<setw(7)<<position <<setw(20)<<sig <<setw(20)<<sig_bg <<setw(20)<<sig_bg-sig <<setw(20)<<bg1<< setw(20)<<bg2 << nline;
1152  }
1153 
1154  fBkgIntegral = bg1+bg2;
1155  fPeakIntegral = sig;
1156  fPeakIntegralError = sig_err;
1157  fSubPeakIntegral = bg_sig;
1158 
1159  fLog << dolog;
1160 }
1161 
1162 
1163 //__________________________________________________________
1164 void Peak1D::FitCombined(TH1 *h, TList* peakList, Option_t* optFit, Option_t* optBkg)
1165 {
1166  fLog.GetProcessMethod() = "FitCombined(TH1 *h, Option_t*, Option_t*)";
1167 
1168  // just in case
1169  if ( peakList->GetSize() == 0 )
1170  return;
1171 
1172  // be sure the signal function is defined for all peaks
1173  Bool_t is_signals = true; Peak1D *peak;
1174 
1175  TIter next(peakList);
1176  while ( (peak = (Peak1D *)next()) ) {
1177  if ( peak->SignalFunction() == 0x0 ) {
1178  is_signals = false;
1179  break;
1180  }
1181  }
1182  if ( fSigFunc == 0x0 )
1183  is_signals = false;
1184  // cannot fit
1185  if ( is_signals == false ) {
1186  fLog << error << "At least one fitting function not defined, use SetSignalFunction() to set it" << dolog;
1187  return;
1188  }
1189 
1190  // deal with bg. if fit of bg required, the histogram is cloned to apply background substraction
1191  TString obkg(optBkg); TF1 *bkg = 0x0; TH1 *lochist = h;
1192  //Int_t nPar = fSigFunc->GetNpar();
1193 
1194  // set all peaks with the same background function
1195  next.Reset();
1196  while ( (peak = (Peak1D *)next() ) ) {
1197  peak->SetBkgFunction(optBkg);
1198  }
1199  bkg = SetBkgFunction(optBkg); // set the background function : no bg means pol0
1200 
1201  Bool_t do_bg_fit = true;
1202  if ( obkg == "" || obkg.Contains("-") ) {
1203  do_bg_fit = false;
1204  bkg->SetParameter(0,0.0); // in principle already done in SetBkgFunction
1205  }
1206  else {
1207  // a clone is needed to realize the background substraction
1208  lochist = (TH1*)h->Clone();
1209  }
1210 
1211  // local histo to perform the fit of the signal part in case of background substraction
1212  if ( do_bg_fit ) {
1213  // real function to fit excluding points outside the limits of the background
1214  bkg = new TF1("Peak_bkg", this, &Peak1D::Bkg, h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax(), fBkgFunc->GetNpar(), "Peak1D", "Bkg");
1215  // todo : add option for bkg fit i.e. optBkg = "lin[opt_fit==>pass to root]"
1216  // cout << " Backgroud fit" << endl;
1217  lochist->Fit(bkg,"RNQ");
1218  // cout << " Backgroud fit" << endl;
1219  // set parameter back to background function and set background level
1220  fBkgFunc->SetParameters(bkg->GetParameters());
1221  // no needed any more
1222  delete bkg;
1223 
1224  // substract bkg so that lochist is ready to be fitted by the signal
1225  lochist->Add(fBkgFunc, -1);
1226  }
1227  // update bkg part
1228  SetBackground(fBkgLeft1, fBkgLeft2, fBkgRight1, fBkgRight2,
1229  fBkgFunc->Eval(fBkgLeft1),
1230  fBkgFunc->Eval(fBkgLeft2),
1231  fBkgFunc->Eval(fBkgRight1),
1232  fBkgFunc->Eval(fBkgRight2));
1233  // update bkg part of other peaks
1234  next.Reset();
1235  while ( (peak = (Peak1D *)next() ) ) {
1236  peak->BkgFunction()->SetParameters(fBkgFunc->GetParameters());
1237  peak->SetBackground(fBkgLeft1, fBkgLeft2, fBkgRight1, fBkgRight2,
1238  fBkgFunc->Eval(fBkgLeft1),
1239  fBkgFunc->Eval(fBkgLeft2),
1240  fBkgFunc->Eval(fBkgRight1),
1241  fBkgFunc->Eval(fBkgRight2));
1242  }
1243 
1244  // functions to fit and draw
1245 
1246  if(fFitCombined != 0x0)
1247  delete fFitCombined;
1248 
1249  Int_t NPeaks = peakList->GetSize()+1;
1250  Int_t NPar = fSigFunc->GetNpar()-2;
1251 
1252  TString FuncName = fSigFunc->GetName();
1253 
1254  TF1 *f = (TF1*)gROOT->FindObjectAny(Form("CombinedPeak_%p",this));
1255  if(f != 0x0) delete f;
1256 
1257  if(FuncName.Contains("gaus"))
1258  fFitCombined = new TF1(Form("CombinedPeak_%p",this),this,&Peak1D::Gaus,fBkgLeft1,fBkgRight2,2+NPar*NPeaks,"Peak1D","Gaus");
1259  else if(FuncName.Contains("land"))
1260  fFitCombined = new TF1(Form("CombinedPeak_%p",this),this,&Peak1D::Landau,fBkgLeft1,fBkgRight2,2+NPar*NPeaks,"Peak1D","Landau");
1261  else if(FuncName.Contains("DSCB"))
1262  fFitCombined = new TF1(Form("CombinedPeak_%p",this),this,&Peak1D::DoubleSidedCrystalBallFct,fBkgLeft1,fBkgRight2,2+NPar*NPeaks,"Peak1D","DoubleSidedCrystalBallFct");
1263  else if(FuncName.Contains("DTGaus"))
1264  fFitCombined = new TF1(Form("CombinedPeak_%p",this),this,&Peak1D::DoubleTailedGaussian,fBkgLeft1,fBkgRight2,2+NPar*NPeaks,"Peak1D","DoubleTailedGaussian");
1265  else
1266  {
1267  fLog << error << "Unkown fit function : "<<FuncName <<" fit ignored !" << dolog;
1268  return;
1269  }
1270  fFitCombined->FixParameter(0,NPeaks);
1271  fFitCombined->FixParameter(1,0); // No Bkg for fit
1272  fFitCombined->SetLineColor(30);
1273 
1274  fFitCombined->SetParName(0,"NumberOfPeaks");
1275  fFitCombined->SetParName(1,"UsingBkg");
1276 
1277 
1278  next.Reset();
1279  Int_t which_peak = 0;
1280 
1281  peak = this;
1282  while ( which_peak ==0 || (peak = (Peak1D *)next() ) )
1283  {
1284  TF1 *PeakFunc = peak->SignalFunction();
1285 
1286  for(int i=0 ; i<NPar ; i++)
1287  {
1288  fFitCombined->SetParameter(2+NPar*which_peak+i,PeakFunc->GetParameter(i+2));
1289  double min, max;
1290  PeakFunc->GetParLimits(i+2,min,max);
1291  fFitCombined->SetParLimits(2+NPar*which_peak+i,min,max);
1292  fFitCombined->SetParName(2+NPar*which_peak+i,Form("%s_%d",PeakFunc->GetParName(i+2),which_peak));
1293  }
1294  which_peak++;
1295  }
1296 
1297  // deal with fit option
1298  TString roptFit(optFit);
1299  if ( ! roptFit.Contains("N") ) {
1300  roptFit += "N";
1301  }
1302  if(!roptFit.Contains("V"))
1303  {
1304  if ( !roptFit.Contains("Q")) {
1305  roptFit += "Q";
1306  }
1307  }
1308  // do fit
1309  SetBinWidth(lochist->GetBinWidth(1));
1310 
1311 
1312  if(FuncName.Contains("DTGaus"))
1313  {
1314 
1315  next.Reset();
1316  which_peak = 0;
1317  peak = this;
1318  while ( which_peak == 0 || (peak = (Peak1D *)next() ) )
1319  {
1320  fFitCombined->FixParameter(2+NPar*which_peak+3,-100);
1321  fFitCombined->FixParameter(2+NPar*which_peak+4,100);
1322 
1323  which_peak++;
1324  }
1325 
1326 // lochist->Fit(fFitCombined, roptFit.Data());
1327 
1328  next.Reset();
1329  which_peak = 0;
1330  peak = this;
1331  while ( which_peak == 0 || (peak = (Peak1D *)next() ) )
1332  {
1333  fFitCombined->SetParameter(2+NPar*which_peak+3,-2);
1334  fFitCombined->SetParLimits(2+NPar*which_peak+3,-50,-0.5);
1335 // fFitCombined->FixParameter(2+NPar*which_peak+3,-100);
1336 // fFitCombined->FixParameter(2+NPar*which_peak+4,100);
1337  fFitCombined->SetParameter(2+NPar*which_peak+4,2);
1338  fFitCombined->SetParLimits(2+NPar*which_peak+4,1.,50.);
1339 
1340  peak->SetPeakFromFit();
1341  which_peak++;
1342  }
1343  }
1344 
1345  lochist->Fit(fFitCombined, roptFit.Data());
1346 
1347 
1348  double IntErr = 0.;
1349  double TotInt = fFitCombined->Integral(fPosition-4*fFWHM, fPosition+4*fFWHM)/fBinWidth;
1350 
1351  if ( roptFit.Contains("E") ) IntErr = 2*fFitCombined->IntegralError(fPosition-4*fFWHM, fPosition+4*fFWHM)/fBinWidth;
1352  else IntErr = 2*sqrt(TotInt);
1353 
1354  fFitCombined->FixParameter(1,1.); // Background activated
1355 
1356  next.Reset();
1357  which_peak = 0;
1358 
1359  peak = this;
1360  while ( which_peak == 0 || (peak = (Peak1D *)next() ) )
1361  {
1362  TF1 *PeakFunc = peak->SignalFunction();
1363  TF1 *BckFunc = peak->BkgFunction();
1364 
1365  for (Int_t i = 0; i < fBkgFunc->GetNpar(); i++)
1366  {
1367  BckFunc->SetParameter(i,fBkgFunc->GetParameter(i));
1368  BckFunc->SetParError(i,fBkgFunc->GetParError(i));
1369  }
1370 
1371  PeakFunc->FixParameter(0,1);
1372  PeakFunc->FixParameter(1,0);
1373 
1374  Int_t NPar = fSigFunc->GetNpar()-2;
1375  for(int i=0 ; i<NPar ; i++)
1376  {
1377  PeakFunc->SetParameter(i+2,fFitCombined->GetParameter(2+NPar*which_peak+i));
1378  PeakFunc->SetParError(i+2,fFitCombined->GetParError(2+NPar*which_peak+i));
1379  }
1380 
1381  peak->SetBinWidth(fBinWidth);
1382  peak->SetPeakFromFit();
1383  peak->Area("f");
1384  peak->SetPeakIntegralError(IntErr*peak->GetPeakIntegral()/TotInt);
1385  peak->EnableFit();
1386  peak->SignalAndBkgFunction();
1387  peak->PeakFunction()->SetLineStyle(2);
1388  peak->PeakFunction()->SetLineColor(kBlue-6);
1389 
1390  // print on log
1391  roptFit = optFit;
1392  if ( !roptFit.Contains("Q") )
1393  peak->Print("f");
1394 
1395  which_peak++;
1396  }
1397 
1398  if ( lochist != h )
1399  delete lochist;
1400 
1401  fFitCombined->Draw("same");
1402  gROOT->GetListOfFunctions()->Remove( fFitCombined );
1403 }
1404 
1405 //__________________________________________________________
1406 void Peak1D::Fit(TH1 *h, Option_t* optFit, Option_t* optBkg)
1407 {
1408  fLog.GetProcessMethod() = "Fit(TH1 *h, Option_t*, Option_t* )";
1409 
1410  // be sure the signal function is defined
1411  if ( fSigFunc == 0x0 ) {
1412  fLog << error << "Fitting function not defined, use SetSignalFunction() to set it" << dolog;
1413  return;
1414  }
1415 
1416  // deal with bg. if fit of bg required, the histogram is cloned to apply background substraction
1417  TString obkg(optBkg); TF1 *bkg = 0x0; TH1 *lochist = h; Int_t nPar = fSigFunc->GetNpar();
1418 
1419  bkg = SetBkgFunction(optBkg); // set the background function : no bg means pol0
1420 
1421  Bool_t do_bg_fit = true;
1422  if ( obkg == "" || obkg.Contains("-") ) {
1423  do_bg_fit = false;
1424  bkg->SetParameter(0,0.0); // in principle already done in SetBkgFunction
1425  }
1426  else {
1427  // a clone is needed to realize the background substraction
1428  lochist = (TH1*)h->Clone();
1429  }
1430 
1431  // total number of parameters
1432  nPar += fBkgFunc->GetNpar();
1433 
1434  // local histo to perform the fit of the signal part in case of background substraction
1435  if ( do_bg_fit ) {
1436  // real function to fit excluding points outside the limits of the background
1437  bkg = new TF1("Peak_bkg", this, &Peak1D::Bkg, h->GetXaxis()->GetXmin(), h->GetXaxis()->GetXmax(), fBkgFunc->GetNpar(), "Peak1D", "Bkg");
1438  // todo : add option for bkg fit i.e. optBkg = "lin[opt_fit==>pass to root]"
1439  // cout << " Backgroud fit" << endl;
1440  lochist->Fit(bkg,"RNQ");
1441  // cout << " Backgroud fit" << endl;
1442  // set parameter back to background function and set background level
1443  fBkgFunc->SetParameters(bkg->GetParameters());
1444  // no needed any more
1445  delete bkg;
1446 
1447  // substract bkg so that lochist is ready to be fitted by the signal
1448  lochist->Add(fBkgFunc, -1);
1449  }
1450 
1451  // update bg part
1452  SetBackground(fBkgLeft1, fBkgLeft2, fBkgRight1, fBkgRight2,
1453  fBkgFunc->Eval(fBkgLeft1),
1454  fBkgFunc->Eval(fBkgLeft2),
1455  fBkgFunc->Eval(fBkgRight1),
1456  fBkgFunc->Eval(fBkgRight2));
1457 
1458  // deal with fit option
1459  TString roptFit(optFit);
1460  if ( ! roptFit.Contains("N") ) {
1461  roptFit += "N";
1462  }
1463  if(!roptFit.Contains("V"))
1464  {
1465  if ( !roptFit.Contains("Q")) {
1466  roptFit += "Q";
1467  }
1468  }
1469 
1470  TString namesave = fSigFunc->GetName();
1471 
1472  if(namesave.Contains("DTGaus") && namesave.Contains("Step"))
1473  fSigFunc->SetRange(fBkgLeft1,fBkgRight2);
1474 
1475  SetBinWidth(lochist->GetBinWidth(1));
1476 
1477 
1478  if(namesave.Contains("DTGaus"))
1479  {
1480  fSigFunc->FixParameter(5,-100);
1481  fSigFunc->FixParameter(6,100);
1482  fSigFunc->FixParameter(7,0.);
1483 
1484  lochist->Fit(fSigFunc, roptFit.Data());
1485  SetPeakFromFit();
1486 
1487  fSigFunc->SetParameter(5, -2);
1488  fSigFunc->SetParLimits(5,-50,-0.1);
1489  fSigFunc->SetParameter(6, 2);
1490  fSigFunc->SetParLimits(6,1.,50.);
1491 // fSigFunc->FixParameter(6,100);
1492 
1493  if(namesave.Contains("Step"))
1494  {
1495  fSigFunc->SetParameter(1, h->GetBinContent(TMath::Min(h->FindBin(fPosition-4*fFWHM),h->FindBin(fPosition+4*fFWHM))));
1496 
1497  fSigFunc->SetParameter(7, 0.);
1498  fSigFunc->SetParLimits(7, 0., 1.);
1499  }
1500  }
1501 
1502  lochist->Fit(fSigFunc, roptFit.Data());
1503 
1504  if(namesave.Contains("DTGaus") && namesave.Contains("Step"))
1505  {
1506  if(fPeakFunc != 0x0)
1507  delete fPeakFunc;
1508 
1509  if(fBkgFunc != 0x0)
1510  delete fBkgFunc;
1511 
1512  fBkgFunc = ExtractBkgFromDTGStep(Form("DTGausStepedBkg_%p",this));
1513 
1514  fPeakFunc = (TF1*) fSigFunc->Clone(Form("SigAndBkg_%p",this));
1515  fPeakFunc->SetBit(TObject::kCannotPick);
1516  gROOT->GetListOfFunctions()->Remove(fPeakFunc);
1517  fPeakFunc->SetLineColor(fgFuncColor);
1518 
1519  fSigFunc->SetParameter(1,0);
1520  fSigFunc->SetParameter(7,0);
1521 
1522  SetBackground(fBkgLeft1, fBkgLeft2, fBkgRight1, fBkgRight2,
1523  fBkgFunc->Eval(fBkgLeft1),
1524  fBkgFunc->Eval(fBkgLeft2),
1525  fBkgFunc->Eval(fBkgRight1),
1526  fBkgFunc->Eval(fBkgRight2));
1527  }
1528 
1529 
1530  // clone not needed anymore
1531  if ( lochist != h )
1532  delete lochist;
1533 
1534  // set signal+bkg function. Takes care of copy of the parameters
1535  SetPeakFromFit();
1536 
1537  if ( roptFit.Contains("E") ) Area("fe");
1538  else Area("f");
1539 
1540  SignalAndBkgFunction();
1541 
1542  // fit performed
1543  EnableFit();
1544 
1545  // print on log
1546  roptFit = optFit;
1547  if ( !roptFit.Contains("Q") )
1548  Print("f");
1549 
1550  // display polyline showing diff between h and fit
1551  /*
1552  TPolyLine *l = PadManager::ShowDiff(h,fPeakFunc,fBkgLeft1,fBkgRight2);
1553  l->SetLineStyle(2);
1554  l->SetLineColor(29); */
1555 
1556  if(gPad != 0x0 && !roptFit.Contains("0"))
1557  {
1558  gPad->Modified();
1559  gPad->Update();
1560  }
1561  fLog << dolog;
1562 }
1563 
1564 //___________________________________________________________
1565 void Peak1D::FitIsWithBkg(TH1* h, Option_t* optFit)
1566 {
1567  // fit with signal function removing bkg
1568 
1569  Int_t nPar = fBkgFunc->GetNpar();
1570  TF1* bkg =new TF1("Bkg", this, &Peak1D::Bkg, fBkgLeft1, fBkgRight2, nPar+2, "Peak1D", "Bkg");
1571  bkg->FixParameter(nPar, fBkgLeft2);
1572  bkg->FixParameter(nPar+1, fBkgRight1);
1573 
1574  // transfert initial parameter to bkg function
1575  for (Int_t i = 0; i < nPar; ++i)
1576  bkg->SetParameter(i, fBkgFunc->GetParameter(i));
1577 
1578  h->Fit(bkg, optFit);
1579 
1580  // set parameter back to background function
1581  fBkgFunc->SetParameters(bkg->GetParameters());
1582  fBkgLevelLeft2 = fBkgFunc->Eval(fBkgLeft2);
1583  fBkgLevelRight1 = fBkgFunc->Eval(fBkgRight1);
1584  delete bkg;
1585 
1586  // substract bkg and fit
1587  TH1* lochist = (TH1*)h->Clone();
1588  lochist->Add(fBkgFunc, -1);
1589  lochist->Fit(fSigFunc, optFit);
1590  delete lochist;
1591 }
1592 
1593 //__________________________________________________________
1594 void Peak1D::Print(Option_t* opt) const
1595 {
1596  //
1597  TString Opt(opt);
1598  if ( Opt.Contains("f") ) {
1599  if ( fSigFunc ) {
1600  cout << GetName() << " Fit results :"<<endl;
1601  for(int i=2 ; i<fSigFunc->GetNpar() ; i++)
1602  {
1603  if(((TString)fSigFunc->GetParName(i)) == "AmplitudeStep")
1604  continue;
1605 
1606  cout<<setw(10)<<fSigFunc->GetParName(i)<<": "<<setprecision(7)<<setw(10)<<fSigFunc->GetParameter(i)<<" ("<<setprecision(7)<<setw(10)<<2*fSigFunc->GetParError(i)<<")"<<endl;
1607  }
1608  cout<<setw(10)<<"Area"<<": "<<setprecision(7)<<setw(10)<<fPeakIntegral<<" ("<<setprecision(7)<<setw(10)<<fPeakIntegralError<<")"<<endl;
1609  if(fEfficiencyGraph != 0x0 || fEfficiencyFunc != 0x0)
1610  {
1611  Float_t Eval = 1;
1612 
1613  if(fEfficiencyFunc != 0x0)
1614  Eval = fEfficiencyFunc->Eval(fSigFunc->GetParameter("Position"));
1615  else
1616  Eval = fEfficiencyGraph->Eval(fSigFunc->GetParameter("Position"));
1617 
1618  Float_t CorrArea = fPeakIntegral/Eval;
1619  Float_t CorrAreaErr = CorrArea*sqrt(fPeakIntegralError/fPeakIntegral*fPeakIntegralError/fPeakIntegral+fEffRelError*fEffRelError);
1620  cout<<"Efficiency corrected Area "<<CorrArea<<" +- "<<CorrAreaErr<<endl;
1621 
1622  if(fRefArea != 1)
1623  {
1624  CorrArea = fPeakIntegral/Eval/fRefArea;
1625  CorrAreaErr = CorrArea*sqrt(fPeakIntegralError/fPeakIntegral*fPeakIntegralError/fPeakIntegral+fRefAreaError/fRefArea*fRefAreaError/fRefArea+fEffRelError*fEffRelError);
1626  cout<<"Efficiency corrected Area relative to ref"<<CorrArea<<" +- "<<CorrAreaErr<<endl;
1627  }
1628  }
1629  }
1630  }
1631  else {
1632  printf("\n");
1633  TNamed::Print(opt);
1634  BasePeak::Print(opt);
1635  printf("Position: %5.1f FWHM: %5.1f Resolution: %6.2f %%\n", fPosition, fFWHM, 100.*fFWHM/fPosition);
1636  if ( IsWithBkg() ) // set bkg
1637  printf("Bkg left1: %5.1f Bkg left2: %5.1f Bkg right1: %5.1f Bkg right1: %5.1f\n",
1638  fBkgLeft1, fBkgLeft2, fBkgRight1, fBkgRight2);
1639  }
1640 
1641  if(fIsAGammaSourceRay)
1642  {
1643  cout<<"Source Information: "<<fGammaSourceRay->Emmiter<<endl;
1644  cout<<"Energy = "<<fGammaSourceRay->Energy->GetValue()<< "("<<fGammaSourceRay->Energy->GetError()<<")"<<endl;
1645  cout<<"BranchingRatio = "<<fGammaSourceRay->Intensity->GetValue()<< "("<<fGammaSourceRay->Intensity->GetError()<<")"<<endl;
1646  }
1647 }
1648 
1649 
1650 //__________________________________________________________
1651 void Peak1D::SetLineColorPeak(Color_t color)
1652 {
1653  // set line color for peak
1654  fPolyLinePeak.SetLineColor(color);
1655 }
1656 
1657 //__________________________________________________________
1658 Color_t Peak1D::GetLineColorPeak()
1659 {
1660  return fPolyLinePeak.GetLineColor();
1661 }
1662 
1663 //__________________________________________________________
1664 void Peak1D::SetLineWidthPeak(Width_t width)
1665 {
1666  // set line width for peak
1667  fPolyLinePeak.SetLineWidth(width);
1668 }
1669 
1670 //__________________________________________________________
1671 void Peak1D::SetFillColorPeak(Color_t color)
1672 {
1673  // set fill color for peak
1674  fPolyLinePeak.SetFillColor(color);
1675 }
1676 
1677 //__________________________________________________________
1678 void Peak1D::SetFillStylePeak(Style_t style)
1679 {
1680  // set fill sryle for peak
1681  fPolyLinePeak.SetFillStyle(style);
1682 }
1683 
1684 
1685 //__________________________________________________________
1686 void Peak1D::SetLineColorBkg(Color_t color)
1687 {
1688  // set line color for peak
1689  fPolyLineBkgL.SetLineColor(color);
1690  fPolyLineBkgR.SetLineColor(color);
1691 }
1692 
1693 //__________________________________________________________
1694 void Peak1D::SetLineWidthBkg(Width_t width)
1695 {
1696  // set line width for peak
1697  fPolyLineBkgL.SetLineWidth(width);
1698  fPolyLineBkgR.SetLineWidth(width);
1699 }
1700 
1701 //__________________________________________________________
1702 void Peak1D::SetFillColorBkg(Color_t color)
1703 {
1704  // set fill color for peak
1705  fPolyLineBkgL.SetFillColor(color);
1706  fPolyLineBkgR.SetFillColor(color);
1707 }
1708 
1709 //__________________________________________________________
1710 void Peak1D::SetFillStyleBkg(Style_t style)
1711 {
1712  // set fill sryle for peak
1713  fPolyLineBkgL.SetFillStyle(style);
1714  fPolyLineBkgR.SetFillStyle(style);
1715 }
1716 
1717 //__________________________________________________________
1718 void Peak1D::SetDrawAs(EDrawAs d)
1719 {
1720  if ( fDrawAs == d )
1721  return;
1722  fDrawAs = d;
1723 
1724  // peak <==> gate
1725  SetModified(fPolyLinePeak);
1726 }
1727 
1728 //__________________________________________________________
1729 void Peak1D::Paint(Option_t* o)
1730 {
1731  TString opt = o;
1732 
1733  // first background
1734  if ( IsWithBkg() ) {
1735  if ( IsModified(fPolyLineBkgL) || IsModified(fPolyLineBkgR) )
1736  SetBkgPoints();
1737  if ( !opt.Contains("!l") ) {
1738  fPolyLineBkgL.Paint("f");
1739  fPolyLineBkgR.Paint("f");
1740  fPolyLineBkgL.Paint();
1741  fPolyLineBkgR.Paint();
1742  }
1743  }
1744 
1745  // paint peak but first recalculate the polyline
1746  if ( IsModified(fPolyLinePeak) )
1747  SetPeakPoints();
1748  if ( fDrawAs == kGate ) {
1749  fPolyLinePeak.Paint("f");
1750  }
1751  if ( !opt.Contains("!l") )
1752  fPolyLinePeak.Paint();
1753 
1754  if ( fDrawAs == kPeak ) {
1755  // if fit, functions
1756  if ( fPeakFunc && FitFlag() ) {
1757  fPeakFunc->Paint("same");
1758  if ( fBkgFunc && FitFlag() && IsWithBkg() )
1759  fBkgFunc->Paint("same");
1760  }
1761  // synchronize markers
1762  TMarker* m = 0x0;
1763  for (Int_t i = 0; i < fgkPointsPeak; ++i) {
1764  m = static_cast<TMarker*> (fMarkerList.At(i));
1765  m->SetX(fPolyLinePeak.GetX()[i]);
1766  m->SetY(fPolyLinePeak.GetY()[i]);
1767  }
1768  if ( !opt.Contains("!m") ) {
1769  fMarkerList.Paint();
1770  }
1771  }
1772 
1773  PeakCreator::HandleRefresh();
1774 }
1775 
1776 //__________________________________________________________
1777 void Peak1D::Draw(Option_t* opt)
1778 {
1779  // draw
1780  if (!gPad) {
1781  fLog.GetProcessMethod() = "Draw(Option_t* )"; fLog << error << "Pad not existing" << dolog;
1782  return;
1783  }
1784  TString tmp(opt);
1785  tmp.ToLower();
1786  if (tmp.Contains("g")) {
1787  SetDrawAs(kGate);
1788  }
1789  if (tmp.Contains("p")) {
1790  SetDrawAs(kPeak);
1791  }
1792 
1793  AppendPad(opt);
1794 }
1795 
1796 void Peak1D::SetDrawOption(Option_t* opt)
1797 {
1798  TString tmp(opt);
1799  tmp.ToLower();
1800  if (tmp.Contains("g")) {
1801  SetDrawAs(kGate);
1802  }
1803  if (tmp.Contains("p")) {
1804  SetDrawAs(kPeak);
1805  }
1806 
1807  TObject::SetDrawOption(opt);
1808 }
1809 
1810 
1811 
1812 //__________________________________________________________
1813 Int_t Peak1D::DistancetoPrimitive(Int_t px, Int_t py)
1814 {
1815  Int_t dP = fPolyLinePeak.DistancetoPrimitive(px, py);
1816  Int_t dL = fPolyLineBkgL.DistancetoPrimitive(px, py);
1817  Int_t dR = fPolyLineBkgR.DistancetoPrimitive(px, py);
1818 
1819  Int_t d1 = TMath::Min(dP, dL);
1820  Int_t d2 = TMath::Min(d1, dR);
1821  return d2;
1822 }
1823 
1824 //___________________________________________________________
1825 void Peak1D::ExecuteEvent(Int_t event, Int_t px, Int_t py)
1826 {
1827  if (!gPad->IsEditable())
1828  return;
1829 
1830  Int_t dP = 0;
1831  Int_t dL = 0;
1832  Int_t dR = 0;
1833 
1834  if (fCPolyline == 0x0) {
1835  dP = fPolyLinePeak.DistancetoPrimitive(px, py);
1836  dL = fPolyLineBkgL.DistancetoPrimitive(px, py);
1837  dR = fPolyLineBkgR.DistancetoPrimitive(px, py);
1838 
1839  if (dL <= dP) {
1840  if (dR <= dL) {
1841  fCPolyline = &fPolyLineBkgR;
1842  } else {
1843  fCPolyline = &fPolyLineBkgL;
1844  }
1845  } else {
1846  if (dR <= dP) {
1847  fCPolyline = &fPolyLineBkgR;
1848  } else {
1849  fCPolyline = &fPolyLinePeak;
1850  }
1851  }
1852  // keep copy of the polyline to be modified
1853  fCopyPolyline.SetPolyLine(fCPolyline->GetN(),fCPolyline->GetX(),fCPolyline->GetY());
1854  // SetModified((*fCPolyline),false);
1855  }
1856 
1857  switch (event) {
1858  case kButton1Down:
1859  case kMouseMotion:
1860  case kButton1Motion:
1861  case kButton1Up:
1862  fCPolyline->ExecuteEvent(event, px, py);
1863  break;
1864  case kButton1Double:
1865  // protection since in case of kButton1Double, down may not be called before up
1866  // in such case it crashed in polyline (static *x = 0x0 [root5.32] --> bug report to root --> Corrected svn 43961)
1867  fCPolyline->ExecuteEvent(kButton1Down, px, py);
1868  break;
1869  default:
1870  break;
1871  }
1872 
1873  Int_t how_many, which, sign; Double_t deltax, deltay, val;
1874  switch (event) {
1875 
1876  case kMouseMotion:
1877  fCPolyline = 0x0;
1878  break;
1879 
1880  case kButton1Up:
1881  // do your test before modifications !
1882  how_many = which = 0;
1883  for (Int_t i = 0; i < fCPolyline->GetN() ; i++ ) {
1884  if ( fCPolyline->GetX()[i] != fCopyPolyline.GetX()[i] ) {
1885  which = i;
1886  how_many++;
1887  }
1888  if ( fCPolyline->GetY()[i] != fCopyPolyline.GetY()[i] ) {
1889  which = i;
1890  how_many++;
1891  }
1892  // printf("%d %d %d : %f -> %f \n",i,which,how_many,fCPolyline->GetX()[i],fCopyPolyline.GetX()[i]);
1893  }
1894 
1895  if ( fCPolyline == &fPolyLinePeak ) {
1896  // the full peak has been moved, change postion & translate the polyline (to avoid re-calculations)
1897  if ( how_many > 2 ) {
1898  // just move on X
1899  deltax = fCPolyline->GetX()[fgkPeakIdx] - fCopyPolyline.GetX()[fgkPeakIdx];
1900  val = GetPosition() + deltax;
1901  // back the previous polyline before application of the shift
1902  fPolyLinePeak.SetPolyLine(fCopyPolyline.GetN(),fCopyPolyline.GetX(),fCopyPolyline.GetY());
1903  // check it does not go beyond background in case there is one
1904  if ( IsWithBkg() ) {
1905  if ( val > fBkgLeft2 && val < fBkgRight1 ) {
1906  fPosition = val;
1907  ShiftPolyline(fPolyLinePeak,deltax,0.0);
1908  SetModified(fPolyLinePeak,false);
1909  }
1910  }
1911  else {
1912  fPosition = val;
1913  ShiftPolyline(fPolyLinePeak,deltax,0.0);
1914  SetModified(fPolyLinePeak,false);
1915  }
1916  }
1917  else { // one of the point has moved only
1918  switch (which) {
1919  // change only the height
1920  case fgkPeakIdx:
1921  // compute modifications
1922  deltay = fCPolyline->GetY()[fgkPeakIdx] - fCopyPolyline.GetY()[fgkPeakIdx];
1923  val = GetHeight() + deltay;
1924  // apply modifications
1925  SetHeight(val);
1926  break;
1927  // width
1928  case fgkPeakIdx-1:
1929  case fgkPeakIdx+1:
1930  // compute modifications
1931  deltax = fCPolyline->GetX()[which] - fCopyPolyline.GetX()[which];
1932  if ( which == fgkPeakIdx-1 )
1933  sign = -1;
1934  else
1935  sign = +1;
1936 
1937  val = GetFWHM() + 2*sign*deltax;
1938  SetFWHM(val);
1939  break;
1940  default:
1941  break;
1942  }
1943  // modifications only of the peak values. polyline are delayed at paint
1944  fPolyLinePeak.SetPolyLine(fCopyPolyline.GetN(),fCopyPolyline.GetX(),fCopyPolyline.GetY());
1945  }
1946  // background
1947  } else {
1948  // the full bg has been moved, just translate the polyline
1949  if ( how_many > 2 ) {
1950  // just move on X
1951  deltax = fCPolyline->GetX()[0] - fCopyPolyline.GetX()[0];
1952  // check if move possible
1953  Bool_t do_move = false;
1954  if ( fCPolyline == &fPolyLineBkgL ) { //
1955  if ( fBkgLeft2 + deltax < fPolyLinePeak.GetX()[fgkPeakIdx-1] )
1956  do_move = true;
1957  }
1958  else {
1959  if ( fBkgRight1 + deltax > fPolyLinePeak.GetX()[fgkPeakIdx+1] )
1960  do_move = true;
1961  }
1962  // back the previous polyline before application of the shift
1963  fCPolyline->SetPolyLine(fCopyPolyline.GetN(),fCopyPolyline.GetX(),fCopyPolyline.GetY());
1964  // check it does not overlap with peak
1965  if ( do_move ) {
1966  // global shift
1967  ShiftPolyline((*fCPolyline),deltax,0.0);
1968  // compute new values of bg from polyline
1969  UpdateBkg();
1970  }
1971  }
1972  else { // one of the point has moved only
1973  switch (which) {
1974  // background level+position
1975  case 1:
1976  case 2:
1977  deltax = fCPolyline->GetX()[which] - fCopyPolyline.GetX()[which];
1978  // check it does not overlap with peak
1979  if ( fCPolyline->GetX()[which] < fPolyLinePeak.GetX()[fgkPeakIdx-1] || fCPolyline->GetX()[which] > fPolyLinePeak.GetX()[fgkPeakIdx+1] ) {
1980  // compute new values of bg from polyline
1981  UpdateBkg();
1982  SetPeakPoints();
1983  }
1984  // back the previous polyline before application of the shift
1985  else
1986  fCPolyline->SetPolyLine(fCopyPolyline.GetN(),fCopyPolyline.GetX(),fCopyPolyline.GetY());
1987  break;
1988  default:
1989  fCPolyline->SetPolyLine(fCopyPolyline.GetN(),fCopyPolyline.GetX(),fCopyPolyline.GetY());
1990  break;
1991  }
1992  }
1993  }
1994  fCPolyline = 0x0;
1995  break;
1996  }
1997 }
1998 
1999 void Peak1D::PaintFor(Double_t /*xmin*/, Double_t /*xmax*/, Double_t ymin, Double_t ymax)
2000 {
2001  if ( fDrawAs == kGate ) {
2002  fPolyLinePeak.GetY()[0] = fPolyLinePeak.GetY()[fgkPointsPeak - 1] = ymin;
2003 
2004  // change only if too low to be seen
2005  if ( fPolyLinePeak.GetY()[fgkPeakIdx] < 0.6*ymax ) {
2006  fPolyLinePeak.GetY()[fgkPeakIdx-1] = fPolyLinePeak.GetY()[fgkPeakIdx] = fPolyLinePeak.GetY()[fgkPeakIdx+1] = 0.85*ymax;
2007  }
2008  }
2009  fPolyLineBkgL.GetY()[0] = fPolyLineBkgL.GetY()[fgkPointsBkg-1] = ymin;
2010  fPolyLineBkgR.GetY()[0] = fPolyLineBkgR.GetY()[fgkPointsBkg-1] = ymin;
2011 }
2012 
2013 
2014 //___________________________________________
2015 Int_t Peak1D::Compare(const TObject *obj) const
2016 {
2017  // sort peak by position
2018  Peak1D* event = (Peak1D*) obj;
2019  return (fPosition > event->GetPosition()) ? 1 : -1;
2020 }
2021 
2022 //___________________________________________________________
2023 void Peak1D::SetPeakPoints(TPolyLine* polyline)
2024 {
2025 
2026  Double_t x[fgkPointsPeak] = {0,0,0,0,0,0,0};
2027  Double_t y[fgkPointsPeak] = {0,0,0,0,0,0,0};
2028 
2029  if (polyline->GetN() == fgkPointsPeak) {
2030  fPolyLinePeak.Copy((TPolyLine &)(*polyline));
2031 
2032  } else if (polyline->GetN() == fgkPointsPeak-2) {
2033  x[fgkPeakIdx] = polyline->GetX()[fgkPeakIdx-1];
2034  y[fgkPeakIdx] = polyline->GetY()[fgkPeakIdx-1];
2035  for (Int_t i = 1; i < 3; ++i) {
2036  x[fgkPeakIdx-i] = polyline->GetX()[fgkPeakIdx-1-i];
2037  y[fgkPeakIdx-i] = polyline->GetY()[fgkPeakIdx-1-i];
2038 
2039  x[fgkPeakIdx+i] = polyline->GetX()[fgkPeakIdx-1+i];
2040  y[fgkPeakIdx+i] = polyline->GetY()[fgkPeakIdx-1+i];
2041  }
2042  x[fgkPeakIdx+3] = x[fgkPeakIdx+2];
2043  y[fgkPeakIdx+3] = y[fgkPeakIdx+2];
2044  x[fgkPeakIdx-3] = x[fgkPeakIdx-2];
2045  y[fgkPeakIdx-3] = y[fgkPeakIdx-2];
2046  fPolyLinePeak.SetPolyLine(fgkPointsPeak, x, y);
2047 
2048  } else {
2049  fLog.SetProcessMethod("SetPeakPoints(TPolyLine* )");
2050  fLog << error << "Wrong number of points" << dolog;
2051  return;
2052  }
2053 
2054  fPolyLinePeak.SetLineWidth(2);
2055 }
2056 
2057 //___________________________________________________________
2058 
2059 void Peak1D::SetModified(TPolyLine &p, Bool_t modif)
2060 {
2061  p.SetBit(0x1000,modif); // in principle this bit is not used by Polyline
2062 }
2063 
2064 //___________________________________________________________
2065 Bool_t Peak1D::IsModified(TPolyLine &p)
2066 {
2067  return p.TestBit(0x1000);
2068 }
2069 
2070 //___________________________________________________________
2071 void Peak1D::SetPeakPoints()
2072 {
2073  // set points of polyline from peak members
2074  Int_t n = fgkPointsPeak;
2075 
2076  Double_t x[n], y[n];
2077  Double_t sigma = fFWHM/SigmaToFWHM; Double_t base = ( fBkgLevelLeft2 + fBkgLevelRight1 )/2.;
2078 
2079  // centroid
2080  x[fgkPeakIdx] = fPosition;
2081  y[fgkPeakIdx] = fHeight+base;
2082 
2083  // width
2084  x[fgkPeakIdx-1] = fPosition - fFWHM/2.;
2085  y[fgkPeakIdx-1] = fHeight/2 + base; // default model is a gaussian shape
2086 
2087  x[fgkPeakIdx+1] = fPosition + fFWHM/2.;
2088  y[fgkPeakIdx+1] = fHeight/2 + base; // default model is a gaussian shape
2089 
2090  if ( fDrawAs == kPeak ) {
2091 
2092  if ( FitFlag() && fSigFunc && fBkgFunc ) {
2093 
2094  // centroid
2095  x[fgkPeakIdx] = fPosition;
2096  y[fgkPeakIdx] = fHeight + fBkgFunc->Eval(x[fgkPeakIdx]);
2097 
2098  x[fgkPeakIdx-1] = fPosition - fFWHM/2.;
2099  y[fgkPeakIdx-1] = fSigFunc->Eval(x[fgkPeakIdx-1]) + fBkgFunc->Eval(x[fgkPeakIdx-1]);
2100 
2101  x[fgkPeakIdx+1] = fPosition + fFWHM/2.;
2102  y[fgkPeakIdx+1] = fSigFunc->Eval(x[fgkPeakIdx+1]) + fBkgFunc->Eval(x[fgkPeakIdx+1]);
2103 
2104  x[fgkPeakIdx-2] = fPosition - 2.5*sigma;
2105  x[fgkPeakIdx+2] = fPosition + 2.5*sigma;
2106  x[fgkPeakIdx-3] = fPosition - 3.5*sigma;
2107  x[fgkPeakIdx+3] = fPosition + 3.5*sigma;
2108 
2109  y[fgkPeakIdx-2] = fSigFunc->Eval(x[fgkPeakIdx-2]) + fBkgFunc->Eval(x[fgkPeakIdx-2]);
2110  y[fgkPeakIdx+2] = fSigFunc->Eval(x[fgkPeakIdx+2]) + fBkgFunc->Eval(x[fgkPeakIdx+2]);
2111  y[fgkPeakIdx-3] = fSigFunc->Eval(x[fgkPeakIdx-3]) + fBkgFunc->Eval(x[fgkPeakIdx-3]);
2112  y[fgkPeakIdx+3] = fSigFunc->Eval(x[fgkPeakIdx+3]) + fBkgFunc->Eval(x[fgkPeakIdx+3]);
2113  }
2114  else { // default model is a gaussian shape
2115 
2116  x[fgkPeakIdx-2] = fPosition - 2.5*sigma;
2117  x[fgkPeakIdx+2] = fPosition + 2.5*sigma;
2118  x[fgkPeakIdx-3] = fPosition - 3.5*sigma;
2119  x[fgkPeakIdx+3] = fPosition + 3.5*sigma;
2120 
2121  y[fgkPeakIdx-2] = fHeight*4.39e-02 + base;
2122  y[fgkPeakIdx+2] = fHeight*4.39e-02 + base;
2123  y[fgkPeakIdx-3] = fHeight*2.187e-03 + base;
2124  y[fgkPeakIdx+3] = fHeight*2.187e-03 + base;
2125  }
2126 
2127  } else if ( fDrawAs == kGate ) {
2128 
2129  x[fgkPeakIdx-3] = x[fgkPeakIdx-2] = x[fgkPeakIdx-1];
2130  x[fgkPeakIdx+3] = x[fgkPeakIdx+2] = x[fgkPeakIdx+1];
2131  y[fgkPeakIdx-1] = y[fgkPeakIdx+1] = y[fgkPeakIdx];
2132  y[fgkPeakIdx-2] = y[fgkPeakIdx+2] = y[fgkPeakIdx] - fHeight;
2133 
2134  } else {
2135  fLog.GetProcessMethod() = "SetPeakPoints()" ;
2136  fLog << error << "No Style for Drawing peak" << dolog;
2137  return;
2138  }
2139  //
2140  fPolyLinePeak.SetPolyLine(n, x, y);
2141  // the polyline is synchronised with it ==> set it to modify false
2142  SetModified(fPolyLinePeak,false);
2143 }
2144 
2145 //___________________________________________________________
2146 void Peak1D::SetBkgPoints()
2147 {
2148  // set points of polyline from peak members
2149  Int_t n = fgkPointsBkg;
2150  if ( !IsWithBkg() )
2151  return;
2152 
2153  Double_t xL[n], yL[n];
2154  Double_t xR[n], yR[n];
2155 
2156  // bkg left
2157  xL[0] = fBkgLeft1;
2158  xL[1] = fBkgLeft1;
2159  xL[2] = fBkgLeft2;
2160  xL[3] = fBkgLeft2;
2161 
2162  yL[0] = yL[3] = 0.;
2163  yL[1] = fBkgLevelLeft1;
2164  yL[2] = fBkgLevelLeft2;
2165 
2166  // bkg right
2167  xR[0] = fBkgRight1;
2168  xR[1] = fBkgRight1;
2169  xR[2] = fBkgRight2;
2170  xR[3] = fBkgRight2;
2171 
2172  yR[0] = yR[3] = 0.;
2173  yR[1] = fBkgLevelRight1;
2174  yR[2] = fBkgLevelRight2;
2175 
2176  fPolyLineBkgL.SetPolyLine(n, xL, yL);
2177  fPolyLineBkgR.SetPolyLine(n, xR, yR);
2178 
2179  // this line could be drawn as it they are consistent with the background definition
2180  SetModified(fPolyLineBkgL,false);
2181  SetModified(fPolyLineBkgR,false);
2182 }
2183 
2184 
2185 //___________________________________________________________
2186 void Peak1D::UpdatePeak()
2187 {
2188  // update peak members from polyline
2189 
2190  // centroid
2191  fPosition = fPolyLinePeak.GetX()[fgkPeakIdx];
2192 
2193  // width
2194  fFWHM = fPolyLinePeak.GetX()[fgkPeakIdx+1] - fPolyLinePeak.GetX()[fgkPeakIdx-1];
2195 
2196  // height
2197  Double_t base = 0;
2198 
2199  if (fBkgFlag) {
2200  base = (fPolyLineBkgL.GetY()[2] + fPolyLineBkgR.GetY()[1])/2.;
2201  }
2202 
2203  fHeight = fPolyLinePeak.GetY()[fgkPeakIdx] - base;
2204 
2205  if (fDrawAs == kPeak) {
2206  // area triangle + trapeze
2207  fIntensity = (2*fFWHM + fPolyLinePeak.GetX()[fgkPeakIdx+2]-fPolyLinePeak.GetX()[fgkPeakIdx-2])*fHeight/4.;
2208  }
2209 }
2210 
2211 //___________________________________________________________
2212 void Peak1D::UpdateBkg()
2213 {
2214  // update peak members from polyline
2215 
2216  fBkgLeft1 = fPolyLineBkgL.GetX()[0] = fPolyLineBkgL.GetX()[1];
2217  fBkgLeft2 = fPolyLineBkgL.GetX()[3] = fPolyLineBkgL.GetX()[2];
2218  fBkgRight1 = fPolyLineBkgR.GetX()[0] = fPolyLineBkgR.GetX()[1];
2219  fBkgRight2 = fPolyLineBkgR.GetX()[3] = fPolyLineBkgR.GetX()[2];
2220 
2221  fBkgLevelLeft1 = fPolyLineBkgL.GetY()[1];
2222  fBkgLevelLeft2 = fPolyLineBkgL.GetY()[2];
2223  fBkgLevelRight1 = fPolyLineBkgR.GetY()[1];
2224  fBkgLevelRight2 = fPolyLineBkgR.GetY()[2];
2225 
2226  // this line could be drawn as it they are consistent with the background definition
2227  SetModified(fPolyLineBkgL,false);
2228  SetModified(fPolyLineBkgR,false);
2229 }
2230 
2231 // private functions
2232 
2233 //___________________________________________________________
2234 void Peak1D::InitPeakMarkers()
2235 {
2236  fMarkerList.SetOwner(true);
2237  fMarkerList.Clear();
2238  for (Int_t i = 0; i < fgkPointsPeak ; i++) {
2239  TMarker *m = new TMarker();
2240  if ( i == fgkPeakIdx )
2241  m->SetMarkerStyle(23);
2242  else
2243  m->SetMarkerStyle(20);
2244  m->SetMarkerColor(fgMarkerColor);
2245 
2246  fMarkerList.Add(m);
2247  }
2248 }
2249 
2251 
2252 double Peak1D::Gaus(double*xx,double*pp)
2253 {
2254  double x = xx[0];
2255 
2256  int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
2257  int UsingBG = (int)pp[1]; //Using BG or not
2258 
2259  double f_tot = 0.;
2260 
2261  if(UsingBG != 0.) f_tot += fBkgFunc->Eval(x);
2262 
2263  Int_t Npar = 3;
2264 
2265  for(int i=0 ; i<NSubPeaks ; i++)
2266  {
2267  double Ampli = pp[2+i*Npar+0];
2268  double Mean = pp[2+i*Npar+1];
2269  double Sigma = pp[2+i*Npar+2]*1./sqrt(8.*log(2.));
2270 
2271  f_tot += Ampli*exp(-0.5*((x-Mean)/Sigma)*((x-Mean)/Sigma));
2272  }
2273 
2274  return f_tot;
2275 }
2276 
2278 
2279 double Peak1D::Landau(double*xx,double*pp)
2280 {
2281  double x = xx[0];
2282 
2283  int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
2284  int UsingBG = (int)pp[1]; //Using BG or not
2285 
2286  double f_tot = 0.;
2287 
2288  if(UsingBG != 0.) f_tot += fBkgFunc->Eval(x);
2289 
2290  Int_t Npar = 3;
2291 
2292  for(int i=0 ; i<NSubPeaks ; i++)
2293  {
2294  double Ampli = pp[2+i*Npar+0];
2295  double Mean = pp[2+i*Npar+1];
2296  double Sigma = pp[2+i*Npar+2]*1./sqrt(8.*log(2.));
2297 
2298  f_tot += 8*Ampli*TMath::Log(2)*TMath::Landau(x,Mean,Sigma);
2299  }
2300 
2301  return f_tot;
2302 }
2303 
2305 
2306 double Peak1D::DoubleSidedCrystalBallFct(double*xx,double*pp)
2307 {
2308  double x = xx[0];
2309 
2310  int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
2311  int UsingBG = (int)pp[1]; //Using BG or not
2312 
2313  double f_tot = 0.;
2314 
2315  if(UsingBG != 0.) f_tot += fBkgFunc->Eval(x);
2316 
2317  Int_t Npar = 7;
2318 
2319  for(int i=0 ; i<NSubPeaks ; i++)
2320  {
2321 
2322  // gaussian core
2323  double N = pp[2+i*Npar+0];//norm
2324  double mu = pp[2+i*Npar+1];//mean
2325  double sig = pp[2+i*Npar+2]*1./sqrt(8.*log(2.));//variance
2326  // transition parameters
2327  double a1 = pp[2+i*Npar+3];
2328  double p1 = pp[2+i*Npar+4];
2329  double a2 = pp[2+i*Npar+5];
2330  double p2 = pp[2+i*Npar+6];
2331 
2332  double u = (x-mu)/sig;
2333  double A1 = TMath::Power(p1/TMath::Abs(a1),p1)*TMath::Exp(-a1*a1/2);
2334  double A2 = TMath::Power(p2/TMath::Abs(a2),p2)*TMath::Exp(-a2*a2/2);
2335  double B1 = p1/TMath::Abs(a1) - TMath::Abs(a1);
2336  double B2 = p2/TMath::Abs(a2) - TMath::Abs(a2);
2337 
2338  double result(N);
2339  if (u<-a1) result *= A1*TMath::Power(B1-u,-p1);
2340  else if (u<a2) result *= TMath::Exp(-u*u/2);
2341  else result *= A2*TMath::Power(B2+u,-p2);
2342 
2343  f_tot += result;
2344  }
2345 
2346  return f_tot;
2347 }
2348 
2350 
2351 double Peak1D::DoubleTailedGaussian(double*xx,double*pp)
2352 {
2353  double x = xx[0];
2354 
2355  int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
2356  int UsingBG = (int)pp[1]; //Using BG or not
2357 
2358  double f_tot = 0.;
2359 
2360  if(UsingBG != 0.) f_tot += fBkgFunc->Eval(x);
2361 
2362  Int_t Npar = 5;
2363 
2364  for(int i=0 ; i<NSubPeaks ; i++)
2365  {
2366  double Ampli = pp[2+i*Npar+0];
2367  double Mean = pp[2+i*Npar+1];
2368  double Sigma = pp[2+i*Npar+2]*1./sqrt(8.*log(2.));
2369  double Lambda = pp[2+i*Npar+3];
2370  double Rho = pp[2+i*Npar+4];
2371  double S = 0.;//pp[2+i*Npar+5];
2372 
2373  double U = (x-Mean)/Sigma;
2374  double f_g = Ampli*TMath::Exp(-U*U*0.5);
2375  double f_lambda = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
2376  double f_rho = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
2377  double f_S = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));
2378 
2379  if(U<Lambda) f_tot += f_lambda;
2380  else if(U>Rho) f_tot += f_rho;
2381  else f_tot += f_g;
2382 
2383  f_tot += f_S;
2384  }
2385 
2386  return f_tot;
2387 }
2388 
2390 
2391 double Peak1D::DoubleTailedStepedGaussian(double*xx,double*pp)
2392 {
2393  double x = xx[0];
2394 
2395  int NSubPeaks = (int)pp[0]; //Number of subpeaks in the peak range
2396 
2397  Float_t f_tot = 0.;
2398 
2399  Int_t Npar = 5;
2400 
2401  Float_t Flat_BackGround = pp[1];
2402 
2403  f_tot += Flat_BackGround;
2404 
2405  for(int i=0 ; i<NSubPeaks ; i++)
2406  {
2407  Float_t Ampli = pp[2+i*Npar+0];
2408  Float_t Mean = pp[2+i*Npar+1];
2409  Float_t Sigma = pp[2+i*Npar+2]*1./sqrt(8.*log(2.));
2410  Float_t Lambda = pp[2+i*Npar+3];
2411  Float_t Rho = pp[2+i*Npar+4];
2412  Float_t S = pp[2+i*Npar+5];
2413 
2414  Float_t U = (x-Mean)/Sigma;
2415  Float_t f_g = Ampli*TMath::Exp(-U*U*0.5);
2416  Float_t f_lambda = Ampli*TMath::Exp(-0.5*Lambda*(2.*U-Lambda));
2417  Float_t f_rho = Ampli*TMath::Exp(-0.5*Rho*(2.*U-Rho));
2418  Float_t f_S = Ampli*S*1./((1+TMath::Exp(U))*(1+TMath::Exp(U)));
2419 
2420  if(U<Lambda) f_tot += f_lambda;
2421  else if(U>Rho) f_tot += f_rho;
2422  else f_tot += f_g;
2423 
2424  f_tot += f_S;
2425  }
2426 
2427  return f_tot;
2428 }
2429 
2430 TF1 *Peak1D::ExtractBkgFromDTGStep(TString Name)
2431 {
2432  TF1 *f = new TF1(Name,"[0]+[1]*[2]*1./((1+TMath::Exp((x-[3])/[4]))*(1+TMath::Exp((x-[3])/[4])))",fBkgLeft1, fBkgRight2);
2433  f->SetParameters(fSigFunc->GetParameter(1),
2434  fSigFunc->GetParameter(2),
2435  fSigFunc->GetParameter(7),
2436  fSigFunc->GetParameter(3),
2437  fSigFunc->GetParameter(4));
2438 
2439  f->SetBit(TObject::kCannotPick);
2440  gROOT->GetListOfFunctions()->Remove(f);
2441 
2442  return f;
2443 }
2444 
2445 ClassImp(Peak1D);
Gw::Peak1D::fPeakFunc
TF1 * fPeakFunc
Definition: Peak1D.h:371
Gw::Measure::GetError
virtual Data_T GetError() const
return the error on the measured value
Definition: Measure.h:109
Gw::Peak1D::ExecuteEvent
virtual void ExecuteEvent(Int_t event, Int_t px, Int_t py)
Handle event in pad.
Definition: Peak1D.cpp:1825
printf
printf("******************************************************************** \n")
PeakCreator.h
header file for PeakCreator.cpp
Gw::Peak1D::PaintFor
virtual void PaintFor(Double_t xmin, Double_t xmax, Double_t ymin, Double_t ymax)
set attribute to paint this gate on the given histogram
Definition: Peak1D.cpp:1999
Gw::Peak1D::FitMenu
void FitMenu(const char *nameFunc="gaus", Option_t *optFit="RN", Option_t *optBkg="pol1")
Fit with given function and background.
Definition: Peak1D.cpp:981
Gw::Peak1D::SetLineColorPeak
virtual void SetLineColorPeak(Color_t color)
Definition: Peak1D.cpp:1651
Gw::Peak1D::fLog
LogMessage fLog
Definition: Peak1D.h:397
Gw::Peak1D::fIsAGammaSourceRay
bool fIsAGammaSourceRay
Definition: Peak1D.h:388
Gw::Peak1D::SetPeak
void SetPeak(Double_t pos, Double_t height, Double_t fwhm, Double_t intensity=0)
Set peak members.
Definition: Peak1D.cpp:519
Gw::error
LogMessage & error(LogMessage &)
Definition: GwLogMessage.cpp:238
__LoadWatchers__::tmp
TString tmp
Definition: GANILLoadWatchers.C:47
Gw::warning
LogMessage & warning(LogMessage &)
Definition: GwLogMessage.cpp:232
Gw::Peak1D::Landau
double Landau(double *xx, double *pp)
Landau Fonction.
Definition: Peak1D.cpp:2279
PadManager.h
header file for PadManager.cpp
Gw::Peak1D::SignalAndBkgFunction
virtual TF1 * SignalAndBkgFunction()
Access to signal+Bkg function. Signa + Bkg must be set first, parameter are copied ! ...
Definition: Peak1D.cpp:483
Gw::Peak1D::WithBkg
void WithBkg(Bool_t with_bg=true)
toggle background (fit and display)
Definition: Peak1D.cpp:513
Gw::Peak1D::fCopyPolyline
TPolyLine fCopyPolyline
current polyline
Definition: Peak1D.h:358
Gw::Peak1D::SetPosition
virtual void SetPosition(const Double_t position, Option_t *axis="X")
Set Position of peak.
Definition: Peak1D.cpp:546
Gw::Peak1D::SetFWHM
virtual void SetFWHM(const Double_t FWHM, Option_t *axis="X")
Set FWHM of peak.
Definition: Peak1D.cpp:555
Gw::Peak1D::fgMarkerColor
static Color_t fgMarkerColor
Definition: Peak1D.h:383
Gw::Peak1D::Gaus
double Gaus(double *xx, double *pp)
Gaus Fonction.
Definition: Peak1D.cpp:2252
Gw::Peak1D::ExtractBkgFromDTGStep
TF1 * ExtractBkgFromDTGStep(TString Name)
Definition: Peak1D.cpp:2430
Gw::BasePeak::fIntensity
Double_t fIntensity
intensity of the peak
Definition: BasePeak.h:130
Gw::Peak1D::IsWithBkg
Bool_t IsWithBkg() const
Test if peak defined with background.
Definition: Peak1D.h:196
Gw::Peak1D::SetFillColorBkg
virtual void SetFillColorBkg(Color_t color)
Set line attribute for bkg.
Definition: Peak1D.cpp:1702
Gw::nline
LogMessage & nline(LogMessage &)
Definition: GwLogMessage.cpp:262
Gw::Peak1D::fgFillStylePeak
static Style_t fgFillStylePeak
Definition: Peak1D.h:380
Gw::Peak1D::fPolyLinePeak
TPolyLine fPolyLinePeak
Definition: Peak1D.h:353
Gw::Peak1D::SetFillStylePeak
virtual void SetFillStylePeak(Style_t style)
Definition: Peak1D.cpp:1678
Gw::PeakCreator::HandleRefresh
static void HandleRefresh(TCanvas *canvas=0x0)
Handle modification of the range of the pad.
Definition: PeakCreator.cpp:387
Gw::Peak1D::fgkPointsPeak
static const Int_t fgkPointsPeak
pointer to signal+ bkg function
Definition: Peak1D.h:375
Gw::BasePeak::gGateColor
static GateColor gGateColor
Definition: BasePeak.h:30
Gw::Peak1D::GetPeakIntegral
Double_t GetPeakIntegral() const
Get integral of peak after calling area method.
Definition: Peak1D.h:261
Gw::BasePeak::fDimension
UShort_t fDimension
dimension of the peak
Definition: BasePeak.h:131
Gw::BasePeak::GetHeight
virtual Double_t GetHeight() const
Get height of peak.
Definition: BasePeak.h:66
Gw::Peak1D::~Peak1D
virtual ~Peak1D()
Definition: Peak1D.cpp:442
Gw::Peak1D::SignalAndBkg
Double_t SignalAndBkg(Double_t *x, Double_t *par)
functions used to really fit
Definition: Peak1D.cpp:472
Gw::GammaSourceRay::Emmiter
TString Emmiter
Definition: NDB.h:46
Gw::BasePeak::EnableFit
virtual void EnableFit()
Set flag on to enable fit of the peak.
Definition: BasePeak.h:54
Gw::BasePeak::Print
virtual void Print(Option_t *opt="") const
print current peak
Definition: BasePeak.cpp:132
Gw::Peak1D::fDrawAs
EDrawAs fDrawAs
Definition: Peak1D.h:361
Gw::Peak1D::fPolyLineBkgL
TPolyLine fPolyLineBkgL
Definition: Peak1D.h:354
Gw::Peak1D::PeakFunction
virtual TF1 * PeakFunction()
Access to peak function ... don't delete it !
Definition: Peak1D.h:167
Gw::Peak1D::fGammaSourceRay
GammaSourceRay * fGammaSourceRay
Definition: Peak1D.h:389
Gw::Peak1D::Bkg
Double_t Bkg(Double_t *x, Double_t *par)
Definition: Peak1D.cpp:455
SigmaToFWHM
#define SigmaToFWHM
Definition: Peak1D.cpp:71
Gw::Peak1D::SetPeakIntegralError
void SetPeakIntegralError(Double_t IntErr)
Definition: Peak1D.h:263
Gw::Peak1D::GetLineColorPeak
virtual Color_t GetLineColorPeak()
Get line peak color.
Definition: Peak1D.cpp:1658
Gw::Peak1D::Fit
virtual void Fit(TH1 *histo, Option_t *optFit="RN", Option_t *optBkg="pol1")
Fit peak with background.
Definition: Peak1D.cpp:1406
Gw::Peak1D::fgkPointsBkg
static const Int_t fgkPointsBkg
Definition: Peak1D.h:376
Gw::Peak1D::GetFWHM
virtual Double_t GetFWHM(Option_t *axis="X") const
Get FWHM of peak.
Definition: Peak1D.cpp:962
Gw::Peak1D::fMarkerList
TList fMarkerList
current polyline before modification
Definition: Peak1D.h:360
Gw::dolog
LogMessage & dolog(LogMessage &)
Definition: GwLogMessage.cpp:250
Gw::BasePeak::fHeight
Double_t fHeight
height of the Peak
Definition: BasePeak.h:129
Gw::Peak1D::Paint
virtual void Paint(Option_t *opt="")
Paint method.
Definition: Peak1D.cpp:1729
Gw::Peak1D::SignalFunction
virtual TF1 * SignalFunction()
Access to signal function ... don't delete it !
Definition: Peak1D.h:164
Gw::Peak1D::GetPosition
virtual Double_t GetPosition(Option_t *axis="X") const
Get position of peak.
Definition: Peak1D.cpp:955
Gw::Peak1D::fBinWidth
Double_t fBinWidth
Definition: Peak1D.h:350
Gw::Peak1D::GetMarker
virtual TMarker * GetMarker(Int_t markerId)
Get marker for a given Id.
Definition: Peak1D.cpp:969
Gw::BasePeak::FitFlag
virtual Bool_t FitFlag() const
Get flag for enable fit.
Definition: BasePeak.h:75
Gw::Peak1D::Print
virtual void Print(Option_t *opt="f") const
Print current peak.
Definition: Peak1D.cpp:1594
Gw::BasePeak::SetIntensity
virtual void SetIntensity(const Double_t intensity)
Set intensity of peak.
Definition: BasePeak.h:51
Gw::Peak1D::BkgFunction
virtual TF1 * BkgFunction()
Access to signal function ... don't delete it !
Definition: Peak1D.h:174
Gw::Peak1D::SetHeight
void SetHeight(const Double_t height)
Set height of peak.
Definition: Peak1D.cpp:567
Gw::Peak1D::SetPeakFromFit
void SetPeakFromFit()
Set peak members + polyline.
Definition: Peak1D.cpp:530
Gw::GammaSourceRay::Energy
Measure< Double_t > * Energy
Definition: NDB.h:52
Gw::Peak1D::SetLineWidthPeak
virtual void SetLineWidthPeak(Width_t width)
Definition: Peak1D.cpp:1664
Gw::Peak1D::fgFuncColor
static Color_t fgFuncColor
Definition: Peak1D.h:384
Gw::Peak1D::fPolyLineBkgR
TPolyLine fPolyLineBkgR
Definition: Peak1D.h:355
Gw::Peak1D::DoubleTailedGaussian
double DoubleTailedGaussian(double *xx, double *pp)
Double tailed gaussian (Dino like)
Definition: Peak1D.cpp:2351
Gw::Peak1D::fgkPeakIdx
static const Int_t fgkPeakIdx
Definition: Peak1D.h:377
Gw::Peak1D::SetMarkerColor
void SetMarkerColor(Int_t color)
Set marker color of the peak polyline.
Definition: Peak1D.cpp:640
Gw::Peak1D::DoubleSidedCrystalBallFct
double DoubleSidedCrystalBallFct(double *xx, double *pp)
Double Sided Crystal Ball Fonction.
Definition: Peak1D.cpp:2306
Peak1D.h
header file for a general 1D peak
Gw::Peak1D::fSigFunc
TF1 * fSigFunc
Definition: Peak1D.h:369
Gw::Peak1D::SetBinWidth
virtual void SetBinWidth(Double_t width)
BinWidth of the fitted histogram for area correction.
Definition: Peak1D.h:200
endl
ADF::LogMessage & endl(ADF::LogMessage &log)
Definition: ADFLogMessage.cpp:265
Gw::Peak1D::Area
virtual void Area(Option_t *opt="h")
to get the Area (+bkg) of this peak
Definition: Peak1D.cpp:1031
Gw::Peak1D::SetFillStyleBkg
virtual void SetFillStyleBkg(Style_t style)
Definition: Peak1D.cpp:1710
Gw::BasePeak
A BasePeak is defined by a height, intensity and a dimension of the peak.
Definition: BasePeak.h:19
Gw::GammaSourceRay::Intensity
Measure< Double_t > * Intensity
Definition: NDB.h:53
Gw::Peak1D::DistancetoPrimitive
virtual Int_t DistancetoPrimitive(Int_t px, Int_t py)
Compute distance to object.
Definition: Peak1D.cpp:1813
Gw::Peak1D::SetLineColorBkg
virtual void SetLineColorBkg(Color_t color)
Definition: Peak1D.cpp:1686
Gw::Peak1D::Compare
virtual Int_t Compare(const TObject *obj) const
Definition: Peak1D.cpp:2015
Gw::Peak1D::fgLineColorBkg
static Color_t fgLineColorBkg
Definition: Peak1D.h:381
Gw::Peak1D
A graphical interface for placing schematic peak onto a 1D histogram with a given position...
Definition: Peak1D.h:79
Gw::Peak1D::SetIntensity
void SetIntensity(const Double_t intensity)
Set Intensity of peak.
Definition: Peak1D.cpp:579
Gw::Peak1D::FitCombined
virtual void FitCombined(TH1 *histo, TList *peakList, Option_t *optFit="RN", Option_t *optBkg="pol1")
Fit multi-defined peak with a common background.
Definition: Peak1D.cpp:1164
Gw::Peak1D::SetDrawAs
virtual void SetDrawAs(EDrawAs d)
Set drawing flag.
Definition: Peak1D.cpp:1718
Gw::Peak1D::fBkgFunc
TF1 * fBkgFunc
Definition: Peak1D.h:370
Gw::LogMessage::SetProcessMethod
virtual void SetProcessMethod(const char *)
To set the current method.
Definition: GwLogMessage.cpp:177
Gw::PadManager::GetHisto
static TH1 * GetHisto(TVirtualPad *pad=0x0, Option_t *op="")
look for an histogram into the pad
Definition: PadManager.cpp:67
Gw::Peak1D::fCPolyline
TPolyLine * fCPolyline
Definition: Peak1D.h:357
ClassImp
ClassImp(Peak1D)
Gw::BasePeak::SetHeight
virtual void SetHeight(const Double_t height)
Set height of peak.
Definition: BasePeak.h:48
Gw::Peak1D::SetLineWidthBkg
virtual void SetLineWidthBkg(Width_t width)
Definition: Peak1D.cpp:1694
Gw::Peak1D::IsPointInBkg
Short_t IsPointInBkg(Double_t x, Double_t=0)
to determine if a point is in bg. 0 likely in peak, 1 in bg, 2 outside bg
Definition: Peak1D.h:248
Gw::Peak1D::SetSignalFunction
virtual TF1 * SetSignalFunction(const char *nameFunc="gaus", TH1 *h=0x0)
Set pre-defined function to fit the signal.
Definition: Peak1D.cpp:716
Gw::LogMessage::GetProcessMethod
virtual std::string & GetProcessMethod()
To get the current method.
Definition: GwLogMessage.h:227
Gw::Peak1D::DoubleTailedStepedGaussian
double DoubleTailedStepedGaussian(double *xx, double *pp)
Double tailed gaussian with step (Dino like)
Definition: Peak1D.cpp:2391
Gw::Peak1D::fFitCombined
TF1 * fFitCombined
pointer to signal+ bkg function
Definition: Peak1D.h:372
Gw::Peak1D::Copy
virtual void Copy(TObject &o) const
Assignment operator.
Definition: Peak1D.cpp:353
Gw::Peak1D::SetBkgFunction
virtual TF1 * SetBkgFunction(const char *nameFunc="-")
Set pre-defined function for background during fit.
Definition: Peak1D.cpp:895
Gw::Peak1D::Draw
virtual void Draw(Option_t *opt="")
Definition: Peak1D.cpp:1777
Gw::Peak1D::SetFillColorPeak
virtual void SetFillColorPeak(Color_t color)
Set line attribute for peak.
Definition: Peak1D.cpp:1671
Gw::Peak1D::fgFillStyleBkg
static Style_t fgFillStyleBkg
Definition: Peak1D.h:382
Gw::Peak1D::SetDrawOption
virtual void SetDrawOption(Option_t *option="")
Definition: Peak1D.cpp:1796
Gw::Data::GetValue
Data_T GetValue() const
get the value, cannot be overloaded
Definition: Data.h:114
Gw::Peak1D::SetBackground
void SetBackground(Double_t bgLeft1, Double_t bgLeft2, Double_t bgRight1, Double_t bgRight2, Double_t bgLevelLeft1, Double_t bgLevelLeft2, Double_t bgLevelRight1, Double_t bgLevelRight2)
Set background limits.
Definition: Peak1D.cpp:591
Gw::BasePeak::Copy
virtual void Copy(TObject &o) const
Definition: BasePeak.h:42