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

//------------------------------------------------------------------------------
std::pair<double,double> fitskwlin( char* hs, double xl=-0.1, double xr=0.1 ) {
  TH1 *h = (TH1*)gDirectory->Get(hs);

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

  int nb = h->GetNbinsX();
  double x1 = h->GetBinCenter(1); // first
  double x9 = h->GetBinCenter(nb); // last

  if( xl > x1 && xl < x9 ) x1 = xl; // left
  if( xr > x1 && xr < x9 ) x9 = xr; // right

  // create a TF1 with the range from x1 to x9 and 3 parameters
  TF1 *tanhFcn = new TF1( "tanhFcn", fitSkw, x1, x9, 2 );
  tanhFcn->SetParName( 0, "dyoff" );
  tanhFcn->SetParName( 1, "dyslp" );

  // set start values:
  tanhFcn->SetParameter( 0, 0 ); // dy off [um]
  tanhFcn->SetParameter( 1, 99 ); // dy slope [um/skw]

  h->Fit( "tanhFcn", "R Q", "p" );// R = range from tanhFcn

  return std::make_pair(tanhFcn->GetParameter(0),tanhFcn->GetParameter(1));
}

Double_t fitgp0( char* hs ) {

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

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

  h->SetMarkerStyle(21);
  h->SetMarkerSize(0.8);
  h->SetStats(1);
  gStyle->SetOptFit(101);

  gROOT->ForceStyle();

  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 4 parameters
  TF1 *gp0Fcn = new TF1( "gp0Fcn", gp0Fit, x1, x9, 4 );

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

  gp0Fcn->SetNpx(500);
  gp0Fcn->SetLineWidth(4);
  gp0Fcn->SetLineColor(kMagenta);
  gp0Fcn->SetLineColor(kGreen);

  // set start values for some parameters:
  gp0Fcn->SetParameter( 0, xmax ); // peak position
  gp0Fcn->SetParameter( 1, 4*dx ); // width
  gp0Fcn->SetParameter( 2, nn ); // N
  gp0Fcn->SetParameter( 3, bg );

  // N: not drawing
  // Q: quiet
  // R: use specified range
  h->Fit( "gp0Fcn", "NQR", "ep" );

  return gp0Fcn->GetParameter(1);

}

void KVCanvas::ProjectionX(TH2* hh)
{
   TString pname = Form("%s_px", hh->GetName());
   Int_t ip = 1;
   while (gROOT->FindObject(pname.Data())) {
      pname = Form("%s_px%d", hh->GetName(), ip);
      ip++;
   }

   TH1* px = hh->ProjectionX(pname.Data());
   if (!px) return;
   Double_t minY = (hh->GetYaxis()->GetXmin());
   Double_t maxY = (hh->GetYaxis()->GetXmax());
   Double_t dY = (maxY - minY) * 0.8;

   Double_t maxH = px->GetBinContent(px->GetMaximumBin());

   TGraph* gg = 0;
   if ((gg = (TGraph*)gROOT->FindObject(Form("%s_gjx", hh->GetName())))) gg->Delete();

   gg = new TGraph;
   for (int i = 0; i < px->GetNbinsX(); i++) {
      gg->SetPoint(i, px->GetBinCenter(i), minY + px->GetBinContent(i)*dY / maxH);
   }

   gg->SetName(Form("%s_gjx", hh->GetName()));
   gg->SetTitle(Form("%s_gjx", hh->GetName()));
   gg->SetLineColor(kBlack);
   gg->SetMarkerColor(kBlack);
   gg->SetMarkerStyle(8);
   gg->Draw("PL");

   Modified();
   Update();
}

RooHistN::RooHistN(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the asymmetry between the specified TH1 objects
  // which may have fixed or variable bin widths, but which must both have
  // the same binning. The asymmetry is calculated as (1-2)/(1+2). Error bars are
  // calculated using Binomial statistics. Prints a warning and rounds
  // any bins with non-integer contents. Use the optional parameter to
  // specify the confidence level in units of sigma to use for
  // calculating error bars. The nominal bin width specifies the
  // default used by addAsymmetryBin(), and is used to set the relative
  // normalization of bins with different widths. If not set, the
  // nominal bin width is calculated as range/nbins.

  initialize();
  // copy the first input histogram's name and title
  SetName(data1.GetName());
  SetTitle(data1.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data1).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
		     data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
  // initialize our contents from the input histogram contents
  Int_t nbin= data1.GetNbinsX();
  if(data2.GetNbinsX() != nbin) {
    coutE(InputArguments) << "RooHistN::RooHistN: histograms have different number of bins" << endl;
    return;
  }
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data1.GetBinCenter(bin);
    if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
      coutW(InputArguments) << "RooHistN::RooHistN: histograms have different centers for bin " << bin << endl;
    }
    Stat_t y1= data1.GetBinContent(bin);
    Stat_t y2= data2.GetBinContent(bin);
    addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac);
  }
  // we do not have a meaningful number of entries
  _entries= -1;
}

  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

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

