//Programma di analisi dei file .root prodotti dai TLG //#include #include #include "TCanvas.h" #include "TH1.h" #include "TF1.h" #include "TSpectrum.h" #include "TVirtualFitter.h" Double_t Background(Double_t *x, Double_t *par); Double_t Signal(Double_t *x, Double_t *par); Double_t Total(Double_t *x, Double_t *par); Int_t npeaks; //<=========== //main(char *name) void peak_finder() { //gROOT->Reset(); //Int_t npeaks; npeaks=5; //TFile *file = new TFile(name); TFile *file = new TFile("peaks_histo.root"); MyWinSpectrumPeaks = new TCanvas("MyWinSpectrumPeaks", "spectrum peaks", 200, 10, 900, 700); //TH1F *h_slope = new TH1F("h_slope","Peaks slopes",100,-0.7,0.7); TH1F *h_slope = (TH1F*)file->Get("h_slope;1"); MyWinSpectrumPeaks->Divide(1,2); MyWinSpectrumPeaks->cd(1); h_slope->Draw(); //Use TSpectrum to find the peak candidates TSpectrum *spc = new TSpectrum(2*npeaks); //Int_t nfound = spc->Search(h_slope,1,"new"); Int_t nfound = spc->Search(h_slope,2,""); printf("Found %d candidate peaks to fit\n",nfound); MyWinSpectrumPeaks->Update(); MyWinSpectrumPeaks->cd(2); //signal Double_t *par; par = new Double_t[3+3*nfound]; npeaks = nfound; //Get peaks position and set initial parameters: Float_t *xpeaks = spc->GetPositionX(); for(Int_t i=0; iGetXaxis()->FindBin(xp); printf("Global bin number corresponding to %f is: %d \n", xp,bin); Float_t yp = h_slope->GetBinContent(bin); printf("Content of bin number %f is: %d \n", xp,yp); par[3*p+3] = yp; par[3*p+4] = xp; //par[3*p+5] = 3; par[3*p+5] = .01; printf("Fit parameters: %f %f %f \n", par[3*p+3],par[3*p+4],par[3*p+5]); } //Fit background using a gaussian TF1 *fbackg = new TF1("fbackg","gaus",-0.7,0.7); //h_slope->Fit("fbackg","q0"); h_slope->Fit("fbackg"); par[0] = fbackg->GetParameter(0); par[1] = fbackg->GetParameter(1); par[2] = fbackg->GetParameter(2); printf("/////////////////////////////////////////////////// \n"); printf("Background fit parameters: %f %f %f \n", par[0],par[1],par[2]); printf("/////////////////////////////////////////////////// \n"); //Subtract estimated background to signal //Creates a temporary histogram TH1F *h_temp = (TH1F*)h_slope->Clone(); h_temp->Reset(); TF1 *eback = h_slope->GetFunction("fbackg"); for(Int_t bin = 1;bin <= 100; bin++) { Float_t x = h_slope->GetBinCenter(bin); Double_t fval = eback->Eval(x); Double_t diff = TMath::Abs(fval - h_slope->GetBinContent(bin)); h_temp->Fill(x,diff); } //Fit signal with my function TF1 *fsignal = new TF1("fsignal",Signal,-0.7,0.7,3*npeaks); h_temp->GetListOfFunctions()->Add(fsignal); h_temp->Fit("fsignal"); //TF1 *esig = h_temp->GetFunction("gaus"); //TF1 *esig = h_temp->GetFunction("fsignal"); //Fit background + signal printf("Now be fitting, be patient \n"); eback->GetParameters(&par[0]); //esig->GetParameters(&par[3]); fsignal->GetParameters(&par[3]); for(Int_t m=0; m<3+3*nfound; m++) { printf("Check parameters: %f \n", par[m]); } //Make the final function for fit: TF1 *theory = new TF1("theory",Total,-0.7,0.7,3+3*nfound); theory->SetParameters(par); //TVirtualFitter::Fitter(h_slope,10+3*npeaks); //we may have more than the default 25 parameters //TVirtualFitter::Fitter(h_slope)->SetFCN("theory"); //theory->SetParameters(par); //theory->SetNpx(1000); //gStyle->SetOptFit(1111); //h_slope->GetListOfFunctions()->Add(theory); h_slope->Fit("theory"); } Double_t Background(Double_t *x, Double_t *par) { //The background function Double_t arg = 0; if(par[2]) arg = (x[0] - par[1])/par[2]; Double_t val = par[0]*TMath::Exp(-0.5*arg*arg)*x[0]*x[0]; return val; } Double_t Signal(Double_t *x, Double_t *par) { //The signal function Double_t result = 0; Double_t e = 0; //Int_t npeaks = 6; for(Int_t p=0;p