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

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

void setRange(RooWorkspace& myws, RooPlot* frame, string dsName, bool setLogScale, double dMuonYmin)
{ 
  // Find maximum and minimum points of Plot to rescale Y axis
  TH1* h = myws.data(dsName.c_str())->createHistogram("hist", *myws.var("ctau"), Binning(frame->GetNbinsX(),frame->GetXaxis()->GetXmin(),frame->GetXaxis()->GetXmax()));
  Double_t YMax = h->GetBinContent(h->GetMaximumBin());
  cout << YMax << endl;
  // Double_t YMin = min( h->GetBinContent(h->FindFirstBinAbove(0.0)), h->GetBinContent(h->FindLastBinAbove(0.0)) );
  Double_t YMin = 1e99;
  for (int i=1; i<=h->GetNbinsX(); i++) if (h->GetBinContent(i)>0) YMin = min(YMin, h->GetBinContent(i));
  
  bool isMC = false;
  if (dsName.find("MC")!=std::string::npos) isMC = true;
    
  Double_t Yup(0.),Ydown(0.);
  if(setLogScale)
  {
    if (isMC) Ydown = YMin*0.3;
    else Ydown = YMin/(TMath::Power((YMax/YMin), (0.1/(1.0-0.1-0.4))));
    Yup = YMax*TMath::Power((YMax/YMin), (0.4/(1.0-0.1-0.4)));
  }
  else
  {
    Ydown = max(YMin-(YMax-YMin)*(0.1/(1.0-0.1-0.4)),0.0);
    Yup = YMax+(YMax-YMin)*(0.4/(1.0-0.1-0.4));
  }
  cout << Ydown << " " << Yup << endl;
  frame->GetYaxis()->SetRangeUser(Ydown,Yup);
  delete h;
  
  // Create line to indicate upper fitting range for MC
  if (isMC)
  {
    if (dsName.find("JPSIP")!=std::string::npos)
    {
      TLine* line(0x0);
      if (dMuonYmin >= 1.6) line = new TLine(3.32,Ydown,3.32,Yup);
      else line = new TLine(3.26,Ydown,3.26,Yup);
      line->SetLineStyle(2);
      line->SetLineColor(1);
      line->SetLineWidth(3);
      
      frame->addObject(line);
    }
    else if (dsName.find("PSI2S")!=std::string::npos)
    {
      TLine* line(0x0);
      if (dMuonYmin >= 1.6) line = new TLine(3.95,Ydown,3.95,Yup);
      else line = new TLine(3.85,Ydown,3.85,Yup);
      line->SetLineStyle(2);
      line->SetLineColor(1);
      line->SetLineWidth(3);
      
      frame->addObject(line);
    }

  }
 
}

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

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

// Returns the RMS including 96% of the histogram entries, cutting the tails:
Double_t getRMS96(char* hs, double truncation=96.) {

  bool debug = false;

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

  // Total entries:
  double integral = h->GetEntries();
  int maxbin = h->GetMaximumBin();
  if(debug) cout << "entries=" << integral << " maxbin=" << maxbin << endl;

  double subrange_integral = h->GetBinContent(maxbin);
  int bin = 0;
  while(subrange_integral < truncation/100*integral) {
    bin++;
    // Add one bin to the left:
    subrange_integral += h->GetBinContent(maxbin-bin);
    // Add one bin to the right:
    subrange_integral += h->GetBinContent(maxbin+bin);
    if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << ": entries=" << subrange_integral << endl;
  }
  if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << " now has " << subrange_integral << " entries, this is " << (100.0*subrange_integral)/integral << "%" << endl;

  // Correct by overshoot bin:
  subrange_integral -= h->GetBinContent(maxbin+bin);
  subrange_integral -= h->GetBinContent(maxbin-bin);
  bin--;

  int binlow = maxbin-bin;
  int binhigh = maxbin+bin;
  if(debug) cout << "subrange " << (maxbin-bin) << "-" << (maxbin+bin) << " now has " << subrange_integral << " entries, this is " << (100.0*subrange_integral)/integral << "%" << endl;

  h->GetXaxis()->SetRange(binlow,binhigh); //to restrict range to bins binlow to binhigh
  double rms96 = h->GetRMS(); //will return the RMS within the axis range

  return rms96;
}

