C++ TSpectrum::Background方法代码示例

void Background_order() {
   Int_t i;
   const Int_t nbins = 4096;
   Double_t xmin     = 0;
   Double_t xmax     = 4096;
   Double_t source[nbins];
   gROOT->ForceStyle();

   TH1F *d1 = new TH1F("d1","",nbins,xmin,xmax);
   TH1F *d2 = new TH1F("d2","",nbins,xmin,xmax);
   TH1F *d3 = new TH1F("d3","",nbins,xmin,xmax);
   TH1F *d4 = new TH1F("d4","",nbins,xmin,xmax);

   TString dir  = gROOT->GetTutorialsDir();
   TString file = dir+"/spectrum/TSpectrum.root";
   TFile *f     = new TFile(file.Data());
   TH1F *back = (TH1F*) f->Get("back2");
   back->SetTitle("Influence of clipping filter difference order on the estimated background");
   back->SetAxisRange(1220,1460);
   back->SetMaximum(3000);
   back->Draw("L");

   TSpectrum *s = new TSpectrum();

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder2,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d1->SetBinContent(i + 1,source[i]);
   d1->SetLineColor(kRed);
   d1->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder4,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d2->SetBinContent(i + 1,source[i]);
   d2->SetLineColor(kBlue);
   d2->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder6,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d3->SetBinContent(i + 1,source[i]);
   d3->SetLineColor(kGreen);
   d3->Draw("SAME L");

   for (i = 0; i < nbins; i++) source[i]=back->GetBinContent(i + 1);
   s->Background(source,nbins,40,TSpectrum::kBackDecreasingWindow,
                 TSpectrum::kBackOrder8,kFALSE,
                 TSpectrum::kBackSmoothing3,kFALSE);
   for (i = 0; i < nbins; i++) d4->SetBinContent(i + 1,source[i]);
   d4->SetLineColor(kMagenta);
   d4->Draw("SAME L");
}

