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

void draw_ToF_Calib_Tower() {

  ofstream fout("ToF_Calib_Tower.txt");
  TFile* f = new TFile("/phenix/plhf/zji/taxi/Run13pp510MinBias/8328/data/total.root");
  TH2I* h_hitmap_tof_energy_PbSc = (TH2I*)f->Get("hitmap_tof_energy_PbSc");
  TH2I* h_hitmap_tof_energy_PbGl = (TH2I*)f->Get("hitmap_tof_energy_PbGl");
  TH2I* h_hitmap_tof_tower = (TH2I*)f->Get("hitmap_tof_tower");

  TCanvas* c = new TCanvas("c", "Canvas", 1800, 600);
  gStyle->SetOptStat(1);
  gStyle->SetOptFit(1);
  c->Divide(3,1);

  c->cd(1);
  h_hitmap_tof_energy_PbSc->Draw("colz");

  c->cd(2);
  h_hitmap_tof_energy_PbGl->Draw("colz");

  c->cd(3);
  fout << "towerid    ToF peak" << endl; 
  for(Int_t id=0; id<25000; id++) {
    TH1D* hp_hitmap_tof_tower = (TH1D*)h_hitmap_tof_tower->ProjectionY("_py",id+1,id+1);
    TSpectrum* peak = new TSpectrum();
    Int_t nfound = peak->Search(hp_hitmap_tof_tower,2.,"nodraw");
    if(nfound) {
      Float_t peakX = *peak->GetPositionX();
      hp_hitmap_tof_tower->Fit("gaus","Q0","",peakX-3.,peakX+3.);
      TF1* f_hitmap_tof_tower = hp_hitmap_tof_tower->GetFunction("gaus");
      Double_t mean_hitmap_tof_tower = f_hitmap_tof_tower->GetParameter(1);
      fout << id << "    " << mean_hitmap_tof_tower << endl;
      f_hitmap_tof_tower->Delete();
    }
    else {
      fout << id << "    " << "-100." << endl;
    }
    delete peak;
    hp_hitmap_tof_tower->Delete();
  }
  fout.close();

  TH1D* hp_hitmap_tof_tower = (TH1D*)h_hitmap_tof_tower->ProjectionY("_py",100,100);
  TSpectrum* peak = new TSpectrum();
  peak->Search(hp_hitmap_tof_tower,2.,"nodraw");
  Float_t peakX = *peak->GetPositionX();
  hp_hitmap_tof_tower->Fit("gaus","","",peakX-3.,peakX+3.);
  TF1* f_hitmap_tof_tower = hp_hitmap_tof_tower->GetFunction("gaus");
  Double_t mean_hitmap_tof_tower = f_hitmap_tof_tower->GetParameter(1);
  cout << "peakX=" << peakX << endl;
  cout << "ToF mean=" << mean_hitmap_tof_tower << endl;

  c->Print("ToF_Calib_Tower.pdf");

}

TGraph* autogain152(TH1 *hist) {

   hist->GetXaxis()->SetRangeUser(200.,16000.);
   TSpectrum *s = new TSpectrum();
   Int_t nfound = s->Search(hist,6,"",0.08); //This will be dependent on the source used.
   printf("Found %d candidate peaks to fit\n",nfound);
   if(nfound > 6)
      nfound = 6;

   std::vector<float> vec;
   for(int x=0;x<nfound;x++)
      vec.push_back(s->GetPositionX()[x]);

   std::sort(vec.begin(),vec.end());

   Float_t energies[] = {121.7830, 244.6920, 344.276, 778.903, 964.131, 1408.011};
   TGraph* slopefit = new TGraph(nfound, &(vec[0]), energies);

   printf("Now fitting: Be patient\n");
   slopefit->Fit("pol1");
   if(slopefit->GetFunction("pol1")->GetChisquare() > 10.) {
      slopefit->RemovePoint(slopefit->GetN()-1);
      slopefit->Fit("pol1");
   }
   TChannel *chan = 0;
   slopefit->Draw("AC*");

   return slopefit;
}

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

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");
}

