C++ TH1::GetFunction方法代码示例

double Fit511Photopeak(TH1* h, double* error=0)
{
  TSpectrum spec(1);
  spec.Search(h);
  h->GetXaxis()->SetTitle("Energy [photoelectrons]");
  h->Draw("e");
  TH1* bg = spec.Background(h);
  TH1* sig = (TH1*)(h->Clone());
  sig->SetLineColor(kGreen);
  
  sig->Add(bg,-1);
  sig->Draw("same e");
  sig->Fit("gaus","m","same e");
  TF1* gaus = sig->GetFunction("gaus");
  if(gaus)
    gaus->SetLineColor(kGreen);
  
  bg->SetLineColor(kRed);
  bg->Draw("same e");
  
  TLine* line = new TLine(gaus->GetParameter(1),0,
			  gaus->GetParameter(1),h->GetMaximum());
  line->SetLineColor(kBlue);
  line->SetLineWidth(2);
  line->Draw();
  

  double yield = spec.GetPositionX()[0]/epeak;
  double err = 0;
  
  cout<<"Results from TSpectrum: \n\t"
      <<"Peak = "<<spec.GetPositionX()[0]<<" p.e.; Light Yield = "
      <<yield<<" p.e./keV"<<endl;
  if(gaus){
    yield = gaus->GetParameter(1)/epeak;
    err = gaus->GetParError(1)/epeak;
    cout<<"Results from BG Subtracted Gaus Fit: \n\t"
	<<"Peak = "<<gaus->GetParameter(1)<<" p.e.; Light Yield = "
	<<yield<<" +/- "<<err<<" p.e./keV"<<endl;
    err = max(err, TMath::Abs(yield-spec.GetPositionX()[0]/epeak));
    
  }
  TLegend* leg = new TLegend(.6,.6,.9,.9);
  leg->AddEntry(h,"Raw Spectrum","lpe");
  leg->AddEntry(bg,"Background","lpe");
  leg->AddEntry(sig,"Signal","lpe");
  char title[20];
  sprintf(title,"Yield = %.2f pe/keV",yield);
  leg->AddEntry(line, title ,"l");
  leg->Draw();

  if(error) *error = err;
  return yield;
}

Double_t fitep0sigma( char* hs, int binlow=-999, int binhigh=999) {

  TH1 *h = (TH1*)gDirectory->Get(hs);
  if( h == NULL ){ cout << hs << " does not exist\n"; return 0; }

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1, x9;
  if(binlow < -900 && binhigh > 900) {
    x1 = h->GetBinCenter(1);
    x9 = h->GetBinCenter(nb);
  }
  else {
    x1 = binlow;
    x9 = binhigh;
  }

  // create a TF1 with the range from x1 to x9 and 5 parameters
  TF1 *ep0Fcn = new TF1( "ep0Fcn", ep0Fit, x1, x9, 5 );

  ep0Fcn->SetParName( 0, "mean" );
  ep0Fcn->SetParName( 1, "sigma" );
  ep0Fcn->SetParName( 2, "pow" );
  ep0Fcn->SetParName( 3, "area" );
  ep0Fcn->SetParName( 4, "BG" );

  // Start values for some parameters:
  ep0Fcn->SetParameter( 0, xmax ); // peak position
  ep0Fcn->SetParameter( 1, 4*dx ); // width
  ep0Fcn->SetParameter( 2, 3.3 ); // pow
  ep0Fcn->SetParameter( 3, nn ); // N
  ep0Fcn->SetParameter( 4, bg );
    
  h->Fit("ep0Fcn", "Q R", "ep" );
  TF1 *fit = h->GetFunction("ep0Fcn");
  return fit->GetParameter(1);
}