std::vector<double> fitskwpol( char* hs, double xl=-0.1, double xr=0.1) {
  TH1 *h = (TH1*)gDirectory->Get(hs);

  if( h == NULL ) {
    cout << hs << " does not exist\n";
    return std::vector<double>();
  }

  int nb = h->GetNbinsX();
  double x1 = h->GetBinCenter(1); // first
  double x9 = h->GetBinCenter(nb); // last

  if( xl > x1 && xl < x9 ) x1 = xl; // left
  if( xr > x1 && xr < x9 ) x9 = xr; // right

  // create a TF1 with the range from x1 to x9 and 3 parameters
  TF1 *tanhFcn = new TF1( "tanhFcn", fitSkwPol, x1, x9, 4 );
  tanhFcn->SetParName( 0, "dyoff" );
  tanhFcn->SetParName( 1, "dyslp" );
  tanhFcn->SetParName( 2, "dypar" );
  tanhFcn->SetParName( 3, "dyhyp" );

  // set start values:
  tanhFcn->SetParameter( 0, 0 ); // dy off [um]
  tanhFcn->SetParameter( 1, 99 ); // dy slope [um/skw]
  tanhFcn->SetParameter( 2, 0 ); // dy parabola
  tanhFcn->SetParameter( 3, 0 ); // dy hyperbola


  h->Fit( "tanhFcn", "R Q", "p" );// R = range from tanhFcn

  std::vector<double> result;
  for(size_t i = 0; i < 4; i++) {
    result.push_back(tanhFcn->GetParameter(i));
  }
  return result;
}

Double_t fitfulllang( char* hs ) {

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

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

  double aa = h->GetEntries();//normalization

  // find peak:
  int ipk = h->GetMaximumBin();
  double xpk = h->GetBinCenter(ipk);
  double sm = xpk / 9; // sigma
  double ns = sm; // noise

  // fit range:
  int ib0 = ipk/2;
  int ib9 = h->GetNbinsX() - 1;

  double x0 = h->GetBinLowEdge(ib0);
  double x9 = h->GetBinLowEdge(ib9) + h->GetBinWidth(ib9);

  // create a TF1 with the range from x0 to x9 and 4 parameters
  TF1 *fitFcn = new TF1( "fitFcn", fitLandauGauss, x0, x9, 4 );

  fitFcn->SetParName( 0, "peak" );
  fitFcn->SetParName( 1, "sigma" );
  fitFcn->SetParName( 2, "area" );
  fitFcn->SetParName( 3, "smear" );

  fitFcn->SetNpx(500);
  fitFcn->SetLineWidth(4);
  fitFcn->SetLineColor(kMagenta);

  // set start values:
  fitFcn->SetParameter( 0, xpk ); // peak position, defined above
  fitFcn->SetParameter( 1, sm ); // width
  fitFcn->SetParameter( 2, aa ); // area
  fitFcn->SetParameter( 3, ns ); // noise

  h->Fit("fitFcn", "NQR", "ep" );// R = range from fitFcn
  return fitFcn->GetParameter(0);
}

//======================================================================
// takes the ID of a graph to fill into a pre-booked histo
//
void fill1DHistoFromGraph(std::string& gid,
			  wTH1 *&wth1)
{
  map<string, wGraph_t *>::const_iterator it=glmap_id2graph.find(gid);
  if (it==glmap_id2graph.end()) {
    cerr<<"Couldn't find graph with id "<<gid<<", define first"<<endl;
    exit(-1);
  }

  if (gl_verbose)
    cout << "Loading histo from graph " << gid << endl;

  TH1 *h = wth1->histo();

  for (int ibin=1; ibin <= h->GetNbinsX(); ibin++)
    h->SetBinContent(ibin,it->second->gr->Eval(h->GetBinCenter(ibin)));

}                                               //  fill1DHistoFromGraph

