#include "TSpectrum" #include "TSpectrumFit" #include "TH1.h" #include "TF1.h" #include "TVirtualFitter.h" #include #include #include Int_t npeaks = 30; Float_t factor=1.e-6*1.e-26/1e4;//factor of 10^-6 due to micro Jy, factor of 1e-26 to conversion to J s-1 cm-2 Hz-1// in order to have cm^-2 its 10^4 cm^2 Float_t c=2.998e8; //speed of light in vaccuum in m/s reftractive index of air would be 1.000277 Float_t planck=6.62606957e-34; //planck constant to calculate enrgy of a photon with given frequency Float_t FoV=12600.;//MAGIC FoV of 3.5 deg expressed in arsec Float_t area=236e4;//Area of one MAGIC telescope in cm^2 //generate n peaks at random Double_t par[3000]; par[0] = 0.8; par[1] = -0.6/1000; Double_t fpeaks(Double_t *x, Double_t *par) { Double_t result = par[0] + par[1]*x[0]; for (Int_t p=0;pSetOptStat(0); Double_t x,y; TNtupleD* data = new TNtupleD("data", "data", "x:y"); data->ReadFile(file); // if (file==NULL){ // printf("File %s not found\n",data); // break; // } TCanvas *ntup=new TCanvas("ntuple", "ntuple"); cout<<"Read "<SetBranchAddress("x",&x); data->SetBranchAddress("y",&y); data->Draw("y:x","","l"); TF1 *fbg=new TF1("fbg","1.42766e-2*x",400,1000); Int_t n = data->GetEntries(); Float_t nu; Double_t e; cout<<"Spectrum contains "<GetEntry(0); Float_t xmin=x; data->GetEntry(n-1); Float_t xmax=x; TGraph *g=(TGraph*)ntup->GetListOfPrimitives()->FindObject("Graph"); TH2F *temp=(TH2F*)ntup->GetListOfPrimitives()->FindObject("htemp"); TH1F *h=new TH1F(); TH1F *h1= new TH1F("h1","h1",n,htemp->GetXaxis()->GetXmin(),htemp->GetXaxis()->GetXmax()); TH1F *h2= new TH1F("h2","h2",n,htemp->GetXaxis()->GetXmin(),htemp->GetXaxis()->GetXmax()); for(Int_t i=0;iGetEvent(i); nu=c/(x*1.e-9); e=nu*planck; h->Fill(x,y); h1->Fill(x,y/e*factor); h2->Fill(x,y/e*factor*area*FoV**2); } TCanvas *c4=new TCanvas(); h->Draw("same l"); break; // c1->Update(); TCanvas *c=new TCanvas("c","c",1000,1000); c->Divide(1,4); c->cd(1); h->GetYaxis()->SetTitle("Intensity [#muJy arcsec^{-2}] "); h->GetXaxis()->SetTitle("#lambda [nm]"); h->Draw(); c->cd(2); h1->GetYaxis()->SetTitle("Intensity [photons s^{-1} cm^{-2} arcsec^{-2}]"); h1->GetXaxis()->SetTitle("#lambda [nm]"); h1->Draw(); c->cd(3); h2->GetYaxis()->SetTitle("Intensity [photons s^{-1}]"); h2->GetXaxis()->SetTitle("#lambda [nm]"); h2->Draw(); c->cd(4); npeaks = TMath::Abs(np); TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks); f->SetNpx(1000); TH1F *h3 = (TH1F*)h->Clone("h3"); //Use TSpectrum to find the peak candidates TSpectrum *s = new TSpectrum(2*npeaks); Int_t nfound = s->Search(h,1,"",0.000001); printf("Found %d candidate peaks to fit\n",nfound); //Estimate background using TSpectrum::Background TH1 *hb = s->Background(h,30,"same"); if (hb) c->Update(); if (np <0) return; //estimate linear background using a fitting method TF1 *fline = new TF1("fline","pol1",0,1000); h3->Fit("fline","qn"); //Loop on all found peaks. Eliminate peaks at the background level par[0] = fline->GetParameter(0); par[1] = fline->GetParameter(1); npeaks = 0; Float_t *xpeaks = s->GetPositionX(); for (Int_t p=0;pGetXaxis()->FindBin(xp); Float_t yp = h3->GetBinContent(bin); if (yp-TMath::Sqrt(yp) < fline->Eval(xp)) continue; par[3*npeaks+2] = yp; par[3*npeaks+3] = xp; par[3*npeaks+4] = 3; npeaks++; } printf("Found %d useful peaks to fit\n",npeaks); printf("Now fitting: Be patient\n"); TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks); //we may have more than the default 25 parameters TVirtualFitter::Fitter(h3,10+3*npeaks); fit->SetParameters(par); fit->SetNpx(1000); h3->Fit("fit"); c->cd(); //c->SaveAs("NSB_Fit.root"); //c->SaveAs("NSB_Fit.C"); }