double fittp0sigma( char* hs ) {

  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ){ cout << hs << " does not exist\n"; return 0; }

  double dx = h->GetBinWidth(1);
  double nmax = h->GetBinContent(h->GetMaximumBin());
  double xmax = h->GetBinCenter(h->GetMaximumBin());
  double nn = 7*nmax;

  int nb = h->GetNbinsX();
  double n1 = h->GetBinContent(1);
  double n9 = h->GetBinContent(nb);
  double bg = 0.5*(n1+n9);

  double x1 = h->GetBinCenter(1);
  double x9 = h->GetBinCenter(nb);

  // create a TF1 with the range from x1 to x9 and 5 parameters
  TF1 *tp0Fcn = new TF1( "tp0Fcn", tp0Fit, x1, x9, 5 );

  tp0Fcn->SetParName( 0, "mean" );
  tp0Fcn->SetParName( 1, "sigma" );
  tp0Fcn->SetParName( 2, "nu" );
  tp0Fcn->SetParName( 3, "area" );
  tp0Fcn->SetParName( 4, "BG" );
   
  // set start values for some parameters:
  tp0Fcn->SetParameter( 0, xmax ); // peak position
  tp0Fcn->SetParameter( 1, 4*dx ); // width
  tp0Fcn->SetParameter( 2, 2.2 ); // nu
  tp0Fcn->SetParameter( 3, nn ); // N
  tp0Fcn->SetParameter( 4, bg );
    
  h->Fit( "tp0Fcn", "Q R", "ep" );
  // h->Fit("tp0Fcn","V+","ep");
  TF1 *fit = h->GetFunction("tp0Fcn");
  return fit->GetParameter(1);
}

