Prototype2/hcalLab/peeks.C

0001 // Illustrates how to find peaks in histograms.
0002 // This script generates a random number of gaussian peaks
0003 // on top of a linear background.
0004 // The position of the peaks is found via TSpectrum and injected
0005 // as initial values of parameters to make a global fit.
0006 // The background is computed and drawn on top of the original histogram.
0007 //
0008 // To execute this example, do
0009 //  root > .x peaks.C  (generate 10 peaks by default)
0010 //  root > .x peaks.C++ (use the compiler)
0011 //  root > .x peaks.C++(30) (generates 30 peaks)
0012 //
0013 // To execute only the first part of the script (without fitting)
0014 // specify a negative value for the number of peaks, eg
0015 //  root > .x peaks.C(-20)
0016 //
0017 //Author: Rene Brun
0018
0019 #include "TCanvas.h"
0020 #include "TMath.h"
0021 #include "TH1.h"
0022 #include "TF1.h"
0023 #include "TRandom.h"
0024 #include "TSpectrum.h"
0025 #include "TVirtualFitter.h"
0026
0027 Int_t npeaks = 30;
0028 Double_t fpeaks(Double_t *x, Double_t *par) {
0029    Double_t result = par[0] + par[1]*x[0];
0030    for (Int_t p=0;p<npeaks;p++) {
0031       Double_t norm  = par[3*p+2];
0032       Double_t mean  = par[3*p+3];
0033       Double_t sigma = par[3*p+4];
0034       result += norm*TMath::Gaus(x[0],mean,sigma);
0035    }
0036    return result;
0037 }
0038 void peaks(Int_t np=10) {
0039    npeaks = TMath::Abs(np);
0040    TH1F *h = new TH1F("h","test",500,0,1000);
0041    //generate n peaks at random
0042    Double_t par[3000];
0043    par[0] = 0.8;
0044    par[1] = -0.6/1000;
0045    Int_t p;
0046    for (p=0;p<npeaks;p++) {
0047       par[3*p+2] = 1;
0048       par[3*p+3] = 10+gRandom->Rndm()*980;
0049       par[3*p+4] = 3+2*gRandom->Rndm();
0050    }
0051    TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
0052    f->SetNpx(1000);
0053    f->SetParameters(par);
0054    TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
0055    c1->Divide(1,2);
0056    c1->cd(1);
0057    h->FillRandom("f",200000);
0058    h->Draw();
0059    TH1F *h2 = (TH1F*)h->Clone("h2");
0060    //Use TSpectrum to find the peak candidates
0061    TSpectrum *s = new TSpectrum(2*npeaks);
0062    Int_t nfound = s->Search(h,2,"",0.10);
0063    printf("Found %d candidate peaks to fit\n",nfound);
0064    //Estimate background using TSpectrum::Background
0065    TH1 *hb = s->Background(h,20,"same");
0066    if (hb) c1->Update();
0067    if (np <0) return;
0068
0069    //estimate linear background using a fitting method
0070    c1->cd(2);
0071    TF1 *fline = new TF1("fline","pol1",0,1000);
0072    h->Fit("fline","qn");
0073    //Loop on all found peaks. Eliminate peaks at the background level
0074    par[0] = fline->GetParameter(0);
0075    par[1] = fline->GetParameter(1);
0076    npeaks = 0;
0077    Double_t *xpeaks = s->GetPositionX();
0078    for (p=0;p<nfound;p++) {
0079       Double_t xp = xpeaks[p];
0080       Int_t bin = h->GetXaxis()->FindBin(xp);
0081       Double_t yp = h->GetBinContent(bin);
0082       if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
0083       par[3*npeaks+2] = yp;
0084       par[3*npeaks+3] = xp;
0085       par[3*npeaks+4] = 3;
0086       npeaks++;
0087    }
0088    printf("Found %d useful peaks to fit\n",npeaks);
0089    printf("Now fitting: Be patient\n");
0090    TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
0091    //we may have more than the default 25 parameters
0092    TVirtualFitter::Fitter(h2,10+3*npeaks);
0093    fit->SetParameters(par);
0094    fit->SetNpx(1000);
0095    h2->Fit("fit");
0096 }
0097
0098
0099
0100
0101
0102 // ;; This buffer is for notes you don't want to save, and for Lisp evaluation.
0103 // ;; If you want to create a file, visit that file with C-x C-f,
0104 // ;; then enter the text in that file's own buffer.
0105
0106 // **************************************************************************
0107 //  assume background is linear what is most likely not true ....
0108 //  par[0]   background constant
0109 //     [1]   background slope
0110 //  assume the signal is the sequence of equally spaced gaussian peaks
0111 //  par[2]   peak to peak separation
0112 //     [3]   normalization constant for the first peak (pedestal)
0113 //     [4]   location of the first peak
0114 //     [5]   rms of the first peak - at this stage we will allow it to vary peak-to-peak
0115 //     then 2 parameters per peak (normalization and RMS)
0116 Double_t sipmpeaks(Double_t *x, Double_t *par) {
0117    Double_t result  = par[0] + par[1]*x[0];
0118    //  add pedestal peak first
0119    result          += par[3]*TMath::Gaus(x[0],par[4],par[5]);
0120    for (Int_t p=0; p<npeaks; p++) {
0121       Double_t norm  = par[6  + p*2];
0122       Double_t mean  = par[4] + par[2]*(p+1);
0123       Double_t sigma = par[8  + p*2];
0124       result += norm*TMath::Gaus(x[0],mean,sigma);
0125    }
0126    return result;
0127 }
0128 /// **************************************************************************
0129 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
0130 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
0131 Int_t spmcalib(TString hName)
0132 {
0133   TH1F * data = (TH1F*)(gSystem->FindObject(hName));
0134   if(!data){
0135     cout<<"spmCalib  Data Histogramm "<<hName<<"  not found"<<endl;
0136     exit();
0137   }
0138   Int_t      maxPeaks(20);        //  maximum number of peaks we have chance to resolve
0139   Double_t   range(150.);         //  range of amplitudes for peak search
0140   Double_t   Double_t norm0(1.),  slope0(-1./range);
0141   Double_t   peak_to_peak(10.);
0142   Double_t   pLoc(0.),            pRMS(1.5);
0143   Double_t   par[6+maxPeaks*2];   //  backgrownd, pedestal and nPeaks of gaussian peaks with constant peak-to-peak separation
0144
0145   par[0] = norm0;
0146   par[1] = slope0;
0147   par[2] = peak_to_peak;       //  peak-to-peak separation
0148   par[3] = norm0;              //  normalization for pedestal peak
0149   par[4] = pLoc;               //  pedestal mean
0150   par[5] = pRMS;               //  pedestal rms
0151   for (int p=0; p< maxPeaks; p++) {
0152       par[6 + p*3] = norm0;                       //  norm of the peak p
0153       par[7 + p*3] = par[4] + par[2]*(p+1);       //  mean of the peak p
0154       par[8 + p*3] = 2.;                          //  rms  of the peak p
0155    }
0156   TString scId   = "sc_"; scId += hName;
0157   TF1  * scF     = (TF1*) (gROOT->FindObject(scId));               if(scF) delete scF;
0158   scF            = new TF1(scId,&sipmpeaks,0,range,6+maxPeaks*2);
0159   TH1F * dataCl  = (TH1F*)(data->Clone("dataCl"));
0160   dataCl        -> SetAxisRange(0., range);
0161   Int_t  binsToUse  = dataCl->GetXaxis()->FindBin(range);
0162   TCanvas * peaks = new TCanvas("peaks","peaks",10,10,1000,900);
0163   peaks         -> Divide(1,2);
0164   peaks         -> cd(1);
0165   dataCl        -> Draw();
0166   TH1F * dataCl2 = (TH1F*)data->Clone("dataCl2");
0167
0168   //  Use TSpectrum to find the peak candidates
0169   TSpectrum *s = new TSpectrum(maxPeaks+1);
0170   Int_t nfound = s->Search(dataCl, 2, "", 0.10);
0171   printf("TSpectrum:  Found %d candidate peaks to fit\n",nfound);
0172
0173   //Estimate background using TSpectrum::Background
0174   TH1 *hb = s->Background(h,20,"same");
0175   if (hb)   peaks->Update();
0176   //  if (np <0) return;
0177   //  function to estimate linear background
0178   peaks           -> cd(2);
0179   TF1 *fline       = new TF1("fline","pol1",0,range);
0180   dataCl          -> Fit("fline","qn");
0181   //Loop on all found peaks. Eliminate peaks at the background level
0182   par[0]           = fline->GetParameter(0);
0183   par[1]           = fline->GetParameter(1);
0184   Int_t     npeaks = 0;
0185   Double_t *xpeaks = s->GetPositionX();
0186   Double_t spacing = 0.;
0187   for (p=0; p<nfound; p++) {
0188     Double_t xp    = xpeaks[p];
0189     Int_t    bin   = dataCl->GetXaxis()->FindBin(xp);
0190     Double_t yp    = dataCl->GetBinContent(bin);
0191     cout<<"Testing Peak "<<p<<"  xp "<<xp<<"  bin "<<bin<<"  yP "<<yp<<"  BCKG "<<fline->Eval(xp)<<endl;
0192     if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
0193     if(npeaks) spacing += xp - par[4 + (npeaks-1)+3];
0194     par[3 + npeaks*3] = yp;
0195     par[4 + npeaks*3] = xp;
0196     par[5 + npeaks*3] = 3.;
0197     npeaks++;
0198   }
0199   printf("Found %d useful peaks to fit\n",npeaks);
0200   printf("Now fitting: Be patient\n");
0201    TF1 *fit = new TF1("fit",&sipmpeaks,0,1000,3+3*npeaks);  // first peak is the pedestal
0202   //we may have more than the default 25 parameters
0203   TVirtualFitter::Fitter(h2,10+3*npeaks);
0204   fit->SetParameters(par);
0205   fit->SetNpx(1000);
0206   h2->Fit("fit");
0207 }
0208
0209
0210
0211
0212
0213 /// **************************************************************************
0214 // do single cell calibration using low end data stored in rv_# and fm_# channel histograms
0215 // the assumption - we have beautiful picture with at least 20 peaks !!!!!!
0216 void sccalib(nPeaks = 20)
0217 {
0218   Double_t   rcpar[2+3+2*nPeaks];   //  backgrownd and nPeaks of uncorrelated gaussian peaks
0219   Double_t   fcpar[2+3+2*nPeaks];   //  backgrownd and nPeaks of uncorrelated gaussian peaks
0220
0221   int nx_c = 2;
0222   int ny_c = NCH;
0223   TString cSC = "Single Cell Calibration (black - raw, red - fit) Run "; cSC += runId.Data();
0224   scCalibDisplay = new TCanvas("scCalibDisplay",cSC,200*nx_c,100*ny_c);
0225   scCalibDisplay->Divide(nx_c, ny_c);
0226   gStyle->SetOptStat(1110);
0227   gStyle->SetOptFit(111);
0228   Int_t  binsToUse(200);     Double_t slope0(-1./(Double_t))binsToUse);  Double_t norm0(1.);
0229   Double_t peak_to_peak(10.);
0230   for (int ich = 0; ich<NCH; ich++){
0231     TH1F *hrv = (TH1F*)rpeak[ich];
0232     hrv       -> SetLineColor(1);
0233     hrv       -> SetAxisRange(0., binsToUse);
0234     TH1F *fmv = (TH1F*)fm[ich];
0235     fmv       -> SetLineColor(2);
0236     fmv       -> SetAxisRange(0., binsToUse);
0237
0238     //  estimate exponential backgrownd we go to 500 ADC counts
0239     rcpar[0] = norm0;
0240     rcpar[1] = slope0;
0241     rcpar[2] = peak_to_peak;       //  peak-to-peak separation
0242     rcpar[3] = norm0;              //  normalization for pedestal peak
0243     rcpar[4] = 2.;                 //  pedestal mean
0244     rcpar[5] = 1.5;                //  pedestal rms
0245     //  we fit peak-to-peak not the individual independent peak positions
0246     for (int p=1; p<= nPeaks; p++) {
0247       rcpar[2*p+4] = norm0;    //  norm of the peak p
0248       rcpar[2*p+5] = 2.;       //  rms of the peak p
0249    }
0250     fcpar[0] = 1.;
0251     fcpar[1] = slope0;
0252     fcpar[2] = peak_to_peak;       //  peak-to-peak separation
0253     fcpar[3] = norm0;              //  normalization for pedestal peak
0254     fcpar[4] = 2.;                 //  pedestal mean
0255     fcpar[5] = 1.5;                //  pedestal rms
0256     //  we fit peak-to-peak not the individual independent peak positions
0257     for (int p=1; p<= nPeaks; p++) {
0258       fcpar[2*p+4] = norm0;    //  norm of the peak p
0259       fcpar[2*p+5] = 2.;    //  rms of the peak p
0260    }
0261
0262     TString rvscId = "rvsc_"; rvscId += ich;
0263     TString fmscId = "fmsc_"; fmscId += ich;
0264     rvsc[ich] = (TF1*) (gROOT->FindObject(rvscId));           if(rvsc[ich]) delete rvsc[ich];
0265     fmsc[ich] = (TF1*) (gROOT->FindObject(fmscId));           if(fmsc[ich]) delete fmsc[ich];
0266     rvsc[ich] = new TF1(rvscId,&scpeaks,0,500,5+2*nPeaks);    fmsc[ich] = new TF1(fmscId,&scpeaks,0,500,5+2*nPeaks);
0267     rvsc[ich] -> SetNpx(binsToUse);                           fmsc[ich] -> SetNpx(binsToUse);
0268     rvsc[ich] -> SetParameters(rcpar);                        fmsc[ich] -> SetParameters(fcpar);
0269     rvsc[ich] -> SetLineColor(1);                             fmsc[ich] -> SetLineColor(2);
0270     TH1F * rh = (TH1F*)(hrv->Clone("rh"));                    TH1F * fh = (TH1F*)(fmv->Clone("fh"));
0271
0272     //  function to estimate linear background
0273     TF1 *fline = new TF1("fline","pol1",0,binsToUse);
0274     //  Use TSpectrum to find the peak candidates
0275     TSpectrum *s = new TSpectrum(2*nPeaks);
0276     scCalibDisplay -> cd(ich*2+1);                            hrv -> Draw();
0277     Int_t      rfound = s->Search(hrv,1,"new");
0278     printf("Found %d candidate peaks in Raw hist ",rfound);
0279     scCalibDisplay -> Update();
0280     hrv -> Fit("fline","qn");
0281     rcpar[0] = fline->GetParameter(0);
0282     rcpar[1] = fline->GetParameter(1);
0283     //Loop on all found peaks. Eliminate peaks at the background level
0284     nrpeaks = 0;
0285     Float_t *xrpeaks = s->GetPositionX();
0286     for (p=0;p<nfound;p++) {
0287       Float_t xp  = xpeaks[p];
0288       Int_t   bin = hrv->GetXaxis()->FindBin(xp);
0289       Float_t yp  = h->GetBinContent(bin);
0290       if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
0291       par[2*npeaks+2] = yp;
0292       par[3*npeaks+3] = xp;
0293       par[3*npeaks+4] = 3;
0294       nrpeaks++;
0295     }
0296
0297
0298
0299
0300     scCalibDisplay -> cd(ich*2+2);                            fmv -> Draw();
0301     int_t ffound = s->Search(fmv,1,"new");
0302     printf("Found %d candidate peaks in Fit hist ",ffound);
0303     scCalibDisplay -> Update();
0304     fmv -> Fit("fline","qn");
0305     fcpar[0] = fline->GetParameter(0);
0306     fcpar[1] = fline->GetParameter(1);
0307
0308     //Loop on all found peaks. Eliminate peaks at the background level
0309    npeaks = 0;
0310    Float_t *xpeaks = s->GetPositionX();
0311    for (p=0;p<nfound;p++) {
0312       Float_t xp = xpeaks[p];
0313       Int_t bin = h->GetXaxis()->FindBin(xp);
0314       Float_t yp = h->GetBinContent(bin);
0315       if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue;
0316       par[3*npeaks+2] = yp;
0317       par[3*npeaks+3] = xp;
0318       par[3*npeaks+4] = 3;
0319       npeaks++;
0320    }
0321
0322
0323
0324
0325
0326
0327
0328
0329     //    rvsc[ich] -> SetParNames ("#mu","g","#sigma","A");         fmsc[ich] -> SetParNames ("#mu","g","#sigma","A");
0330     scCalibDisplay   -> cd(ich+1); hrv->SetAxisRange(0., 200); hrv ->Draw();  fmv->Draw("same");
0331     hrv       -> Fit(rvsc[ich],"rnQ");                         fmv -> Fit(fmsc[ich],"rnq");
0332     rvsc[ich] -> Draw("same");                                 fmsc[ich] -> Draw("same");
0333   }
0334   scCalibDisplay -> Update();
0335 }