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

int e428ana1::MakePara(Double_t thr){
    int npeak,nlap;
    double threshol,thrmin,thrmax;
    TSpectrum *spec;
    for (int i=0; i<4; i++) {
        for (int j=0; j<127; j++) {
            npeak=0;nlap=0;
            threshol=thr;
            thrmax=1;thrmin=0;
            while (npeak!=4 && nlap<100){
                spec = new TSpectrum(10);
                npeak = spec->Search(hdssd[i][j],2,"",threshol);
                if (npeak >3) {
                    thrmin=threshol;
                }else thrmax=threshol;
                threshol=(thrmax-thrmin)/2;
                nlap++;
        //        printf("Npeak   %d   %d    %d    %f\n",i,j,npeak,threshol);
            }
            float *xpeak=spec->GetPositionX();
            float *ypeak=spec->GetPositionY();
            printf("Npeak   %d   %d    %d",i,j,npeak);
            for (int i1=0; i1<npeak; i1++) {
                printf("    %f   %f",xpeak[i1],ypeak[i1]);
            }
            printf("\n");
        }
    }    
    return 0;
}

////    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;
}

I2GFvalues I2GFmainLoop(TH1F *htemp, int N_iter, float N_sigma_range, bool ShowFit)
//Arguments: (histo to be fit, N iterations to find peak using gaus fit, fit range param., do or do not plot fit on histo)
{
  I2GFvalues myI2GFvalues;

  //Set initial values...(in case fit fails)
  myI2GFvalues.rchi2 = -100;
  myI2GFvalues.mean = -100;
  myI2GFvalues.mean_err = -100;
  myI2GFvalues.sigma = -100;
  myI2GFvalues.sigma_err = -100;


  TSpectrum *s = new TSpectrum(); //TSpectrum(1,1)->Argument: (Number of peaks to find, Distance to neighboring peak: "1"-->3sigma)
  Int_t NPeaks;
  Float_t *Peak;                     //TSpectrum *s = new TSpectrum(); --> No warning message 
  Float_t *PeakAmp;                    
  float peak_pos = 0;
  float peak_pos_amp = 0;  //Initial value, assuming positive going peaks
  int peak_pos_bin = 0;
  
  int binMaxCnt = 0;
  float binMaxCnt_value = 0;
  int binMaxCnt_counts = 0;

  int Nbins = 0;
  Int_t zero_value_bin = 0;
  float low_limit = 0;
  float high_limit = 0;
  float peak_pos_count = 0;
  float zero_bin_value = 0;
  float max_bin_value = 0;
 

  TF1 *func;
  TF1 *func1;
  TF1 *func2;
  TF1 *func3;
  float Chi2;
  int NDF = 1;
   

  float f_RChi2 = 1; 
  float f_const = 1;
  float f_mean = 1;
  float f_mean_err = 1;
  float f_sigma = 1;
  float f_sigma_err = 1;
  
  float peak = 1;  

  float f_const2 = 1;
  float f_mean2 = 1;
  float f_mean_err2 = 1;
  float f_sigma2 = 1;
  float f_sigma_err2 = 1;

  //---------Basic Histo Peak Finding Parameters----------
  binMaxCnt = htemp->GetMaximumBin();  //Finds bin w/ max counts (can be dubious, so use TSpectrum)
  binMaxCnt_value = (Double_t) htemp->GetXaxis()->GetBinCenter(binMaxCnt); //if the bin number is known and the bin value (in x-axis units) is wanted
  binMaxCnt_counts = (Int_t) htemp->GetBinContent(binMaxCnt); //Finds counts within particular bin

  //---------TSpectrum Peak Finding Parameters--------
  if (ShowFit) NPeaks = s->Search(htemp, 2, "goff", 0.5); //opens a canvas (one time in a loop), even with:  s->Search(htemp, 2, "nodraw", 0.9);
  //else  NPeaks = s->Search(htemp, 2, "", 0.5);  //s->Search(htemp, 2, "nodraw", 0.9);

  //Npeaks = s->GetNPeaks(); //If using this, need pointer in declaration above: Int_t *NPeaks
  //N_peaks =  Npeaks[0];

  Peak = s->GetPositionX();
  PeakAmp = s->GetPositionY();

  for (int i=0; i<NPeaks; i++)
    {
      if (peak_pos_amp < PeakAmp[i])
	{
	    peak_pos_amp = PeakAmp[i]; //TSpectrum finds peak counts
	    peak_pos = Peak[i]; //TSpectrum finds pos. of peak in x-axis units
	}
    }
  peak_pos_bin = htemp->GetXaxis()->FindBin(peak_pos); //if the bin value (in x-axis units) is known and the bin number is wanted
  peak_pos_count = htemp->GetBinContent(peak_pos_bin);  //counts in particular bin

  //------------------------------------------------------------------------------------------------------------------  

  zero_value_bin = htemp->GetXaxis()->FindBin(0.0); //if the bin value (in x-axis units) is known and the bin number is wanted
  Nbins = htemp->GetSize() - 2; //total number of bins in histo
  zero_bin_value =  htemp->GetXaxis()->GetBinCenter(0); //if the bin number is known and the bin value (in x-axis units) is wanted.
  max_bin_value =  htemp->GetXaxis()->GetBinCenter(Nbins); //if the bin number is known and the bin value (in x-axis units) is wanted.
 
  int TS = 0;
  if (peak_pos >= zero_bin_value  &&  peak_pos <= max_bin_value)  //Make sure that TSpectrum peak is within histo range
    {                                                             //if not, use Par initial values from Basic Histo Peak Find
             
      TS=1; //for cout below                        
      low_limit = peak_pos - (0.1 * abs(max_bin_value-zero_bin_value)); //peakpos-10% of full range of histo
      high_limit = peak_pos + (0.1 * abs(max_bin_value-zero_bin_value)); //peakpos+10% of full range of histo
                                                       
      func = new TF1("func", "gaus");
      //func->FixParameter(1,0);
//.........这里部分代码省略.........

void fit(){
  int range = 1e4;

  TFile *f = new TFile("hgrroot.root");
  TCanvas *ca = new TCanvas("ca","ca",0,0,1400,900);
  ca->Divide(4,4);
  TH1F* h[7][4];
  TF1* fu[7][4][2];
  TF1* fus[7][4][2][3];
  double res[7*4];
  double det[7*4];
  double slope[7*4];
  double offset[7*4];
  double resolution =5;
  TEnv* output = new TEnv("corecal.dat");

  for(int d=0;d<7;d++){
    for(int c=0;c<4;c++){
      h[d][c] = (TH1F*)f->Get(Form("htraceen_b%d_c%d_cr%d_d%d",0,9,c,d+1));
      ca->cd(1+c*4);
      //h[d][c]->GetXaxis()->SetRangeUser(100,4000);
      TSpectrum *sp = new TSpectrum(2,resolution);
      sp->SetResolution(resolution);
      Int_t nfound = 0;
      nfound = sp->Search(h[d][c],resolution,"nobackground",0.5);
      Float_t *xpeaks = sp->GetPositionX();
      Float_t *ypeaks = sp->GetPositionY();
//       for(int p=0;p<nfound;p++){
// 	cout << xpeaks[p] << "\t" << ypeaks[p] << endl;
//       }
      if(nfound!=2){
	cout << "Found " << nfound << " peaks in spectrum, too many, aborting" << endl;
	continue;
      }
      h[d][c]->DrawCopy();
      //check if first peak is lower in energy, otherwise swap them
      if(xpeaks[0]>xpeaks[1]){
	Float_t temp = xpeaks[1];
	xpeaks[1] = xpeaks[0];
	xpeaks[0] = temp;
	temp = ypeaks[1];
	ypeaks[1] = ypeaks[0];
	ypeaks[0] = temp;
	
      }

      for(int p=0;p<nfound;p++){
	ca->cd(1+c*4+1+p);
	h[d][c]->GetXaxis()->SetRangeUser(xpeaks[p]-range,xpeaks[p]+range);
	h[d][c]->DrawCopy();
	fu[d][c][p] = new TF1(Form("fcore_d%d_c%d_p%d",d,c,p),fgammagaussbg,xpeaks[p]-range,xpeaks[p]+range,6);
	fu[d][c][p]->SetLineColor(3);
	fu[d][c][p]->SetLineWidth(1);
	fu[d][c][p]->SetParameter(0,0);//bg const
	fu[d][c][p]->SetParameter(1,0);//bg slope
	fu[d][c][p]->SetParameter(2,h[d][c]->Integral(xpeaks[p]-range,xpeaks[p]+range));//norm
	fu[d][c][p]->SetParameter(3,xpeaks[p]);//mean
	fu[d][c][p]->SetParLimits(3,xpeaks[p]-500,xpeaks[p]+500);//mean
	fu[d][c][p]->SetParameter(4,100);//sigma
	fu[d][c][p]->SetParLimits(4,0.001,1000);//sigma
	fu[d][c][p]->SetParameter(5,h[d][c]->GetBinContent(h[d][c]->FindBin(xpeaks[p]-range)));//step
	
	h[d][c]->Fit(fu[d][c][p],"Rqn");
	fu[d][c][p]->Draw("same");


	fus[d][c][p][0] = new TF1(Form("fcore_d%d_c%d_p%d_bg",d,c,p),fgammabg,xpeaks[p]-range,xpeaks[p]+range,6);
	fus[d][c][p][1] = new TF1(Form("fcore_d%d_c%d_p%d_st",d,c,p),fgammastep,xpeaks[p]-range,xpeaks[p]+range,6);
	fus[d][c][p][2] = new TF1(Form("fcore_d%d_c%d_p%d_ga",d,c,p),fgammagaus,xpeaks[p]-range,xpeaks[p]+range,6);
	

	fus[d][c][p][0]->SetLineColor(5);
	fus[d][c][p][1]->SetLineColor(4);
	fus[d][c][p][2]->SetLineColor(2);
	for(int k=0;k<3;k++){
	  fus[d][c][p][k]->SetLineWidth(1);
	  for(int l=0;l<6;l++)
	    fus[d][c][p][k]->SetParameter(l,fu[d][c][p]->GetParameter(l));
	  fus[d][c][p][k]->Draw("same");
	}

      }//peaks
      //res[d*4+c] = 2.35*fu[d][c][1]->GetParameter(4)*(1332.492-1173.228)/(fu[d][c][1]->GetParameter(3)-fu[d][c][0]->GetParameter(3));

      slope[d*4+c] = (1332.492-1173.228)/(fu[d][c][1]->GetParameter(3)-fu[d][c][0]->GetParameter(3));
      offset[d*4+c] = (1332.492+1173.228)-slope[d*4+c]*(fu[d][c][1]->GetParameter(3)+fu[d][c][0]->GetParameter(3));
      offset[d*4+c]*=0.5;
      
      //cout << fu[d][c][0]->GetParameter(3)*slope[d*4+c] + offset[d*4+c] << "\t" << fu[d][c][1]->GetParameter(3)*slope[d*4+c] + offset[d*4+c] << endl;
      
      res[d*4+c] = fu[d][c][1]->GetParameter(2);
      det[d*4+c] = d*5+c;

      output->SetValue(Form("Slope.d%d.c%d",d,c),slope[d*4+c]);
      output->SetValue(Form("Offset.d%d.c%d",d,c),offset[d*4+c]);


      //cout << det[d*4+c] <<"\t" << res[d*4+c] << endl;
      //cout << fu[d][c][0]->GetParameter(4) <<"\t" << fu[d][c][1]->GetParameter(4) << endl;
      //ca->cd(1+c*4+3);
//.........这里部分代码省略.........

//------------------------------------
//function to fit energy distributions
energyRes* get_res(int snum, Double_t energy, bool do_pion, bool use_f_pion, bool do_fit, bool do_show, bool do_print=false, bool do_batch=false){
	Sample* sp = sample_map[snum];
	if(!sp) { std::cout << "Sample " << snum << " is not loaded." << std::endl; energyRes* theRes = new energyRes(0,0); return theRes; }
	
	//select correct file
	std::string fpre = sp->fpre;
	if(do_pion) fpre += "_pion";
	else fpre += "_elec";

	//make filenames
	std::stringstream drawname, fname, piname;
	fname << sp->dir << "/" << fpre << "_" << energy << "gev_10k.root";
	if(do_pion) piname << "#pi^{-} " << energy << " GeV";
	else piname << "e^{-} " << energy << " GeV";

	//open file and tree
	TFile* _file;
	_file = TFile::Open((fname.str()).c_str());
	TTree* totalTree = (TTree*)_file->Get("Total");

	//default histo settings
	//double Emin = 0.1*energies[num]; //lower cut to remove weird peaks near E=zero
	double Emin = 0;
	double Emax = 2*energy;
	int nbins = 100;
	
	//ecal & hcal energies need to be calibrated
	get_sampling_factors(snum);

	//make tree drawing expressions
	//define mip as ecal < 1 gev = 1000 mev
	if(use_f_pion) drawname << sp->sam_pion;
	else drawname << sp->sam_elec;
	
	if(sp->det==Hcal) drawname << "*(hcal+" << sp->zeroWt << "*zero)/1000";
	else drawname << "*ecal/1000";

	drawname << ">>htemp(" << nbins << "," << Emin << "," << Emax << ")";
	//std::cout << drawname.str() << std::endl;

	TH1F* h_res; //to store histos drawn from tree
	TF1* gfit;
	TF1* gsnL;
	TF1* gsnR;

	//plotting variables
	TCanvas* can;
	TPad* pad;
	TLegend* leg;
	TPaveText* pave;
	TPaveText* pave_par;
	TLine *aLline;
	TLine *aRline;

	//create instance of energyRes object
	energyRes* theRes = new energyRes(energy,2);

	//draw w/ appropriate cut
	totalTree->Draw((drawname.str()).c_str(),"","hist goff");
	h_res = (TH1F*)gDirectory->Get("htemp");
	h_res->SetTitle("");
	h_res->GetXaxis()->SetTitle("Energy [GeV]");
	h_res->SetLineWidth(2);
	h_res->SetLineColor(kBlack);

	//names
	std::string ofit;
	if(do_fit) ofit = "fit";
	else ofit = "nofit";
	std::stringstream oname;
	oname << pdir << "/" << fpre;
	if(use_f_pion) oname << "_fpion";
	oname << "_response_" << ofit << "_" << energy << "gev";
	
	//get values from histo
	Double_t m = h_res->GetMean();
	Double_t me = h_res->GetMeanError();
	//Double_t m = h_res->GetBinCenter(h_res->GetMaximumBin()); //peak
	Double_t s = h_res->GetRMS();
	Double_t se = h_res->GetRMSError();
	Int_t N = h_res->GetEntries();
	
	std::vector<Double_t> stats(3,0);
	std::vector<Double_t> stat_err(3,0);
	stats[0] = N;
	stat_err[0] = 0;
	stats[1] = m;
	stat_err[1] = me;
	stats[2] = s;
	stat_err[2] = se;

	//find peak
	TSpectrum *spec = new TSpectrum(5);
	if(nbins < 100) spec->Search(h_res,6,"nobackground nodraw goff"); //turn off background removal when nbins too small
	else spec->Search(h_res,6,"nodraw goff");
	Float_t* xpos = spec->GetPositionX();
	Float_t* ypos = spec->GetPositionY();
	Double_t p = xpos[0];
//.........这里部分代码省略.........

//--------------------------------------
//function to calculate sampling factors
std::pair<Double_t,Double_t> g4_sample(int snum, Double_t energy, bool do_pion, bool do_show, bool do_print=false, bool set_val=true){
	Sample* sp = sample_map[snum];
	if(!sp) { std::cout << "Sample " << snum << " is not loaded." << std::endl; return std::pair<Double_t,Double_t>(0.,0.); }

	//select correct file
	std::string fpre = sp->fpre;
	if(do_pion) fpre += "_pion";
	else fpre += "_elec";

	//make filenames
	std::stringstream drawname, fname, piname;
	fname << sp->dir << "/" << fpre << "_" << energy << "gev_10k.root";
	if(do_pion) piname << "#pi^{-} " << energy << " GeV";
	else piname << "e^{-} " << energy << " GeV";

	//open file and tree
	TFile* _file;
	_file = TFile::Open((fname.str()).c_str());
	TTree* totalTree = (TTree*)_file->Get("Total");

	//get histo from tree (no display)
	//define mip as sam_ecal*ecal < 1 gev = 1000 mev (for pions in HCAL)
	if(sp->det==Hcal) drawname << "(hcal+" << sp->zeroWt << "*zero)/1000>>hsam(200)";
	else drawname << "(ecal)/1000>>hsam(200)";
	
	totalTree->Draw((drawname.str()).c_str(),"","hist goff");
	TH1F* hsam = (TH1F*)gDirectory->Get("hsam");
	
	//use parameters from histo to start fit
	TSpectrum* spec = new TSpectrum(5);
	spec->Search(hsam,5,"nodraw goff");
	Float_t* xpos = spec->GetPositionX();
	Float_t* ypos = spec->GetPositionY();

	Double_t m = xpos[0];
	Double_t me = hsam->GetMeanError();
	Double_t N = hsam->GetEntries();
	std::stringstream s_mean;
	s_mean.precision(3);
	Double_t f = energy/m;
	Double_t f_err = energy*(me/(m*m));
	s_mean << f << " #pm " << f_err;

	TPolyMarker* pm = new TPolyMarker(1, xpos, ypos);
	hsam->GetListOfFunctions()->Add(pm);
	pm->SetMarkerStyle(23);
	pm->SetMarkerColor(kRed);
	pm->SetMarkerSize(1.3);

	std::cout.precision(6);
	std::cout << "f_" << (do_pion ? "pion" : "elec") << " = " << f << " +/- " << f_err << std::endl;
	
	//plotting and printing
	if (do_show){
		TCanvas* can = new TCanvas("sample","sample",700,500);
		can->cd();
		TPad* pad = new TPad("graph","",0,0,1,1);
		pad->SetMargin(0.12,0.05,0.15,0.05);
		pad->Draw();
		pad->cd();
		
		//formatting
		hsam->SetTitle("");
		hsam->GetXaxis()->SetTitle("Energy [GeV]");
		//hsam->SetStats(kTRUE);
		//gStyle->SetOptStat("mr");
		hsam->SetLineWidth(2);
		hsam->SetLineColor(kBlack);
		hsam->GetYaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetYaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetTitleSize(32/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetLabelSize(28/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetYaxis()->SetTickLength(12/(pad->GetWh()*pad->GetAbsHNDC()));
		hsam->GetXaxis()->SetTickLength(12/(pad->GetWh()*pad->GetAbsHNDC()));
		
		hsam->Draw();
		
		std::stringstream Nname;
		Nname << "N = " << N;
		
		//determine placing of pave
		Double_t xmin;
		if (m/((hsam->GetXaxis()->GetXmax() + hsam->GetXaxis()->GetXmin())/2) < 1) xmin = 0.65;
		else xmin = 0.2;
		
		//legend
		TPaveText *pave = new TPaveText(xmin,0.65,xmin+0.2,0.85,"NDC");
		pave->AddText((piname.str()).c_str());
		pave->AddText((Nname.str()).c_str());
		pave->AddText("Peak sampling factor:");
		pave->AddText((s_mean.str()).c_str());
		pave->SetFillColor(0);
		pave->SetBorderSize(0);
		pave->SetTextFont(42);
		pave->SetTextSize(0.05);
		pave->Draw("same");

		if(do_print) {
//.........这里部分代码省略.........