void HistToText () {
	
	TFile	*f = new TFile("/music/bpt1/bpt/bdn/shane/137i/rootfiles/137i07.root", "READ");
//	TFile	*f = new TFile("/music/bpt1/bpt/bdn/shane/135sb/rootfiles/135sb08.root", "READ");
//	TFile	*f = new TFile("/music/bpt1/bpt/bdn/shane/136sb/rootfiles/136sb01.root", "READ");
//	TFile	*f = new TFile("/music/bpt1/bpt/bdn/shane/140i/rootfiles/140i0.root", "READ");
	
	TH1		*h = (TH1*)f->Get("h_En");
	Int_t nBins = h->GetXaxis()->GetNbins();
	cout << nBins << endl;
	
	ofstream outfile;
	outfile.open("HistToTextOutput.txt", std::ofstream::out);
	
	Double_t x, y;
	for (Int_t i = 0; i<=nBins; i++) {
		x = h->GetBinCenter(i);
		y = h->GetBinContent(i);
		outfile << x << "," << y << endl;
	}
	outfile.close();
}

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

    tree->SetBranchAddress("bad",  &bad);
    tree->SetBranchAddress("centrality",   &centrality);
    tree->SetBranchAddress("trk_Q",   trk_Q   );
    tree->SetBranchAddress("trk_Psi", trk_Psi );
    tree->SetBranchAddress("cal_Q",   cal_Q   );
    tree->SetBranchAddress("cal_Psi", cal_Psi );

    float Trk_Q  [RDET][NHAR];
    float Trk_Psi[RDET][NHAR];
    float Cal_Q  [CDET][NHAR];
    float Cal_Psi[CDET][NHAR];

    Init(cent0, cent1);

    cout<<"Run Centrality "<<cent0<<"   "<<cent1<<endl;
    int nevents = tree->GetEntries();
    cout<<"nevents  "<<nevents<<endl;
       // for(int iev=0; iev<2000000; iev++){
    float maxqcut[CDET][NHAR] = {{0}};
        for(int ihar=0; ihar<NHAR; ihar++){ for(int idet=0; idet<CDET; idet++){ maxqcut[idet][ihar] = 0; }}
        for(int iev=0; iev<nevents; iev++){
        tree->GetEntry(iev);
        if(iev%1000000==0) cout<<iev<<endl;

        if(centrality>=cent1 || centrality<cent0) continue;

        hcent->Fill(centrality);
        for(int idet = 0; idet<RDET; idet++){
            for(int ihar=0; ihar<NHAR; ihar++){
                Trk_Q[idet][ihar]   = trk_Q[idet][ihar];
                Trk_Psi[idet][ihar] = trk_Psi[idet][ihar];
            }
        }

        for(int idet = 0; idet<CDET; idet++){
            for(int ihar=0; ihar<NHAR; ihar++){
                Cal_Q[idet][ihar]   = cal_Q[idet][ihar];
                Cal_Psi[idet][ihar] = cal_Psi[idet][ihar];
            }
        }
        double dQ[]   = {0,0,0,0,0,0,0};
        double dPsi[] = {0,0,0,0,0,0,0};

        for(int ih=0; ih<NHAR;ih++){
            dQ[ih]   = Cal_Q[0][ih] - Cal_Q[1][ih];
            dPsi[ih] = (Cal_Psi[0][ih] - Cal_Psi[1][ih]);
            dPsi[ih] = atan2(sin(dPsi[ih]), cos(dPsi[ih]));
        }

        for(int ihar=0; ihar<NHAR; ihar++){
            for(int idet=0; idet<CDET; idet++){
                if(Cal_Q[idet][ihar] >maxqcut[idet][ihar]) maxqcut[idet][ihar] = Cal_Q[idet][ihar];
                hqFcal[idet][ihar] ->Fill(Cal_Q[idet][ihar]);
            }
            hQPsi[NA-1][ihar] ->Fill(dPsi[ihar],dQ[ihar]);
        }

    }//end of events

    TH1* htmp;
    float steps[] = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.01};
    float qcut[CDET][NHAR][NQ];
    for(int idet=0; idet<CDET; idet++){
        for(int ihar=0; ihar<NHAR; ihar++){
            sprintf(name,"scale_%d_%d", idet, ihar);
            htmp = (TH1*)hqFcal[idet][ihar] ->Clone(name);
            htmp-> Scale(1.0/htmp->GetEntries());
            htmp-> SetTitle(name);
            double sum = 0;
            int cnt = 0;
            for(int ib=1; ib<=htmp->GetNbinsX(); ib++){
                sum+= htmp->GetBinContent(ib);
                if(sum>steps[cnt]) { qcut[idet][ihar][cnt] = htmp->GetBinCenter(ib); cnt++;}
            }
            qcut[idet][ihar][9] = maxqcut[idet][ihar] + 0.00001;
            cout<<idet<<"  "<<ihar<<"  "<<cnt<<"   "<<sum<<endl;
        }
    }

    ofstream tout;
    sprintf(name,"cutInfo1_%d_%d.txt", cent0, cent1);
    tout.open(name);
   tout<<"cent "<<cent0<<"--"<<cent1<<endl;
    for(int idet=0; idet<CDET; idet++){
        tout<<"IDET = "<<idet<<endl;
        for(int ihar=0; ihar<NHAR; ihar++){
            for(int iq=0; iq<NQ; iq++){
                if(iq==0) tout<<"{"<<qcut[idet][ihar][iq]<<", ";
                else if(iq<9) tout<<qcut[idet][ihar][iq]<<", ";
                else if(iq ==9) tout<<qcut[idet][ihar][iq]<<"}, "<<endl;
            }
        }
    }
    tout.close();
    fileout->Write();
    fileout->Close();

    }

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

  trend_canvas->Print(trend_openpdfname.c_str());



  // Now we have the histogram names, so loop though them
  for (int iHist = 0; iHist < ordered_histogram_names.size(); ++iHist) {

    std::string histname = ordered_histogram_names.at(iHist);
    std::string fullhistname = "DataQuality_LowLevel/" + histname;

    // Check to see if we actually want this plot
    if ( !WantPlot(histname)) {
      continue;
    }


    // Create the trend plot
    std::string trendplotname = histname + "_TrendPlot";
    std::string trendplottitle = "Trend plot of " + histname;

    // Get some useful information from the histogram in the first run file (for the y-axis range of the trend plot)
    TH1F* hist = (TH1F*) files[0]->Get(fullhistname.c_str());

    TH2F* trend_plot = new TH2F(trendplotname.c_str(), trendplottitle.c_str(), n_runs, first_run, first_run+n_runs, hist->GetNbinsX(), hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());

    trend_plot->GetYaxis()->SetTitle(hist->GetXaxis()->GetTitle());
    trend_plot->GetYaxis()->SetLabelSize(0.03);

    trend_plot->GetXaxis()->SetTitle("Run Number");
    trend_plot->GetXaxis()->SetLabelSize(0.03);

    trend_plot->GetZaxis()->SetTitle(hist->GetYaxis()->GetTitle());
    trend_plot->GetZaxis()->SetLabelSize(0.03);
    trend_plot->GetZaxis()->SetTitleOffset(0.85);
    trend_plot->GetZaxis()->SetTitleSize(0.03);
    trend_plot->SetStats(false);

    // Set the fraction trend plots to have a maximum of 1 so that red = bad
    if (histname.find("Fraction") != std::string::npos) {
      trend_plot->GetZaxis()->SetRangeUser(0,1);
    }
    // Copy the axis labels so that they show the detector and channel names
    else if (histname.find("hDQ_IslandRate") != std::string::npos) {
      for (int jBin = 1; jBin < hist->GetNbinsX(); ++jBin) {
	trend_plot->GetYaxis()->SetBinLabel(jBin, hist->GetXaxis()->GetBinLabel(jBin));
      }
      trend_plot->GetYaxis()->SetTitle(false);
      trend_plot->GetYaxis()->SetLabelSize(0.02);
      trend_plot->GetZaxis()->SetLabelSize(0.02);
    }

    // Check to see if we want this histogram as an individual or a trend plot
    bool want_trend_plot = WantAsTrendPlot(histname);

    // Loop through the runs
    for (int iRun = 0; iRun < n_runs; ++iRun) {
      
      // Make sure this file exists (there may be some missing in our golden run ranges)
      if (files[iRun]->IsZombie())
	continue;

      // Get the histogram
      TH1* hist = (TH1*) files[iRun]->Get(fullhistname.c_str());
      
      if (want_trend_plot == false) {

	// Print the histogram out directly to the individual picture book
	individual_canvases[iRun]->cd();
	ZoomIndividualPlot(hist);
	hist->Draw("COLZ");
	individual_canvases[iRun]->SetLogz(WantLogZ(histname)); // change to a log-z scale if we want to
	individual_canvases[iRun]->Update();
	individual_canvases[iRun]->Print(individual_basepdfname[iRun].str().c_str());
      
      }
      else {
	// Fill this histogram into the trend plot
	for (int iBin = 1; iBin <= hist->GetNbinsX(); ++iBin) {
	  trend_plot->Fill(first_run + iRun, hist->GetBinCenter(iBin), hist->GetBinContent(iBin)); // (x = run #, y = time stamp, z = N_TPI)
	} // end for loop (filling trend plot)
      } // end if (trend plot)
    } // end for loop (through runs)

    // Now export the trend plot to PDF (if applicable)
    if (want_trend_plot) {
      trend_canvas->cd();
      ZoomTrendPlot(trend_plot);
      trend_plot->Draw("COLZ");
      trend_canvas->Print(trend_basepdfname.str().c_str());
    }
  }


  // Now go through and clean up
  for (int iRun = 0; iRun < n_runs; ++iRun) {
    individual_canvases[iRun]->Print(individual_closepdfname[iRun].c_str());
  }
  trend_canvas->Print(trend_closepdfname.c_str());
}