TGraph* autogain(TH1 *hist,TNucleus *nuc) {    //Display The fits on a TPad  

   if(!hist || !nuc)
      return 0;

   nuc->SetSourceData();

   if(nuc->GetA() == 152) {
      return autogain152(hist);
   }

// Search
   hist->GetXaxis()->SetRangeUser(200.,16000.);
   TSpectrum *s = new TSpectrum();
   Int_t nfound = s->Search(hist,6,"",0.1); //This will be dependent on the source used.
   printf("Found %d candidate peaks to fit\n",nfound);
// Match

   nuc->TransitionList.Sort();

   std::vector<float> engvec;
   TIter iter(&(nuc->TransitionList));
   TObject* obj;
   while(obj = iter.Next()) {
      if(!obj->InheritsFrom("TGRSITransition"))
         continue;
      TGRSITransition *tran = (TGRSITransition*)obj;

      engvec.push_back(static_cast<float>(tran->energy));
      if(engvec.size() == nfound)
         break;
   }

   if(nfound != engvec.size())
      return 0;

   Float_t *posPeaks = s->GetPositionX();
   Float_t *energies = &(engvec[0]);

   for(int x=0;x<nfound;x++) {
      printf("posPeaks[%i] = %f\t\tenrgies[%i] = %f\n",x,posPeaks[x],x,energies[x]);
   }

   TGraph *slopefit = new TGraph(nfound,posPeaks,energies ); 

   printf("Now fitting: Be patient\n");
   slopefit->Fit("pol1");
   slopefit->Draw("AC*");

   return slopefit;

}

//______________________________________________________________________________
void Fit(Int_t run)
{
    // Perform fit.
    
    Char_t tmp[256];
    
    // delete old function
    if (gFitFunc) delete gFitFunc;

    // create fitting function
    if (gFitFunc) delete gFitFunc;
    sprintf(tmp, "fGauss_%d", run);
    gFitFunc = new TF1(tmp, "expo(0)+gaus(2)");
    gFitFunc->SetLineColor(2);
    
    // estimate peak position
    TSpectrum s;
    s.Search(gH, 10, "goff nobackground", 0.05);
    Double_t peak = TMath::MaxElement(s.GetNPeaks(), s.GetPositionX());
    
    // prepare fitting function
    gFitFunc->SetRange(gMin, gMax);
    gFitFunc->SetParameter(2, gH->GetXaxis()->FindBin(peak));
    gFitFunc->SetParLimits(2, 0, 100000);
    gFitFunc->SetParameter(3, peak);
    gFitFunc->SetParameter(4, 20);
    gFitFunc->SetParLimits(4, 18, 100);
    
    for (Int_t i = 0; i < 10; i++)
        if (!gH->Fit(gFitFunc, "RBQ0")) break;
    
    // get peak
    peak = gFitFunc->GetParameter(3);

    // indicator line
    gLine->SetX1(peak);
    gLine->SetX2(peak);
    gLine->SetY1(0);
    gLine->SetY2(gH->GetMaximum());

    // draw 
    gCFit->cd();
    gH->GetXaxis()->SetRangeUser(0, 2000);
    gH->Draw();
    gFitFunc->Draw("same");
    gLine->Draw("same");

    // fill overview histogram
    gHOverview->SetBinContent(run+1, peak);
    gHOverview->SetBinError(run+1, 0.0001);
}

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