//.........这里部分代码省略.........
			mean = hCell->GetMean();
			rms = hCell->GetRMS();
			hCellSub = (TH1F *) hCellSubList->At(h * DOMINO_NCELL + ch);
			meansub = hCellSub->GetMean();
			rmssub = hCellSub->GetRMS();
			fgauss->SetParameter(0, (Double_t) nevt / 4.);
			fgauss->SetParameter(1, mean);
			fgauss->SetParameter(2, rms);
			//            hCell->Fit("fgauss","QN0");
			grData->SetPoint(ch, ch, mean);
			grData->SetPointError(ch, 0, rms);
			grDataSub->SetPoint(ch, ch, meansub);
			//            grDataSub->SetPointError(ch,0.5,rmssub);
			grDataSub->SetPointError(ch, 0.5, 2.1);
			hRMSdist->Fill(rms);
			hRMSFitdist->Fill(fgauss->GetParameter(2));
			//           cout << "cell:" << ch << " index:" << h*DOMINO_NCELL+ch << " Mean,RMS:"<< mean << "," << rms<<endl;
		}
	}

	Double_t x, y, chtmp, x1, x2, y1, y2;

	/*TList *grCellCalibList = OpenCalibFile("CalibrationData1000events.root");

	TGraphErrors *grCellCalib;
	TGraphErrors *grDataSubCalib = new TGraphErrors(DOMINO_NCELL);
	grDataSubCalib->SetTitle("Data after calibration correction");
	grDataSub = (TGraphErrors *) grDataSubList->At(anaChannel);


	for(ch = 0; ch < DOMINO_NCELL; ch++) {
		grCellCalib = ((TGraphErrors *) grCellCalibList->At(ch));
		grCellCalib->Fit("pol3", "Q");
		TF1 *pol3fit = ((TF1 *) grCellCalib->GetFunction("pol3"));
		grDataSub->GetPoint(ch, x, y);
		chtmp = y - (Double_t)(pol3fit->Eval(y/3.25));
		grDataSubCalib->SetPoint(ch, x, chtmp);
	}

	TCanvas *cGrTest = new TCanvas("grTest", "test per vedere i dati", 1000,1000);

	grDataSubCalib->Draw("APEL");*/


	TString Title = "Charge Distribution per channel";
	gStyle->SetOptFit(111);
	TCanvas *cdistch = new TCanvas("cdistch", Title, 1000, 1000);
	cdistch->Divide(3, 3);
	for (int i = 0; i < DOMINO_NCH; i++) {
		cdistch->cd(i + 1);
		TH1 *dhist = (TH1 *) DistChList->At(i);
		dhist->DrawCopy();
		dhist->SetLineWidth(1);
		dhist->Fit("gaus", "Q");
		dhist->GetFunction("gaus")->SetLineColor(4);
		dhist->GetFunction("gaus")->SetLineWidth(2);
	}

	TString Title = "Charge Distribution Pedestals Subtracted per channel";
	TCanvas *cdistchsub = new TCanvas("cdistchsub", Title, 1000, 1000);
	cdistchsub->Divide(3, 3);
	for (int i = 0; i < DOMINO_NCH; i++) {
		cdistchsub->cd(i + 1);
		TH1 *dsubhist = (TH1 *) DistChSubList->At(i);
		dsubhist->DrawCopy();
		dsubhist->SetLineWidth(1);

//------------------------------------------------------------------------------
void PlotAlignmentValidation::plotSubDetResiduals(bool plotNormHisto,unsigned int subDetId)
{
  setNiceStyle();
 
  gStyle->SetOptStat(11111);
  gStyle->SetOptFit(0000);

  TCanvas *c = new TCanvas("c", "c", 600,600);
  c->SetTopMargin(0.15);
  TString histoName= "";
  if (plotNormHisto) {histoName= "h_NormXprime";}
  else histoName= "h_Xprime_";
  switch (subDetId){
  case 1 : histoName+="TPBBarrel_0";break;
  case 2 : histoName+="TPEendcap_1";break;
  case 3 : histoName+="TPEendcap_2";break;
  case 4 : histoName+="TIBBarrel_0";break;
  case 5 : histoName+="TIDEndcap_1";break;
  case 6 : histoName+="TIDEndcap_2";break;
  case 7 : histoName+="TOBBarrel_3";break;
  case 8 : histoName+="TECEndcap_4";break;
  case 9 : histoName+="TECEndcap_5";break;
  }
  int tmpcounter = 0;
  TH1 *sumHisto = 0;
  for(std::vector<TkOfflineVariables*>::iterator it = sourceList.begin();
      it != sourceList.end(); ++it) {
    if (tmpcounter == 0 ) {
      TFile *f= (*it)->getFile();
      sumHisto =(TH1*) f->FindKeyAny(histoName)->ReadObj();//FindObjectAny(histoName.Data());
      sumHisto->SetLineColor(tmpcounter+1);
      sumHisto->SetLineStyle(tmpcounter+1);
      sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
      sumHisto->Draw();
      
      //get statistic box coordinate to plot all boxes one below the other
      //gStyle->SetStatY(0.91);
      //gStyle->SetStatW(0.15);
      //gStyle->SetStatBorderSize(1);
      //gStyle->SetStatH(0.10);
      
      
      tmpcounter++;
    } else {
      sumHisto = (TH1*) (*it)->getFile()->FindObjectAny(histoName);
      sumHisto->SetLineColor(tmpcounter+1);
      sumHisto->SetLineStyle(tmpcounter+1);
      sumHisto->GetFunction("tmp")->SetBit(TF1::kNotDraw);
      //hstack->Add(sumHisto);
      
      c->Update();
      tmpcounter++;  
    }
    TObject *statObj = sumHisto->GetListOfFunctions()->FindObject("stats");
    if (statObj && statObj->InheritsFrom(TPaveStats::Class())) {
      TPaveStats *stats = static_cast<TPaveStats*>(statObj);
      stats->SetLineColor(tmpcounter+1);
      stats->SetTextColor(tmpcounter+1);
      stats->SetFillColor(10);
      stats->SetX1NDC(0.91-tmpcounter*0.1);
      stats->SetX2NDC(0.15);
      stats->SetY1NDC(1);
      stats->SetY2NDC(0.10);
      sumHisto->Draw("sames");
    }
  }
  //hstack->Draw("nostack");
  char PlotName[1000];
  sprintf( PlotName, "%s/%s.eps", outputDir.c_str(), histoName.Data() );
  
  c->Print(PlotName);
  //delete c;
  //c=0;
    
}

//.........这里部分代码省略.........

	gProgress->Reset();
	gProgress->SetMax(nevt);

	gSystem->ProcessEvents();

	while (ievt <= nevt && !flagEnd) {
		fread((void *) &event_data, 1, sizeof(event_data), fdata);
		if (feof(fdata))
			flagEnd = 1;

		p = (struct channel_struct *) &event_data.ch[0]; // read bunch of data
		dep = (struct channel_struct *) &event_data.ch[1]; // read bunch of data
		// goes to channel to analyze
		p += anaChannel;

		// read data, subtract pedestals values and save results in grAnaChDataTemp graph with
		// fixed error for each point (x = 0.5 and y = 2.1). Also generate an array with Domino
		// X and Y values

		TGraphErrors *grAnaChDataTemp = new TGraphErrors(DOMINO_NCELL);

		for (int ch = 0; ch < DOMINO_NCELL; ch++) {
			// Read pedestal value for this cell
			grPed->GetPoint(ch, itmp, PedVal);
			chtmp = (Double_t)(p->data[ch]); // data value
			chtmp = chtmp - PedVal;
			grAnaChDataTemp->SetPoint(ch, (Double_t) ch, chtmp);
			grAnaChDataTemp->SetPointError(ch, 0.5, 2.1);
		}
		// create fit functions
		TF1 *fsin = new TF1("fsin", sigSin, 0., 1024., 4);
		fsin->SetParameters(600., 255., 150., 150.);
		fsin->SetParNames("amplitude", "Period", "Phase", "DC-Offset");

		grAnaChDataTemp->Fit("fsin", "Q");
		TF1 *fsinFit = grAnaChDataTemp->GetFunction("fsin");
		fsinFit->SetParNames("amplitude", "Period", "Phase", "DC-Offset");

		// debug
		cfitTest->cd(ievt);
		grAnaChDataTemp->SetMarkerStyle(20);
		grAnaChDataTemp->SetMarkerSize(0.3);
		grAnaChDataTemp->GetYaxis()->SetLabelSize(0.12);
		grAnaChDataTemp->GetXaxis()->SetLabelSize(0.12);
		grAnaChDataTemp->Draw("APE");

		Double_t fitPeriod, fitAmplitude, chisquare;
		fitPeriod = fsinFit->GetParameter("Period");
		fitAmplitude = TMath::Abs(fsinFit->GetParameter("amplitude"));
		chisquare = fsinFit->GetChisquare();

		cout << "period: " << fitPeriod << " amplitude: " << fitAmplitude << " chisquare: " << chisquare << endl;

		if(chisquare > 0.1e+06) {
			gProgress->Increment(1);
			gSystem->DispatchOneEvent(kTRUE);
			ievt++;
			continue;
		}

		gProgress->Increment(1);
		gSystem->DispatchOneEvent(kTRUE);

		hPeriod->Fill(fitPeriod);
		fitusati++;

		ievt++;

	}

	cout << "fit scartati :" << nevt - fitusati << endl;
	//draw
	TString Title = "Period distribution for nevt events";
	TCanvas *cPeriod = new TCanvas("cPeriod", Title, 700, 700);
	hPeriod->Draw();
	hPeriod->Fit("gaus");

	TF1 *fgausFit = hPeriod->GetFunction("gaus");
	//mean = fgausFit->GetParameter(1);
	//	rms = fgausFit->GetParameter(2);

	mean = hPeriod->GetMean();
	rms = hPeriod->GetRMS();

	TString OutFile = "Period";
	OutFile += nevt;
	OutFile += "events.dat";
	FILE *f = fopen(OutFile.Data(), "w");
	fwrite(&mean, sizeof(mean), 1, f);
	fwrite(&rms, sizeof(rms), 1, f);

	((TGMainFrame *) gProgress->GetParent())->CloseWindow();
	fclose(f);

	cout << "mean: " << mean << " rms: " << rms << endl;

	fclose(fdata);

}