void peaks(Int_t np=10) {
   npeaks = TMath::Abs(np);
   TH1F *h = new TH1F("h","test",500,0,1000);
   //generate n peaks at random
   Double_t par[3000];
   par[0] = 0.8;
   par[1] = -0.6/1000;
   Int_t p;
   for (p=0;p<npeaks;p++) {
      par[3*p+2] = 1;
      par[3*p+3] = 10+gRandom->Rndm()*980;
      par[3*p+4] = 3+2*gRandom->Rndm();
   }
   TF1 *f = new TF1("f",fpeaks,0,1000,2+3*npeaks);
   f->SetNpx(1000);
   f->SetParameters(par);
   TCanvas *c1 = new TCanvas("c1","c1",10,10,1000,900);
   c1->Divide(1,2);
   c1->cd(1);
   h->FillRandom("f",200000);
   h->Draw();
   TH1F *h2 = (TH1F*)h->Clone("h2");
   //Use TSpectrum to find the peak candidates
   TSpectrum *s = new TSpectrum(2*npeaks);
   Int_t nfound = s->Search(h,2,"",0.10);
   printf("Found %d candidate peaks to fit\n",nfound);
   //Estimate background using TSpectrum::Background
   TH1 *hb = s->Background(h,20,"same");
   if (hb) c1->Update();
   if (np <0) return;

   //estimate linear background using a fitting method
   c1->cd(2);
   TF1 *fline = new TF1("fline","pol1",0,1000);
   h->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;
   Double_t *xpeaks = s->GetPositionX();
   for (p=0;p<nfound;p++) {
      Double_t xp = xpeaks[p];
      Int_t bin = h->GetXaxis()->FindBin(xp);
      Double_t yp = h->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(h2,10+3*npeaks);
   fit->SetParameters(par);
   fit->SetNpx(1000);
   h2->Fit("fit");
}

void FindPeak(TH1 *hm, int * i, char * namefile){
  int np =5, p, max = 0;
  int npeaks = TMath::Abs(np);
  double par[3000];
  par[0] = 0.8;
  par[1] = -0.6/1000;

  for (p=0;p<npeaks;p++) {
    par[3*p+2] = 1;
    par[3*p+3] = 10+gRandom->Rndm()*980;
    par[3*p+4] = 3+2*gRandom->Rndm();
  }
 
  TSpectrum *s = new TSpectrum(2*npeaks,1);
  int nfound = s->Search(hm,2,"",0.10);
  printf("Found %d candidate peaks to fit\n",nfound);

  TH1 *hb = s->Background(hm,20,"same");
  if (np <0) return;

  // loope over peaks
  TF1 *fline = new TF1("fline","pol1",0,1000);
  hm->Fit("fline","qn");
  par[0] = fline->GetParameter(0);
  par[1] = fline->GetParameter(1);
  npeaks = 0;
  float *xpeaks = s->GetPositionX();
  for (p=0;p<nfound;p++) {
    float xp = xpeaks[p];
    int bin = hm->GetXaxis()->FindBin(xp);
    float yp = hm->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");
  if (max < npeaks) max = npeaks;
  TF1 *fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
  TVirtualFitter::Fitter(hm,10+3*npeaks);
  fit->SetParameters(par);
  fit->SetNpx(1000);
  hm->Fit("fit"); 
}

////    This version needs thought
void SiCalibrator::FindPeaks() {
	Load();
	int iterate = 25;
	int nbins;
	TSpectrum* spec = new TSpectrum(10);
	avesigma = 0;
	int nsigma = 0;
	for (int src=0; src < CalData.size(); src++) {
		if (CalData[src].hSource != 0) {
			for (int ch=0; ch<CalData[src].sourcedata.size(); ch++) {
				TH1D* hbi = CalData[src].hSource->ProjectionY("hb",ch+1,ch+1);
				nbins = hbi->GetNbinsX();
				TH1D* hbk1 = (TH1D*) spec->Background(hbi,iterate);
				hbi->Add(hbk1,-1);
				//Estimate parameters
				double lim = 5;
				double max = 0;
				double peak = 0;
				int i = nbins - 1;
				while (hbi->GetBinContent(i) < lim && i > 0) {
					max = i;
					i--;
				}
				double sigma = (max/CalibSource::sourcelist[src].betas.back().E)*2; //2 keV
				if (sigma > 1) {
					avesigma += sigma;
					nsigma++;
				}
				gErrorIgnoreLevel = kError;
				int npeaks = spec->Search(hbi,sigma,"nodraw ",0.001);
				float* adc = spec->GetPositionX();
				float* amp = spec->GetPositionY();
				for (int i=0; i<npeaks; i++) {
					CalData[src].sourcedata[ch].ADC.push_back(adc[i]);
					CalData[src].sourcedata[ch].Amp.push_back(amp[i]);
				}
			}
		}
	}
	avesigma = avesigma/nsigma;
}

//________________________________________________________________________________
void peaks(TH1* h, Int_t np=10, Int_t ij=0) {
  if (! h) return;
  npeaks = TMath::Abs(np);
  if (! c1)  c1 = new TCanvas();
  else       c1->Clear();
  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
  c1->Divide(1,2);
  c1->cd(1);
  h->Draw();
  TH1F *h2 = (TH1F*)h->Clone("h2");
  //Use TSpectrum to find the peak candidates
  TSpectrum *s = new TSpectrum(2*npeaks);
  Int_t nfound = s->Search(h,5,"",0.05);
  printf("Found %d candidate peaks to fit\n",nfound);
  if (! nfound) return;
  //Estimate background using TSpectrum::Background
  TH1 *hb = s->Background(h,20,"same");
  hb->Draw("same");
  c1->Update();
  if (c2 && ij > 0) {c2->cd(ij); h->Draw(); c2->Update();}
  if (np <0) return;

  //estimate linear background using a fitting method
  c1->cd(2);
  TF1 *fline = new TF1("fline","pol1",0,1000);
  fline->FixParameter(1,0.);
  h->Fit("fline","qnlw");
  if (c2 && ij > 0) {c2->cd(ij+1); h->Draw(); c2->Update(); c1->cd(2);}
  //Loop on all found peaks. Eliminate peaks at the background level
  Double_t par[3000];
  par[0] = fline->GetParameter(0);
  par[1] = fline->GetParameter(1);
  npeaks = 0;
  Float_t *xpeaks = s->GetPositionX();
  Float_t ymax = 0;
  for (Int_t p=0;p<nfound;p++) {
    Float_t xp = xpeaks[p];
    Int_t bin = h->GetXaxis()->FindBin(xp);
    Float_t yp = h->GetBinContent(bin);
    if (yp-3*TMath::Sqrt(yp) < fline->Eval(xp)) continue;
    par[3*npeaks+2] = yp;
    if (yp > ymax) ymax = yp;
    par[3*npeaks+3] = xp;
    par[3*npeaks+4] = 3;
    npeaks++;
  }
  cout << "Found " << npeaks << " useful peaks to fit" << endl;
  Int_t NP = 0;
  Int_t nbad = 0;
  TString LineH("  {\"");  LineH +=  h->GetName(); LineH += "\"";
  TString Line("");
  struct ParErr_t {Double_t par, err;};
  ParErr_t parErr[10];
  TF1 *fit = 0;
  if (ymax > 2*par[0]) {
    cout << "Now fitting: Be patient" << endl;
    fit = new TF1("fit",fpeaks,0,1000,2+3*npeaks);
    TVirtualFitter::Fitter(h2,10+3*npeaks); //we may have more than the default 25 parameters
    fit->SetParameter(0,par[0]);
    fit->FixParameter(1,0.);
    for (Int_t p = 0; p < npeaks; p++) {
      fit->SetParName(3*p+2,Form("A%i",p));
      fit->SetParLimits(3*p+2,0,1e6);
      fit->SetParameter(3*p+2,par[3*p+2]);
      fit->SetParName(3*p+3,Form("#mu%i",p));
      fit->SetParLimits(3*p+3,TMath::Max(0.,par[3*p+3]-2), TMath::Min(240.,par[3*p+3]+2));
      fit->SetParameter(3*p+3,par[3*p+3]);
      fit->SetParName(3*p+4,Form("#sigma%i",p));
      fit->SetParLimits(3*p+4,0.01,20);
      fit->SetParameter(3*p+4,par[3*p+4]);
    }
    fit->SetNpx(1000);
    h2->SetStats(0);
    h2->Fit("fit","l");
    if (c2 && ij > 0) {c2->cd(ij); h2->Draw("same"); c2->Update(); c1->cd(2);}
    fit->GetParameters(par);
    for (Int_t p = 0; p<np;p++) {
      if (p < npeaks && par[3*p+2] > 5*fit->GetParError(3*p+2) && 
	  par[3*p+2] > par[0]) {
	if (TMath::Abs(par[3*p+4]) > 5.0) nbad++;
	//      Line += Form(",%f,%f,%7.2f,%5.2f",par[3*p+2],fit->GetParError(3*p+2),par[3*p+3],TMath::Abs(par[3*p+4]));
	parErr[NP].par = par[3*p+3];
	parErr[NP].err = TMath::Abs(par[3*p+4]);
	for (Int_t i = 0; i < NP; i++) {
	  if (parErr[i].par > parErr[NP].par) {
	    ParErr_t temp = parErr[i];
	    parErr[i] = parErr[NP];
	    parErr[NP] = temp;
	  }
	}
	NP++;
      } 
    }
  }
  for (Int_t p = 0; p < np; p++) {
    if (p < NP) Line += Form(",%7.2f,%5.2f",parErr[p].par,parErr[p].err);
    else  Line += ",0,0";
  }
  Line += "},";
//.........这里部分代码省略.........

//.........这里部分代码省略.........
			/**************************************
		     * Peaks estimation
		     **************************************
		    */
			for (int i = 0; i < nbins; ++i){
				g->SetPoint(i, i, caen.trace[channel][i]);
			}

			Double_t y_max = TMath::MaxElement(g->GetN(),g->GetY());
			Float_t * source = new Float_t[nbins];
			Float_t * dest = new Float_t[nbins];

			for (int i = 0; i < nbins; ++i){
				source[i]=y_max-caen.trace[channel][i];
				g->SetPoint(i, i, source[i]);
			}

			//Use TSpectrum to find the peak candidates
			TSpectrum *s = new TSpectrum();
		   	Int_t nfound = s->SearchHighRes(source, dest, nbins, 3, 2, kTRUE, 2, kFALSE, 5);

		    /**************************************
		     * Background estimation
		     **************************************
		    */
		    Int_t  	ssize = nbins;
		    Int_t  	numberIterations = 20;
		    Int_t  	direction = s->kBackIncreasingWindow;
		    Int_t  	filterOrder = s->kBackOrder2;
		    bool  	smoothing = kFALSE;
		    Int_t  	smoothWindow = s->kBackSmoothing3;
		    bool  	compton = kFALSE;
		    for (int i = 0; i < nbins; ++i) baseline[i] = source[i];
		    s->Background(baseline, ssize, numberIterations, direction, filterOrder, smoothing, smoothWindow, compton);

		    /**************************************
		     * Peaks and integral estimation
		     **************************************
		    */
		    Double_t px[2], py[2];
		    for (int i = 0; i < nbins; ++i) dest[i] = source[i]-baseline[i];
		    if(nfound==2){
			   bin = s->GetPositionX()[0];
			   fPositionX1 = bin;
			   fPositionY1 = dest[bin];
			   px[0] = bin;
			   py[0] = dest[bin];
			   bin = s->GetPositionX()[1];
			   fPositionX2 = bin;
			   fPositionY2 = dest[bin];
			   px[1] = bin;
			   py[1] = dest[bin];
		    }
			int posxa=6;
		    int posxb=9;
		    switch (filenum){
		    	case 1081:
		    		posxa=6;
		    		posxb=9;
		    		break;

		    	case 1082:
		    		posxa=5;
		    		posxb=7;
		    		break;

void Play()
{

    //gSystem->Load("libRooFit");
    //using namespace RooFit;

    TChain* S = new TChain("mjdTree");
    //Look at COPPI calibration data from Nov 2013
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001923.root");
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001924.root");
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001925.root");
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001926.root");
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001927.root");
    S->Add("/project/projectdirs/majorana/data/mjd/surfprot/data/gatified/P3AKF/mjd_run40001928.root");
    S->Print();
    //S->Add("COPPIsA__run1131.root");
    
    TH1D *hc146 = new TH1D("hc146","Channel 146", 2000, 0, 800E3 );
    TH1D *hc147 = new TH1D("hc147","Channel 147", 2000, 0, 800E3 );

   TCanvas* c1 = new TCanvas("c1","",1500,800);
   TCanvas* c2 = new TCanvas("c2","",1500,800);
   TCanvas* c3 = new TCanvas("c3","",1500,800);
   c1->Divide(1,2);
   c2->Divide(1,2);
   c3->Divide(1,2);
    
   // --- Find Peaks in Spectrum ---
   c1->cd(1);
   S->Draw("energy>>hc146","channel==146");
   c1->cd(2);
   S->Draw("energy>>hc147","channel==147");

      // we calibrate with 133Ba and 60Co so let's build TSpectrum and look for those peaks
      // let's focus on the highest intensity lines:
      // 30.973 0.628
      // 356.0129 0.6205
      // 30.625 0.34
      // 80.9979 0.329
      //
      // 1332.492 0.999826
      // 1173.228 0.9985
      // 8.26-8.33 0.0000136
      // 7.46-7.48 0.000098

   //Use a gaussian to fit each of the peaks.
   TF1* tgaus = new TF1("tgaus","gaus",0,900E3);
   
   c1->cd(1);
   TSpectrum *s = new TSpectrum(12);
   Int_t nfound = s->Search(hc146,3,"",0.05);
   vector<calData> x146; //TSpectrum guess
   printf("Found %d candidate peaks to fit in channel 146 spectrum:\n",nfound);
   TH1 *hb = s->Background(hc146,20,"same");
   TH1D* hc146bf = (TH1D*)hc146->Clone("hc146bf");
   hc146bf->Add(hb,-1);
   if (hb) c1->Update();
   c2->cd(1);
   hc146bf->Draw();;
   Float_t *xpeaks = s->GetPositionX();
   calData d;
   for (int i = 0; i < nfound; i++) {
      //printf("%f : %f \n",s->GetPositionX()[i],s->GetPositionY()[i]);
      d.adc=s->GetPositionX()[i];
      x146.push_back(d);
   }
   sort(x146.begin(),x146.end(),CompareByadc);
   for(std::vector<calData>::iterator it=x146.begin(); it!=x146.end(); ++it)
   {
      tgaus->SetParameter(1,(*it).adc);
      TFitResultPtr r = hc146bf->Fit(tgaus,"SQ+","",(*it).adc-0.02*((*it).adc),(*it).adc+0.02*((*it).adc));
      (*it).fadc=r->Parameter(1);
      (*it).efadc=r->ParError(1);
   }
   cout << " Ts X \t\t Fit X \t\t\t err(x) " << endl;
   for(int i = 0; i < nfound; i++)
   {
      printf("%f \t %f \t +/- \t %f \n",x146[i].adc,x146[i].fadc,x146[i].efadc);
   }
   
   c1->cd(2);
   nfound = s->Search(hc147,3,"",0.05);
   vector<calData> x147; //TSpectrum guess
   printf("Found %d candidate peaks to fit in channel 147 spectrum:\n",nfound);
   TH1 *hb147 = s->Background(hc147,20,"same");
   TH1D* hc147bf = (TH1D*)hc147->Clone("hc147bf");
   hc147bf->Add(hb147,-1);
   if (hb147) c1->Update();
   c2->cd(2);
   hc147bf->Draw();;
   xpeaks = s->GetPositionX();
   for (int i = 0; i < nfound; i++) {
      //printf("%f : %f \n",s->GetPositionX()[i],s->GetPositionY()[i]);
      d.adc=s->GetPositionX()[i];
      x147.push_back(d);
   }
   sort(x147.begin(),x147.end(),CompareByadc);
   for(std::vector<calData>::iterator it=x147.begin(); it!=x147.end(); ++it)
   {
      tgaus->SetParameter(1,(*it).adc);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	TFile *infile = new TFile(filename);
	TTree *tree= (TTree*)infile->Get("acq_tree_0");
	TBranch *branch = tree->GetBranch(Form("acq_ch%d",channel));
	branch->SetAddress(&inc_data.timetag);	

	TH1F *h_raw = new TH1F("h_raw","Acquisizione",NBINS,0,MAXHISTOCHN);
	UInt_t toentry=branch->GetEntries();
	printf("getHistoFromFile: There are %d entries in the branch\n",toentry);
	for(i=0; i<toentry; i++) {
		branch->GetEntry(i);
		h_raw->Fill(inc_data.qlong);
	}

  h_raw->Draw();
    
	TSpectrum	*s = new TSpectrum(10);
	Int_t		e, nPeaks, bTemp, bFirstPeak = 9999;
	Float_t	*xPeaks;
    
	// trovo il primo picco
	nPeaks = s->Search(h_raw->Rebin(2, "h_raw_rebinned"));
	if (nPeaks>0) {
		xPeaks = s->GetPositionX();
		// loop sui picchi per trovare il primo
		for (i=0;i<nPeaks;i++) {
			bTemp = h_raw->GetXaxis()->FindBin(xPeaks[i]);
			if (bTemp<bFirstPeak) { bFirstPeak = bTemp; }
		}
	} else { bFirstPeak = 0; }
    // sottraggo il fondo Compton
	Float_t 		*bgBins = new Float_t[NBINS];
	TH1F			*h_filtered = new TH1F("h_filtered","Picchi",NBINS,0,MAXHISTOCHN);
	/*for (i = 0; i < bFirstPeak; i++) { bgBins[i]=h_raw->GetBinContent(i+1); }
	s->Background(bgBins,bFirstPeak,20,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder2,kFALSE,TSpectrum::kBackSmoothing15,kFALSE);
	for (i = 0; i < bFirstPeak; i++) { h_filtered->SetBinContent(i+1, (h_raw->GetBinContent(i+1)-bgBins[i])); }
	for (i = 0; i < NBINS-bFirstPeak; i++) { bgBins[i]=h_raw->GetBinContent(bFirstPeak+i+1); }
	s->Background(bgBins,NBINS-bFirstPeak,20,TSpectrum::kBackDecreasingWindow,TSpectrum::kBackOrder2,kFALSE,TSpectrum::kBackSmoothing15,kFALSE);
	for (i = bFirstPeak; i < NBINS; i++) { h_filtered->SetBinContent(i+1, (h_raw->GetBinContent(i+1)-bgBins[i-bFirstPeak])); }
	//TCanvas * background = new TCanvas("background","Estimation of bg",10,10,1000,700);*/
	
  for (i = 0; i < NBINS; i++) { bgBins[i]=h_raw->GetBinContent(i+1); }

	s->Background(bgBins, NBINS, 20, TSpectrum::kBackDecreasingWindow, TSpectrum::kBackOrder2, kFALSE, TSpectrum::kBackSmoothing15, kFALSE);

	for (i = 0; i < NBINS; i++) { h_filtered->SetBinContent(i+1, (h_raw->GetBinContent(i+1)-bgBins[i])); }
	
  h_raw->Draw("L");
	h_filtered->SetLineColor(kRed);
	h_filtered->Draw("SAME L");
  c0->Update();
		// trovo i picchi e calibro
    
    
	nPeaks = s->Search(h_filtered->Rebin(2, "h_filtered_rebinned"),2,"",0.05);
    xPeaks = s->GetPositionX();
	//	sigma_fall[0] = 100;


  timer->TurnOn();
  timer->Reset();
  timer->TurnOff();



	if (nPeaks<NUMENERGIES) {
		// trovati troppi pochi picchi