/*
  Finds peak in timing difference spectrum
*/
double DiffPeak(TTree *&tree, TString pidcut = "PID>1000")
{
	TSpectrum *spec = new TSpectrum(10); // peak finder
	TH1D *diff = new TH1D("diff","diff",200,-5,5);

	tree->Project("diff","diff",pidcut.Data());
	spec->Search(diff,2,"nodraw",0.7);
	double diffpeak = spec->GetPositionX()[0];

	// testing
	//TCanvas *c_tmp = new TCanvas();
	//diff->Draw();
	//c_tmp->WaitPrimitive();
	//delete c_tmp;
	

	delete diff;
	delete spec;

	return diffpeak;
} // end DiffPeak()

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);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  // Define fit ranges
  double xLowBin[nPMT][nBinsXY][nBinsXY];
  double xHighBin[nPMT][nBinsXY][nBinsXY];
  double xLow[nPMT][nBinsXY][nBinsXY];
  double xHigh[nPMT][nBinsXY][nBinsXY];
  int maxBin[nPMT][nBinsXY][nBinsXY];
  double maxCounts[nPMT][nBinsXY][nBinsXY];
  double binCenterMax[nPMT][nBinsXY][nBinsXY];
  double meanVal[nPMT][nBinsXY][nBinsXY]; // Holds the mean of the distribution as defined by first fitting the peak
  double fitMean[nPMT][nBinsXY][nBinsXY]; // Holds the mean of the low energy peak
  double fitSigma[nPMT][nBinsXY][nBinsXY]; // Holds the sigma of the low energy peak
  double endpoint[nPMT][nBinsXY][nBinsXY]; // Holds the endpoint


  /////// First determine roughly where the low energy peak is
  TSpectrum *spec;

  for (int p=0; p<nPMT; p++) {
    for (int i=0; i<nBinsXY; i++) {
      for (int j=0; j<nBinsXY; j++) {	
	
	double r = sqrt(power(posmap.getBinCenter(j),2)+power(posmap.getBinCenter(i),2));
	
        // Find bin with maximum content
	hisxy[p][i][j]->GetXaxis()->SetRange(2,nBinHist);
        maxBin[p][i][j] = hisxy[p][i][j]->GetMaximumBin();
        maxCounts[p][i][j] = hisxy[p][i][j]->GetBinContent(maxBin[p][i][j]);
        binCenterMax[p][i][j] = hisxy[p][i][j]->GetBinCenter(maxBin[p][i][j]);
	
	  
	if (r<=(50.+2*xyBinWidth))
	      {
		spec = new TSpectrum(20);
		Int_t npeaks = spec->Search(hisxy[p][i][j],1.5,"",0.5);
		
		if (npeaks==0)
		  {
		    cout << "No peaks identified at PMT" << p << " position " << posmap.getBinCenter(i) << ", " << posmap.getBinCenter(j) << endl;
		  }
		else
		  {
		    Float_t *xpeaks = spec->GetPositionX(); // Note that newer versions of ROOT return a pointer to double...
		    TAxis *xaxis = (TAxis*)(hisxy[p][i][j]->GetXaxis());
		    Int_t peakBin=0;
		    Double_t BinSum=0.;
		    Double_t BinSumHold = 0.;
		    Int_t maxPeak=0.;
		    for (int pk=0;pk<npeaks;pk++) {
		      peakBin = xaxis->FindBin(xpeaks[pk]);
		      //Sum over 3 center bins of peak and compare to previos BinSum to see which peak is higher
		      BinSum = hisxy[p][i][j]->GetBinContent(peakBin) + hisxy[p][i][j]->GetBinContent(peakBin-1) + hisxy[p][i][j]->GetBinContent(peakBin+1);
		      if (BinSum>BinSumHold) {
			BinSumHold=BinSum;
			maxPeak=pk;
		      }
		    }
		    binCenterMax[p][i][j] = xpeaks[maxPeak];
		  }
		delete spec;
	      }
	xLow[p][i][j] = binCenterMax[p][i][j]*2./3.;
	xHigh[p][i][j] = 1.5*binCenterMax[p][i][j];

	
      }
    }

//------------------------------------------------------------------------
void DrawWalk()
{
  Int_t npeaks = 20;
  Int_t sigma=3.;
  Bool_t down = false;

  Int_t index[20];
  Char_t buf1[10], buf2[10];
  
  TCanvas *c1 = new TCanvas("c1", "c1",0,48,1280,951);
  gStyle->SetOptStat(0);
  c1->Divide(4,3);
  
  for (Int_t i=0; i<12; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"T0_C_%i_CFD",i+1);
      sprintf(buf2,"CFDvsQTC%i",i+1);
      cout<<buf1<<" "<<buf2<<endl;
      TH2F *qtc_cfd = (TH2F*) gFile->Get(buf2);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      //       cfd->Draw();
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.05);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0){
	Float_t *xpeak = s->GetPositionX();
  	TMath::Sort(nfound, xpeak, index,down);
  	Float_t xp = xpeak[index[0]];
	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
	Float_t yp = cfd->GetBinContent(xbin);
	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
	Float_t hmax = xp+10*sigma;
	Float_t hmin = xp-10*sigma;
	cout<<hmin<< " "<<hmax<<endl;
	//	    cfd->GetXaxis()->SetRange(hmin,hmax);
	//	    TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
	// cfd->Fit("g1","R");
	Int_t hminbin=  qtc_cfd->GetXaxis()->GetFirst();
	Int_t hmaxbin=  qtc_cfd->GetXaxis()->GetLast();
	Int_t nbins= qtc_cfd->GetXaxis()->GetNbins();
	cout<<"  qtc_cfd "<<hminbin<<" "<<hmaxbin<<" "<<nbins<<endl;
	//  qtc_cfd->Draw();
	TProfile *pr_y = qtc_cfd->ProfileX();
	Float_t maxHr=pr_y->GetBinCenter(hmaxbin);
	
	pr_y->SetMaximum(hmax);
	pr_y->SetMinimum(hmin);
	Int_t np=nbins/20;
	Double_t *xx = new Double_t[np];
	Double_t *yy = new Double_t[np];
	Int_t ng=0;
	Double_t yg=0;
	for (Int_t ip=1; ip<nbins; ip++)
	  {
	    if(ip%20 != 0 ) {
	      if (pr_y->GetBinContent(ip) !=0)
		yg +=pr_y->GetBinContent(ip);
	      //	cout<<ng<<" "<<pr_y->GetBinContent(ip)<<" "<<yg<<endl;
	      ng++;}
	    else {
	      xx[ip/20] = Float_t (pr_y->GetBinCenter(ip));
	      yy[ip/20] = yg/ng;
	      yg=0;
	      ng=0;
	     	      cout<<ip<<" "<<ip/20<<" "<< xx[ip/20]<<" "<< yy[ip/20]<<endl;
	    }
	  }
	TH2F *hr = new TH2F("hr"," ",np,0,maxHr, np, hmin, hmax);
	hr->Draw();
	TGraph *gr = new TGraph(np,xx,yy);
	gr->SetMinimum(hmin);
	gr->SetMaximum(hmax);
	gr->SetMarkerStyle(20);
	gr->Draw("P");
	//	delete [] xx;
	//	delete [] yy;
	
	// pr_y->Rebin(10);
		   //  pr_y->Draw();
	    
      }
      
    }
  

}