void printCalibStat(Int_t run, const char * fname,  TTreeSRedirector * pcstream){

  //
  // Dump the statistical information about all histograms in the calibration files 
  //    into the statistical tree, print on the screen (log files) as well 
  //
  //
  // 1. Default dump for all histograms
  //    Information to dump:
  //    stat =Entries, Mean, MeanError,  RMS, MaxBin
  //    Branch naming convention:
  //    <detName>_<hisName><statName>
  //
  // 2. Detector statistical information  - to be implemented by expert
  //                                      - First version implemented by MI 
  //  
  // 

  TFile *fin = TFile::Open(fname);
  if (!fin) return;
  const Double_t kMaxHis=10000;
  
  TList * keyList = fin->GetListOfKeys();
  Int_t nkeys=keyList->GetEntries();
  Double_t *hisEntries = new Double_t[kMaxHis];
  Double_t *hisMean = new Double_t[kMaxHis];
  Double_t *hisMeanError = new Double_t[kMaxHis];
  Double_t *hisRMS = new Double_t[kMaxHis];
  Double_t *hisMaxBin = new Double_t[kMaxHis];
  Int_t counter=0;
  
  if (pcstream) (*pcstream)<<"calibStatAll"<<"run="<<run;
  for (Int_t ikey=0; ikey<nkeys; ikey++){
    TObject * object = fin->Get(keyList->At(ikey)->GetName());
    if (!object) continue;
    if (object->InheritsFrom("TCollection")==0) continue;
    TSeqCollection *collection  = (TSeqCollection*)object; 
    Int_t nentries= collection->GetEntries();
    for (Int_t ihis=0; ihis<nentries; ihis++){
      TObject * ohis = collection->At(ihis);
      if (!ohis) continue;
      if (ohis->InheritsFrom("TH1")==0) continue;
      TH1* phis = (TH1*)ohis;
      hisEntries[counter]=phis->GetEntries();	
      Int_t idim=1;
      if (ohis->InheritsFrom("TH2")) idim=2;
      if (ohis->InheritsFrom("TH3")) idim=3;        
      hisMean[counter]=phis->GetMean(idim);	
      hisMeanError[counter]=phis->GetMeanError(idim);	
      hisRMS[counter]=phis->GetRMS(idim);	
      hisMaxBin[counter]=phis->GetBinCenter(phis->GetMaximumBin());	
      if (pcstream) (*pcstream)<<"calibStatAll"<<
		      Form("%s_%sEntries=",keyList->At(ikey)->GetName(), phis->GetName())<<hisEntries[counter]<<	
		      Form("%s_%sMean=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMean[counter]<<	
		      Form("%s_%sMeanError=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMeanError[counter]<<	
		      Form("%s_%sRMS=",keyList->At(ikey)->GetName(), phis->GetName())<<hisRMS[counter]<<	
		      Form("%s_%sMaxBin=",keyList->At(ikey)->GetName(), phis->GetName())<<hisMaxBin[counter];	
      //printf("Histo:\t%s_%s\t%f\t%d\n",keyList->At(ikey)->GetName(), phis->GetName(), hisEntries[counter],idim);
      counter++;
    }
    delete object;
  }    
  
  //
  // Expert dump example (MI first iteration):
  //
  // 0.)  TOF dump
  //

  Int_t tofEvents=0;
  Int_t tofTracks=0;
  TList * TOFCalib = (TList*)fin->Get("TOFHistos");      
  if (TOFCalib) {
    TH1 *histoEvents = (TH1*)TOFCalib->FindObject("hHistoVertexTimestamp");
    TH1 *histoTracks = (TH1*)TOFCalib->FindObject("hHistoDeltatTimestamp");
    if (histoEvents && histoTracks){
      tofEvents = TMath::Nint(histoEvents->GetEntries());
      tofTracks = TMath::Nint(histoTracks->GetEntries());
    }
    delete TOFCalib;
  }
  printf("Monalisa TOFevents\t%d\n",tofEvents);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFevents="<<tofEvents;
  printf("Monalisa TOFtracks\t%d\n",tofTracks);
  if (pcstream) (*pcstream)<<"calibStatAll"<<"TOFtracks="<<tofTracks;

  //
  // 1.)  TPC  dump - usefull events/tracks  for the calibration
  //
  Int_t tpcEvents=0;
  Int_t tpcTracks=0;
  TObject* obj = dynamic_cast<TObject*>(fin->Get("TPCCalib"));
  TObjArray* array = dynamic_cast<TObjArray*>(obj);
  TDirectory* dir = dynamic_cast<TDirectory*>(obj);
  AliTPCcalibTime  * calibTime = NULL;
  if (dir) {
    calibTime = dynamic_cast<AliTPCcalibTime*>(dir->Get("calibTime"));
  }
  else if (array){
    calibTime = (AliTPCcalibTime *)array->FindObject("calibTime");
//.........这里部分代码省略.........

TH1* GetCentK(TDirectory* top, Double_t c1, Double_t c2, Int_t s,
	      TLegend* l)
{
  TString dname; dname.Form("cent%06.2f_%06.2f", c1, c2);
  dname.ReplaceAll(".", "d");
  TDirectory* d = top->GetDirectory(dname);
  if (!d) {
    Warning("GetCetnK", "Directory %s not found in %s",
	    dname.Data(), top->GetName());
    return;
  }

  TDirectory* det = d->GetDirectory("details");
  if (!det) {
    Warning("GetCetnK", "Directory details not found in %s",
	    d->GetName());
    d->ls();
    return;
  }

  TObject* o = det->Get("scalar");
  if (!o) {
    Warning("GetCetnK", "Object scalar not found in %s",
	    det->GetName());
    return;
  }

  if (!o->IsA()->InheritsFrom(TH1::Class())) {
    Warning("GetCetnK", "Object %s is not a TH1, but a %s",
	    o->GetName(), o->ClassName());
    return;
  }
  TH1* h = static_cast<TH1*>(o->Clone());
  Color_t col = cc[(s-1)%10];
  h->SetLineColor(col);
  h->SetMarkerColor(col);
  h->SetFillColor(col);
  h->SetFillStyle(1001);
  // h->SetTitle(Form("%5.2f-%5.2f%% #times %d", c1, c2, s));
  h->SetTitle(Form("%2.0f-%2.0f%% + %d", c1, c2, s-1));
  TF1* f = new TF1("", "[0]",-2.2,2.2);
  f->SetParameter(0,s-1);
  f->SetLineColor(col);
  f->SetLineStyle(7);
  f->SetLineWidth(1);
  // h->Scale(s);
  h->Add(f);
  h->GetListOfFunctions()->Add(f);
  f->SetParameter(0,s);
  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
    if (TMath::Abs(h->GetBinCenter(i)) > 2) {
      h->SetBinContent(i,0);
      h->SetBinError(i,0);
    }
  }
  
  TLegendEntry* e = l->AddEntry(h, h->GetTitle(), "f");
  e->SetFillColor(col);
  e->SetFillStyle(1001);
  e->SetLineColor(col);

  return h;
}

double leptonic_fitter_algebraic::likeliest_scale( const TH1& ITF )
{
  Int_t ibin = ITF.GetMaximumBin();
  return ITF.GetBinCenter( ibin );
}