//
//----------------------------------------------------------------------
//
int fittp0( char* hs ) {

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

  if( h == NULL ){

    cout << hs << " does not exist\n";

  }

  else{
   
    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);
    cout << hs << ": " << x1 << " - " << x9 << endl;

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

    tp0Fcn->SetNpx(500);
    tp0Fcn->SetLineWidth(4);
    tp0Fcn->SetLineColor(kMagenta);
    tp0Fcn->SetLineColor(kGreen);
   
    // set start values for some parameters:

    cout << hs << " " << dx << ", " << nn << ", " << xmax << endl;

    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", "R", "ep" );
    // h->Fit("tp0Fcn","V+","ep");

    h->Draw("histepsame");  // data again on top
  }
}

//.........这里部分代码省略.........
				TH1* hgen = (TH1*) (outRootFile->Get(gen_hist.str().c_str()));
				//->Clone(gen_hist_name.str().c_str()));
				if (hgen == NULL) { cout << "hgen = " << gen_hist.str() << " was not found!" << endl; assert(false); } 
				hgen->SetDirectory(0);

				hreco->Sumw2();
				hsmear->Sumw2();
				hgen->Sumw2();

				const int rebin = hist2print.at(ihist).Rebin();
				const string title = hist2print.at(ihist).Title();
				const double xmin = hist2print.at(ihist).Xmin();
				const double xmax = hist2print.at(ihist).Xmax();

				if (rebin>1)
				{
					hreco->Rebin(rebin);
					hsmear->Rebin(rebin);
					hgen->Rebin(rebin);
				}
				if (title.length()>0)
				{
					hreco->SetTitle(title.c_str());
					hsmear->SetTitle(title.c_str());
					hgen->SetTitle(title.c_str());
				}
				if (xmin != LargeNegNum || xmax != LargeNegNum)
				{
					hreco->GetXaxis()->SetRangeUser(xmin,xmax);
					hsmear->GetXaxis()->SetRangeUser(xmin,xmax);
					hgen->GetXaxis()->SetRangeUser(xmin,xmax);
				}

				const double reco_max_y  = hreco->GetBinContent(hreco->GetMaximumBin());
				const double smear_max_y = hsmear->GetBinContent(hsmear->GetMaximumBin());
				const double y_max = max(reco_max_y, smear_max_y);
				double y_min = 9999.0;
				for (unsigned bin=1; bin<hreco->GetNbinsX(); ++bin)
				{
					const double v1 = hreco->GetBinContent(bin);
					const double v2 = hsmear->GetBinContent(bin);
					const double minv = min(v1,v2);
					if (minv != 0 && minv < y_min) y_min = minv;
					
				}

				cout << hreco->GetName() << "->ymin/max = " << y_min << "(" << y_min/2.0 << ")/" << y_max << "(" << y_max*2.0 << ")" << endl;
				hreco->GetYaxis()->SetRangeUser(y_min/2.0, y_max*2.0);
				hsmear->GetYaxis()->SetRangeUser(y_min/2.0, y_max*2.0);


				hgen->SetLineColor(kBlue);
				hgen->SetMarkerColor(kBlue);
				hgen->SetMarkerStyle(24);
				hgen->SetLineWidth(2);
				hsmear->SetLineColor(kRed);
				hsmear->SetMarkerColor(kRed);
				hsmear->SetMarkerStyle(24);
				hsmear->SetLineWidth(2);
				hreco->SetLineWidth(2);
				hreco->SetMarkerStyle(kDot);
				hreco->SetLineColor(kBlack);
				hreco->SetMarkerColor(kBlack);
				//hreco->GetXaxis()->SetRangeUser(0,300);
				//hsmear->GetXaxis()->SetRangeUser(0,300);

