]>
Commit | Line | Data |
---|---|---|
8215f572 | 1 | /* Author: Anders G. Knospe, The University of Texas at Austin |
2 | Created: 28 January 2014 | |
3 | ||
4 | This macro is a simple implementation of a Voigtian peak fit atop a background. A Voigtian peak is a convolution of a Breit-Wigner peak and a Gaussian. Here, it is implemented using the TMath::Voigt() function in ROOT. | |
5 | ||
6 | Input Parameters: | |
7 | h: input histogram containing a peak | |
8 | formula: character string containing the functional form of the background, as it would be implemented in a TF1. | |
9 | For example, a quadratic background can be implemented as "pol2(0)" or "[0]+[1]*x+[2]*x*x". | |
10 | par: array containing peak parameters {mass,resolution,width} | |
11 | fix: array containing three flags {fix mass,fix resolution,fix width}, which allow any of these peak parameters to be fixed during the fit | |
12 | range: array containing four values {R1,R2,R3,R4} | |
13 | R1 (R4) is the lower (upper) limit of the fit range | |
14 | R2 (R3) is the lower (upper) limit of the peak range. The peak range is excluded from the fit initially so that the background can be estimated. The peak range is included in all subsequent fits. | |
15 | file: file to which results may be saved (optional) | |
16 | */ | |
17 | ||
18 | int VoigtianFit(TH1D* h,char* formula,double* par,int* fix,double* range,TFile* file) | |
19 | { | |
20 | if(!h){cerr<<"Error in VoigtianFit(): missing histogram"<<endl; return 1;} | |
21 | ||
22 | char name[500]; sprintf(name,h->GetName()); | |
23 | cout<<"using VoigtianFit():\n h="<<name<<"\n background formula="<<formula<<endl; | |
24 | cout<<" mass="<<par[0]; | |
25 | if(fix[0]) cout<<" (fixed)"<<endl; | |
26 | else cout<<" (free)"<<endl; | |
27 | cout<<" resolution="<<par[1]; | |
28 | if(fix[1]) cout<<" (fixed)"<<endl; | |
29 | else cout<<" (free)"<<endl; | |
30 | cout<<" width="<<par[2]; | |
31 | if(fix[2]) cout<<" (fixed)"<<endl; | |
32 | else cout<<" (free)"<<endl; | |
33 | cout<<" range="<<range[0]<<" "<<range[1]<<" "<<range[2]<<" "<<range[3]<<endl; | |
34 | if(file) cout<<" file="<<file->GetName()<<endl; | |
35 | else cout<<" output not saved"<<endl; | |
36 | ||
37 | int j; | |
38 | ||
39 | //create copy of histogram h with peak removed | |
40 | TH1D* a=(TH1D*) h->Clone(Form("%s_nopeak",name)); | |
41 | for(j=h->GetXaxis()->FindBin(1.000001*range[1]);j<=h->GetXaxis()->FindBin(0.999999*range[2]);j++){ | |
42 | a->SetBinContent(j,0.); | |
43 | a->SetBinError(j,0.); | |
44 | } | |
45 | ||
46 | //get initial estimate of background | |
47 | TF1* fb=new TF1(Form("%s_back",name),formula,range[0],range[3]); | |
48 | a->Fit(fb,"RQN"); | |
49 | ||
50 | //define peak fit function | |
51 | int vp=fb->GetNpar(); | |
52 | TF1* fp=new TF1(Form("%s_peak",name),Form("%s+[%i]*TMath::Voigt(x-[%i],[%i],[%i])",formula,vp,vp+1,vp+2,vp+3),range[0],range[3]); | |
53 | ||
54 | //set initial parameter values, only the peak height is free | |
55 | for(j=0;j<vp;j++){fp->SetParameter(j,fb->GetParameter(j)); fp->FixParameter(j,fb->GetParameter(j));} | |
56 | fp->SetParameter(vp,h->GetBinContent(h->GetXaxis()->FindBin(0.5*(range[2]-range[1])))-fb->Eval(0.5*(range[2]-range[1]))); | |
57 | for(j=0;j<3;j++){fp->SetParameter(vp+j+1,par[j]); fp->FixParameter(vp+j+1,par[j]);} | |
58 | ||
59 | h->Fit(fp,"RQN"); | |
60 | ||
61 | if(!fix[2]){//release width | |
62 | fp->ReleaseParameter(vp+3); | |
63 | fp->SetParError(vp+3,0.1*fp->GetParameter(vp+3)); | |
64 | h->Fit(fp,"RQN"); | |
65 | } | |
66 | ||
67 | if(!fix[1]){//release resolution | |
68 | fp->ReleaseParameter(vp+2); | |
69 | fp->SetParError(vp+2,0.1*fp->GetParameter(vp+2)); | |
70 | h->Fit(fp,"RQN"); | |
71 | } | |
72 | ||
73 | if(!fix[0]){//release mass | |
74 | fp->ReleaseParameter(vp+1); | |
75 | fp->SetParError(vp+1,0.1*fp->GetParameter(vp+1)); | |
76 | h->Fit(fp,"RQN"); | |
77 | } | |
78 | ||
79 | //release background constant parameter | |
80 | fp->ReleaseParameter(0); | |
81 | fp->SetParError(0,fb->GetParError(0)); | |
82 | h->Fit(fp,"RQN"); | |
83 | ||
84 | //release other background parameters | |
85 | for(j=1;j<vp;j++){ | |
86 | fp->ReleaseParameter(j); | |
87 | fp->SetParError(j,fb->GetParError(j)); | |
88 | } | |
89 | h->Fit(fp,"RQN"); | |
90 | ||
91 | //final fit | |
92 | cerr<<"doing final fit"<<endl; | |
93 | h->Fit(fp,"RNI"); | |
94 | ||
95 | //save output | |
96 | if(file){ | |
97 | file->cd(); | |
98 | h->Write(); | |
99 | a->Write(); | |
100 | fb->Write(); | |
101 | fp->Write(); | |
102 | } | |
103 | ||
104 | return 0; | |
105 | } |