Prototype2/hcalLab/fitter.C

0001 #include "TCanvas.h"
0002 #include "TH1.h"
0003 #include "TF1.h"
0004 #include "TRandom.h"
0005 #include "TSpectrum.h"
0006 #include "TVirtualFitter.h"
0007
0008 Int_t npeaks = 30;
0009 Double_t fpeaks(Double_t *x, Double_t *par) {
0010    Double_t result = par[0] + par[1]*x[0];
0011    for (Int_t p=0;p<npeaks;p++) {
0012       Double_t norm  = par[3*p+2];
0013       Double_t mean  = par[3*p+3];
0014       Double_t sigma = par[3*p+4];
0015       result += norm*TMath::Gaus(x[0],mean,sigma);
0016    }
0017    return result;
0018 }
0019 void peaks(Int_t np=10) {
0020    npeaks = np;
0021    //  create histogram with x[0,1000] with randomly positioned np peaks each sigma of 4
0022    //  total number of parameters to fit are 2(linear backgrownd) +3*np
0023    TH1F *h = new TH1F("h","test",500,0,1000);
0024    //generate n peaks at random
0025    Double_t par[3000];
0026    par[0] = 0.8;
0027    par[1] = -0.6/1000;
0028    Int_t p;
0029    for (p=0;p<npeaks;p++) {
0030       par[3*p+2] = 1;
0031       par[3*p+3] = 10+gRandom->Rndm()*980;
0032       par[3*p+4] = 3+2*gRandom->Rndm();
0033    }
0034    TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
0035    f->SetNpx(1000);
0036    f->SetParameters(par);
0037    TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
0038    c1->Divide(1,2);
0039    c1->cd(1);
0040    h->FillRandom("f",200000);
0041    h->Draw();
0042    //  now extracting those pareameters from MC generated histogram
0043    TH1F *h2 = (TH1F*)h->Clone("h2");
0044    //Use TSpectrum to find the peak candidates
0045    TSpectrum *s = new TSpectrum(2*npeaks);
0046    Int_t nfound = s->Search(h,1,"new");
0047    printf("Found %d candidate peaks to fitn",nfound);
0048    c1->Update();
0049    c1->cd(2);
0050
0051    //estimate linear background
0052    TF1 *fline = new TF1("fline","pol1",0,1000);
0053    h->Fit("fline","qn");
0054    //Loop on all found peaks. Eliminate peaks at the background level
0055    par[0] = fline->GetParameter(0);
0056    par[1] = fline->GetParameter(1);
0057    npeaks = 0;
0058    Float_t *xpeaks = s->GetPositionX();
0059    for (p=0;p<nfound;p++) {
0060       Float_t xp = xpeaks[p];
0061       Int_t bin = h->GetXaxis()->FindBin(xp);
0062       Float_t yp = h->GetBinContent(bin);
0063       if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
0064       par[3*npeaks+2] = yp;
0065       par[3*npeaks+3] = xp;
0066       par[3*npeaks+4] = 3;
0067       npeaks++;
0068    }
0069
0070    printf("Found %d useful peaks to fitn",npeaks);
0071    printf("Now fitting: Be patientn");
0072    TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
0073    TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
0074    fit->SetParameters(par);
0075    fit->SetNpx(1000);
0076    h2->Fit("fit");
0077 }
0078