//______________________________________________________________________________
void CheckTime(const Char_t* loc)
{
    // Check time.

    Char_t t[256];

    // number of runs
    Int_t nRuns = gFiles->GetSize();

    // create arrays
    const Int_t nCh = 352;
    Double_t** pedPos = new Double_t*[nCh];
    Double_t* runNumbersD = new Double_t[nRuns];
    for (Int_t i = 0; i < nCh; i++) pedPos[i] = new Double_t[nRuns];

    // open the output files
    TFile* fROOTout = new TFile("/tmp/tagger_time.root", "RECREATE");

    // create directories
    for (Int_t i = 0; i < nCh; i++)
    {
        sprintf(t, "%03d", i);
        fROOTout->mkdir(t);
    }
    
    TF1* func = new TF1("func", "gaus(0)+pol1(3)", -1000 , 1000);

    // loop over runs
    for (Int_t i = 0; i < nRuns; i++)
    {   
        // get the file
        TFile* f = (TFile*) gFiles->At(i);

        // extract run number
        Int_t runNumber;
        sprintf(t, "%s/ARHistograms_CB_%%d.root", loc);
        sscanf(f->GetName(), t, &runNumber);
        runNumbersD[i] = (Double_t)runNumber;

        printf("Processing run %d (%d/%d)\n", runNumber, i+1, nRuns);

        fROOTout->cd();
        
        TH2* h2 = (TH2*) f->Get("CaLib_Tagger_Time");

        // loop over channels
        for (Int_t j = 0; j < nCh; j++)
        {
            // load histogram
            sprintf(t, "%03d", j);
            fROOTout->cd(t);
            sprintf(t, "%d_%d", i, j);
            TH1* h = h2->ProjectionX(t, j+1, j+1);
            h->Rebin(2);
            
            sprintf(t, "Run_%d", runNumber);
            TCanvas* c = new TCanvas(t, t);
            TLine* tline = 0;
            
            // check entries
            if (h->GetEntries())
            {
                // fit gaussian to peak
                Double_t maxPos = h->GetXaxis()->GetBinCenter(h->GetMaximumBin());
                func->SetParameters(1, maxPos, 0.5, 1, 0.1);
                func->SetRange(maxPos - 2, maxPos + 2);
                func->SetParLimits(0, 0, 1000);
                for (Int_t k = 0; k < 10; k++)
                    if (!h->Fit(func, "RBQ")) break;

                // save position in file and memory
                maxPos = func->GetParameter(1);
                pedPos[j][i] = maxPos;

                h->GetXaxis()->SetRangeUser(maxPos - 10, maxPos + 10);
                h->Draw();
                
                tline = new TLine(maxPos, 0, maxPos, 10000);
                tline->SetLineColor(kRed);
                tline->SetLineWidth(2);
                tline->Draw();
            }
            else
            {
                h->Draw();
            }

            c->Write(c->GetName(), TObject::kOverwrite);
    
            delete h;
            delete c;
            if (tline) delete tline;
        }
        
        delete h2;
    }
    
    // create pedestal evolution graphs
    fROOTout->cd();
//.........这里部分代码省略.........

//________________________________________________________________
void KVZALineFinder::FindALine(Int_t zz, Int_t width)
{
    fLinearHisto->SetAxisRange(zz-0.5,zz+0.5,"Y");

    KVIDLine* line = (KVIDLine*)fGrid->GetIdentifier(zz,2*zz+1); // A=2*zz+1 : dummy, A is ignored in this case
    if(!line)
    {
        int i=1;
        while(!(line = (KVIDLine*)fGrid->GetIdentifier(zz+i,2*zz+1))) i++;
    }
    if(!line) return;

    Double_t lX, lY;
    line->GetStartPoint(lX,lY);
    Int_t xbmin = 0;//fLinearHisto->GetYaxis()->FindBin(lX);
    line->GetEndPoint(lX,lY);
    Int_t xbmax = fLinearHisto->GetXaxis()->FindBin(lX);
    
    // create lines
    TList Lines;
    KVSpiderLine* tmp = 0;

    fLinearHisto->SetAxisRange(fLinearHisto->GetXaxis()->GetBinCenter(50),lX,"X");//fLinearHisto->GetXaxis()->GetXmax(),"X");
    TH1* tmph = fLinearHisto->ProjectionX(Form("tmph%d",zz));
    Int_t startBin = (Int_t)(tmph->GetMaximumBin()*0.95);
    delete tmph;

    TH1* projey = 0;
    if(startBin)
    {
        projey = fLinearHisto->ProjectionY("ProjectionAfterLin",startBin-width*3,startBin+width*3);
        int nfound = fSpectrum.Search(projey,0.05,"goff",0.0001);
        Info("FindALine","%d peack found...",nfound);
#if ROOT_VERSION_CODE > ROOT_VERSION(5,99,01)
        Double_t* xpeaks = fSpectrum.GetPositionX();
        Double_t* ypeaks = fSpectrum.GetPositionY();
#else
        Float_t* xpeaks = fSpectrum.GetPositionX();
        Float_t* ypeaks = fSpectrum.GetPositionY();
#endif
        for(int p=0;p<nfound;p++)
        {
            if(p>8) break;
            if(ypeaks[p]<10) continue;
            Double_t xline = fLinearHisto->GetXaxis()->GetBinCenter(startBin);
            Double_t yline = xpeaks[p];
            KVSpiderLine* tmp = 0;
            TIter next(&Lines);
            while((tmp=(KVSpiderLine*)next()))
            {
                Info("FindALine","line found but I don't know why...");
                if(TMath::Abs(tmp->GetY()-yline)<0.05) continue;
            }
            tmp = new KVSpiderLine(zz,-1);
            Lines.AddLast(tmp);
            tmp->AddPoint(xline,yline);
            fPoints->SetPoint(fNPoints,xline,yline);
            fNPoints++;
        }
        if(projey) delete projey;
    }
    else Error("FindALine","not starting bin indicated...");
    SortLines(&Lines);

    Int_t nLines = Lines.GetSize();
    tmp = 0;
    for(int xx=startBin-width; xx>xbmin; xx-=width)
    {
        projey = fLinearHisto->ProjectionY("ProjectionAfterLin",xx-width/2,xx+width/2);
        int nfound = fSpectrum.Search(projey,0.05,"goff",0.02);
#if ROOT_VERSION_CODE > ROOT_VERSION(5,99,01)
        Double_t* xpeaks = fSpectrum.GetPositionX();
        Double_t* ypeaks = fSpectrum.GetPositionY();
#else
        Float_t* xpeaks = fSpectrum.GetPositionX();
        Float_t* ypeaks = fSpectrum.GetPositionY();
#endif
        for(int p=0;p<nfound;p++)
        {
            if(p>=nLines+1) continue;
            if(ypeaks[p]<5) continue;
            Double_t xline = fLinearHisto->GetXaxis()->GetBinCenter(xx);
            Double_t yline = xpeaks[p];
            KVSpiderLine* tmp = 0;
            TIter next(&Lines);
            while((tmp=(KVSpiderLine*)next()))
            {
                if((TMath::Abs(tmp->GetY()-yline)<0.05)) break;
            }
            if(tmp)
            {
                if((TMath::Abs(tmp->GetX()-xline)<10*width)) tmp->AddPoint(xline,yline);
            }
        }
        if(projey) delete projey;
    }

    TIter nextli(&Lines);
    while((tmp=(KVSpiderLine*)nextli()))tmp->Sort(true);
//.........这里部分代码省略.........

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

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

void ptBestFit(float BIN_SIZE=5.0,bool BLIND=false,TString MASS,TString NAME)
{
  gROOT->ProcessLine(".x ../../common/styleCMSTDR.C");
  gSystem->Load("libHiggsAnalysisCombinedLimit.so");
  gROOT->ForceStyle();
  gStyle->SetOptStat(0);
  gStyle->SetOptTitle(0);
  gROOT->SetBatch(1);
  gStyle->SetPadRightMargin(0.04);
  gStyle->SetPadLeftMargin(0.16);
  gStyle->SetPadTopMargin(0.06);
  gStyle->SetPadBottomMargin(0.10);
  gStyle->SetTitleFont(42,"XY");
  gStyle->SetTitleSize(0.0475,"XY");
  gStyle->SetTitleOffset(0.9,"X");
  gStyle->SetTitleOffset(1.5,"Y");
  gStyle->SetLabelSize(0.0375,"XY");

  RooMsgService::instance().setSilentMode(kTRUE);
  for(int i=0;i<2;i++) {
    RooMsgService::instance().setStreamStatus(i,kFALSE);
  }
  float XMIN = 80;
  float XMAX = 200; 

  TFile *f1 = TFile::Open("datacards/datacard_m"+MASS+"_"+NAME+".root");
  TFile *f2 = TFile::Open("combine/mlfit.vbfHbb_"+NAME+"_mH"+MASS+".root");
  TFile *f3 = TFile::Open("root/sig_shapes_workspace_B80-200.root");
  TFile *f4 = TFile::Open("root/data_shapes_workspace_"+NAME+".root");

  RooWorkspace *w = (RooWorkspace*)f1->Get("w");
  //w->Print();
  RooAbsPdf *bkg_model = (RooAbsPdf*)w->pdf("model_s");
  RooFitResult *res_s  = (RooFitResult*)f2->Get("fit_s"); 
  RooFitResult *res_b  = (RooFitResult*)f2->Get("fit_b");
  RooRealVar *rFit     = dynamic_cast<RooRealVar *>(res_s->floatParsFinal()).find("r");
  RooDataSet *data     = (RooDataSet*)w->data("data_obs");
  
  int nparS=0,nparB=0;
  cout << res_s->floatParsFinal().getSize() << endl;
  cout << res_b->floatParsFinal().getSize() << endl;
  nparS = res_s->floatParsFinal().getSize();
  nparB = res_b->floatParsFinal().getSize();  
  float chi2sumS = 0.;
  float chi2sumB = 0.;
  int nparsum = 0;
//  if (BLIND) {
//    res_b->Print();
//  }
//  else {
//    res_s->Print();
//  }
  
  w->allVars().assignValueOnly(res_s->floatParsFinal());
//  w->Print();
//  w->allVars()->Print();

  RooWorkspace *wSig = (RooWorkspace*)f3->Get("w"); 
  RooWorkspace *wDat = (RooWorkspace*)f4->Get("w"); 

  const RooSimultaneous *sim = dynamic_cast<const RooSimultaneous *> (bkg_model);
  const RooAbsCategoryLValue &cat = (RooAbsCategoryLValue &) sim->indexCat();
  TList *datasets = data->split(cat,true);
  TIter next(datasets);
  //int count = 0; 
  for(RooAbsData *ds = (RooAbsData*)next();ds != 0; ds = (RooAbsData*)next()) {
	 //if (count > 0) return 0;
	 //count++;
    RooAbsPdf *pdfi = sim->getPdf(ds->GetName());
    RooArgSet *obs = (RooArgSet*)pdfi->getObservables(ds);
    RooRealVar *x = dynamic_cast<RooRealVar *>(obs->first());

    RooRealVar *yield_vbf = (RooRealVar*)wSig->var("yield_signalVBF_mass"+MASS+"_"+TString(ds->GetName()));
    RooRealVar *yield_gf  = (RooRealVar*)wSig->var("yield_signalGF_mass"+MASS+"_"+TString(ds->GetName()));
    TString ds_name(ds->GetName());
    //----- get the QCD normalization -----------
    RooRealVar *qcd_norm_final = dynamic_cast<RooRealVar *>(res_s->floatParsFinal()).find("CMS_vbfbb_qcd_norm_"+ds_name);
    RooRealVar *qcd_yield      = (RooRealVar*)wDat->var("yield_data_"+ds_name);

    float Nqcd  = exp(log(1.5)*qcd_norm_final->getVal())*qcd_yield->getVal();
    float eNqcd = log(1.5)*qcd_norm_final->getError()*Nqcd;
    cout<<"QCD normalization = "<<Nqcd<<" +/- "<<eNqcd<<endl;
    
    TH1 *hCoarse = (TH1*)ds->createHistogram("coarseHisto_"+ds_name,*x);
    float norm = hCoarse->Integral();
  
	 int rebin = BIN_SIZE/hCoarse->GetBinWidth(1);
    hCoarse->Rebin(rebin);

    float MIN_VAL = TMath::Max(0.9*hCoarse->GetBinContent(hCoarse->GetMinimumBin()),1.0);
    float MAX_VAL = 1.3*hCoarse->GetBinContent(hCoarse->GetMaximumBin());
    RooDataHist ds_coarse("ds_coarse_"+ds_name,"ds_coarse_"+ds_name,*x,hCoarse);

    TH1F *hBlind = (TH1F*)hCoarse->Clone("blindHisto_"+ds_name);
    for(int i=0;i<hBlind->GetNbinsX();i++) {
      double x0 = hBlind->GetBinCenter(i+1);
      if (x0 > 100 && x0 < 150) {
        hBlind->SetBinContent(i+1,0);
        hBlind->SetBinError(i+1,0);
      }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
       if (xlowVec != 0) hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xlowVec);
       //fixed bin histo
       else hist[theHistCounter] = new TH1F(thisHistName,thisHistName,nBins,xLow,xHigh);
       hist[theHistCounter] -> Sumw2();
       hist[theHistCounter] -> SetFillColor(col[theHistCounter]);
       hist[theHistCounter] -> SetFillStyle(1001);
       theLegends.push_back(*listOfSamples.at(iSample)->Legend());
    }

    //fill the histogram
    TString thisScale = Form("%f *", *(listOfSamples.at(iSample)->Scale()));
    if (isFirstOfSerie) tree[iSample] -> Draw(myVar + " >> " + TString(hist[theHistCounter] -> GetName()),thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    else tree[iSample] -> Draw(myVar + " >>+ " + TString(hist[theHistCounter] -> GetName()),thisScale + "evt_weight*kf_weight*pu_weight" + myCut);
    
    //add the histogram to the stack if the last of the series:
    //either last sample or ~ sample followed by non ~ sample
    if (isLastOfSerie) {
       if (xlowVec != 0) {
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinError  (iBin,hist[theHistCounter]->GetBinError(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
         for (int iBin = 1; iBin <= nBins; iBin++) hist[theHistCounter]->SetBinContent(iBin,hist[theHistCounter]->GetBinContent(iBin)/hist[theHistCounter]->GetBinWidth(iBin));
       }
       hs -> Add(hist[theHistCounter]);
       theHistCounter++;
    }
    
  }//end loop on samples

  //Fix the legend
  for (int iHisto = theHistCounter-1; iHisto >= 0; iHisto--) {
    leg -> AddEntry(hist[iHisto], theLegends[iHisto], "f");   
  }
  
  //get the maximum to properly set the frame
  float theMax = hdata -> GetBinContent(hdata -> GetMaximumBin()) + hdata -> GetBinError(hdata -> GetMaximumBin());
  TH1* theMCSum = (TH1*) hs->GetStack()->Last();
  float theMaxMC = theMCSum->GetBinContent(theMCSum->GetMaximumBin()) + theMCSum->GetBinError(theMCSum->GetMaximumBin());
  if (theMaxMC > theMax) theMax = theMaxMC;
  
  //prepare the ratio band and plot
  TH1* theMCRatioBand = makeRatioBand(theMCSum);
  TH1* theRatioPlot = makeRatioPlot(hdata,theMCSum);
    
  TCanvas* can = new TCanvas();
  can -> SetLogy(isLog);
  
  TPad *pad1 = new TPad("pad1","top pad",0,0.30,1,1);
  pad1->SetBottomMargin(0.02);
  pad1->SetLeftMargin(0.13);
  pad1->Draw();
  TPad *pad2 = new TPad("pad2","bottom pad",0,0.0,1,0.30);
  pad2->SetTopMargin(0.02);
  pad2->SetLeftMargin(0.13);
  pad2->SetBottomMargin(0.4);
  pad2->SetGridy();
  pad2->Draw();
  
  pad1->cd();
  hs->Draw("hist");
  hdata->Draw("same,pe");
  if (drawSignal) hsignal->Draw("same,hist");
  if (drawLegend) leg->Draw("same");
  //hs->GetXaxis()->SetTitle(myAxisNameX);
  hs->GetYaxis()->SetTitle(myAxisNameY);
  hs->GetXaxis()->SetLabelSize(0.04);
  hs->GetYaxis()->SetLabelSize(0.04);
  hs->GetXaxis()->SetLabelOffset(0.025);

//______________________________________________________________________________
void CalibrateFit(Int_t nPar, Double_t time_limit, Double_t bad_frac)
{
    // Calibrate using iterative fitting.

    const Double_t convergence_factor = 0.05;

    // iterate calibrations
    for (Int_t i = 0; i >= 0; i++)
    {
        // create new histogram
        TH2* h = CreateHisto("fit_histo", nPar, gOff);

        // user info
        printf("Calibration iteration %d\n", i+1);

        // loop over elements
        Int_t n = 0;
        Int_t outside_n = 0;
        for (Int_t j = 0; j < nPar; j++)
        {
            // create projection
            TH1* hProj = (TH1*) h->ProjectionX(TString::Format("Proj_%d", j).Data(), j+1, j+1, "e");

            // check for filled histograms
            if (hProj->GetEntries())
            {
                // create fitting function
                TF1* func = new TF1(TString::Format("Func_%d", j).Data(), "gaus", -5, 5);
                func->SetParameter(0, 1);
                func->SetParameter(1, hProj->GetXaxis()->GetBinCenter(hProj->GetMaximumBin()));
                func->SetParameter(2, 0.1);

                // fit histogram
                hProj->GetXaxis()->SetRangeUser(-5, 5);
                hProj->Fit(func, "0Q");
                Double_t mean = func->GetParameter(1);

                // update offset
                gOff[j] = gOff[j] + convergence_factor * mean / gGain[j];
                n++;
                if (TMath::Abs(mean) > time_limit)
                    outside_n++;

                // clean-up
                delete func;
            }

            // clean-up
            delete hProj;
        }

        // clean-up
        delete h;

        // user info
        Double_t outside_frac = (Double_t)outside_n/(Double_t)n;
        printf("Element outside limits: %.1f%%\n", outside_frac*100.);
        if (outside_frac < bad_frac)
            break;
    }
}