//--------------------------------------
//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) {
//.........这里部分代码省略.........

//________________________________________________________________________________
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 += "},";
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  for(Int_t i=0;i<8;i++) {
    sprintf(canvastitle,"%s%i","Raw_FEBEX_300_",i);  	
    canv_300[i] = new TCanvas(canvastitle,canvastitle,0,0,1220,900);
    canv_300[i]->SetFillColor(0);
    canv_300[i]->SetFrameFillColor(0);
    canv_300[i]->Divide(4,4);
  }
  
  TF1** myfit_30 = new TF1*[512];
  TF1** myfit_300 = new TF1*[512];
 
  Int_t nfound;
  TSpectrum *s = new TSpectrum(2*npeaks);
  
  //TH1F** hEnergy2_30 = new TH1F*[128];
  //TH1F** hEnergy2_300 = new TH1F*[128];
  TF1 *fline = new TF1("fline","pol1",0,1000);
  Float_t *xpeaks;
  Float_t xp=0;
  Int_t bin=0;
  Float_t yp=0;
  
  for(Int_t i=0;i<128;i++) { 
    //sprintf(canvastitle,"%s%i","hEnergy2_30_",i);  	   
    //TH1F *hEnergy2_30 = (TH1F*)hEnergy_30[i]->Clone(canvastitle);
    //sprintf(canvastitle,"%s%i","hEnergy2_300_",i);  	   
    //TH1F *hEnergy2_300 = (TH1F*)hEnergy_300[i]->Clone(canvastitle);

    if(i<16) canv_30[0]->cd(i+1);
    else if(i>15 &&i<32) canv_30[1]->cd(i-15);
    else if(i>31 &&i<48) canv_30[2]->cd(i-31);
    else if(i>47 &&i<64) canv_30[3]->cd(i-47);
    else if(i>63 &&i<80) canv_30[4]->cd(i-63);
    else if(i>79 &&i<96) canv_30[5]->cd(i-79);
    else if(i>95 &&i<112) canv_30[6]->cd(i-95);
    else if(i>111 &&i<128) canv_30[7]->cd(i-111);
    
    hEnergy_30[i]->Draw();
    //Use TSpectrum to find the peak candidates
    nfound = s->Search(hEnergy_30[i],20,"new");
    xpeaks = s->GetPositionX();
    for (Int_t poo=0;poo<nfound;poo++) {
      Float_t xp = xpeaks[poo];
      cout << " pos[" << poo << "]=" << xp; 
      sprintf(canvastitle,"%s%i%s%i","myfit_30_",i,"_",poo);  	
      myfit_30[i*poo] = new TF1(canvastitle,"[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2]) + [3]",xpeaks[poo]-25,xpeaks[poo]+25);
      myfit_30[i*poo]->SetParameter(0, 500);
      myfit_30[i*poo]->SetParameter(1, xpeaks[poo]);
      myfit_30[i*poo]->SetParameter(2, 20);
      myfit_30[i*poo]->SetParameter(3, 0);  
      myfit_30[i*poo]->SetLineColor(2);
      myfit_30[i*poo]->SetLineWidth(1);
      hEnergy_30[i]->Fit(canvastitle,"R");  
      fitParFile << myfit_30[i*poo]->GetParameter(0) << " " 
		 << myfit_30[i*poo]->GetParameter(1) << " "
		 << myfit_30[i*poo]->GetParameter(2) << " "
		 << myfit_30[i*poo]->GetParameter(3) << endl;
    }
    cout << endl;
    
    if(i<16) canv_300[0]->cd(i+1);
    else if(i>15 &&i<32) canv_300[1]->cd(i-15);
    else if(i>31 &&i<48) canv_300[2]->cd(i-31);
    else if(i>47 &&i<64) canv_300[3]->cd(i-47);
    else if(i>63 &&i<80) canv_300[4]->cd(i-63);
    else if(i>79 &&i<96) canv_300[5]->cd(i-79);
    else if(i>95 &&i<112) canv_300[6]->cd(i-95);
    else if(i>111 &&i<128) canv_300[7]->cd(i-111);
    
    hEnergy_300[i]->Draw();
    //Use TSpectrum to find the peak candidates
    nfound = s->Search(hEnergy_300[i],20,"new");
    xpeaks = s->GetPositionX();
    for (poo=0;poo<nfound;poo++) {
      Float_t xp = xpeaks[poo];
      cout << " pos[" << poo << "]=" << xp; 
      sprintf(canvastitle,"%s%i%s%i","myfit_300_",i,"_",poo);  	
      myfit_300[i*poo] = new TF1(canvastitle,"[0] / sqrt(2.0 * TMath::Pi()) / [2] * exp(-(x-[1])*(x-[1])/2./[2]/[2]) + [3]",xpeaks[poo]-20,xpeaks[poo]+20);
      myfit_300[i*poo]->SetParameter(0, 500);
      myfit_300[i*poo]->SetParameter(1, xpeaks[poo]);
      myfit_300[i*poo]->SetParameter(2, 20);
      myfit_300[i*poo]->SetParameter(3, 0);  
      myfit_300[i*poo]->SetLineColor(2);
      myfit_300[i*poo]->SetLineWidth(1);
      hEnergy_300[i]->Fit(canvastitle,"R");  
      fitParFile << myfit_300[i*poo]->GetParameter(0) << " " 
		 << myfit_300[i*poo]->GetParameter(1) << " "
		 << myfit_300[i*poo]->GetParameter(2) << " "
		 << myfit_300[i*poo]->GetParameter(3) << endl;
    }
    cout << endl;
    
  }
  for(Int_t i=0;i<8;i++)   {
    canv_30[i]->Update();
    canv_300[i]->Update();
  }
  fitParFile.close();
  
}