C++ TH1F::SetMinimum方法代码示例

void produceHisto()
{ 
  TH1F * hmostimp =new TH1F("number most important",
			"Mean value number of over threshold digits produced by \
most important  particle",
                        20,0.,23.);
  gStyle->SetOptStat(0); 
  hmostimp->Fill(1.,33.51);
  hmostimp->Fill(5.,34.32);
  hmostimp->Fill(10.,35.78);
  hmostimp->Fill(15.,36.85);
  hmostimp->Fill(20.,37.61);
  hmostimp->Write();
  hmostimp->SetMinimum(32.);
  delete hmostimp;
  
  TH1F * hall =new TH1F("number all3",
			"Mean value Number of over threshold digits produced by \
alll three stored particles",
                        20,0.,23.);
  gStyle->SetOptStat(0); 
  hall->Fill(1.,77.05);
  hall->Fill(5.,76.905);
  hall->Fill(10.,71.9);
  hall->Fill(15.,69.7);
  hall->Fill(20.,65.15);
  hall->Write();
  hall->SetMinimum(32.);    
}

void EflowCASTOR()
{

  NCanvas(1,2,"data");

  TH1F * HDat;
  TH1F * HZ2;

  //  TFile *fzee = new TFile("ZDiffOutputfile_73Weight2010.root");
  TFile *fzee = new TFile("ZDiffOutputfile.root");
  // TFile *fzee = new TFile("ZDiffOutputfile.root");
  HDat  =  (TH1F*)fzee->Get("NVTX1_EnergyCastorModule_DATA");
  HZ2  =  (TH1F*)fzee->Get("NVTX1_EnergyCastorModule_Z2PY6");
  data->cd(1);
  NHSetMarker(HDat,2,20,0.8);
  HDat->SetMinimum(10);
  HDat->Draw();
  NStat(HDat,0);
  NSetTitle(HDat,"CASTOR Modules", "Entries");

  HDat->SetTitle("CASTOR Module Energy 2010");
  HZ2->Draw("SAME HIST");
  NHSetMarker(HDat,2,20,0.4);
  TLegend *legend = new TLegend(0.6,0.7,0.9,0.8);
  legend->SetTextFont(72);
  legend->SetTextSize(0.06);
  legend->SetBorderSize(0);


  legend->AddEntry(HDat,"Data2010","p");
  legend->AddEntry(HZ2,"Z2","l");
  legend->Draw();


  TH1F * Hra = (TH1F * ) HDat->Clone();
  //  Hra->GetXaxis->SetRange(0,12);
  NStat(Hra,0);
  NHSetMarker(Hra,4,20,0.8);
  Hra->SetMinimum(1);
  Hra->Divide(HZ2);
  NLogAxis(0,0);
  data->cd(2);
  Hra->Draw();
  Hra->SetTitle("Ratio DATA/MC");
  Hra->Fit("pol0");
  NSetTitle(Hra,"CASTOR Modules", "Ratio");
  NLogAxis(0,1);

}

TH1F* CreateEffHisto(TH1F* hGen, TH1F* hRec)
{
// create an efficiency histogram

  Int_t nBins = hGen->GetNbinsX();
  TH1F* hEff = (TH1F*) hGen->Clone("hEff");
  hEff->SetTitle("");
  hEff->SetStats(kFALSE);
  hEff->SetMinimum(0.);
  hEff->SetMaximum(110.);
  hEff->GetYaxis()->SetTitle("#epsilon [%]");

  for (Int_t iBin = 0; iBin <= nBins; iBin++) {
    Double_t nGen = hGen->GetBinContent(iBin);
    Double_t nRec = hRec->GetBinContent(iBin);
    if (nGen > 0) {
      Double_t eff = nRec/nGen;
      hEff->SetBinContent(iBin, 100. * eff);
      Double_t error = sqrt(eff*(1.-eff) / nGen);
      if (error == 0) error = 0.0001;
      hEff->SetBinError(iBin, 100. * error);
    } else {
      hEff->SetBinContent(iBin, -100.);
      hEff->SetBinError(iBin, 0);
    }
  }

  return hEff;
}

void DrawInsituCorr_vs_Eta(Str jetAlgo) {
  Nbins=sizeof(Ebins)/sizeof(double);
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);

  TH1F *h = new TH1F("","",100,-4.8,4.8);
  h->SetXTitle("Jet #eta_{det}");
  h->SetYTitle("Residual correction for data only");
  h->SetMaximum(1.1); h->SetMinimum(0.85);

  h->Draw();
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(-4.65,1,4.65,1);
  
  for (int i=0;i<Nbins;++i) {
    Graph *g = GetInsituGraph(Ebins[i]);
    FormatGraph(g,i); g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("p_{T} = %.0f GeV",Ebins[i]),x,y,i);
  }

  l->SetLineWidth(2);
  l->SetLineColor(kGray+2);
  l->SetLineStyle(1);

  tex->SetNDC(); tex->SetTextAlign(12);
  //tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,GetJetDesc(jetAlgo));
}

void DrawInsituCorr_vs_Pt(Str jetAlgo) {
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);
  TH1F *h = new TH1F("","",100,15,2500);
  h->SetXTitle("Jet p_{T}");
  h->SetYTitle("Residual correction for data only");
  h->SetMaximum(1.1); h->SetMinimum(0.85);
  h->Draw();
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(16,1,2400,1);
  for (int i=0;i<6;++i) {
    double eta=4.0*i/5;
    Graph *g = GetInsituGraphVsPt(eta);
    FormatGraph(g,i); g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("#eta = %.1f",eta),x,y,i);
  }

  l->SetLineWidth(2);
  l->SetLineColor(kGray+2);
  l->SetLineStyle(1);

  tex->SetNDC(); tex->SetTextAlign(12);
  //tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,GetJetDesc(jetAlgo));
}

void DrawJMS_vs_Eta(Str jetAlgo) {
  double massEbins[] = {50, 100, 250, 500, 750, 1500};
  double massNbins = sizeof(massEbins)/sizeof(double);
  myJES = new JetCalibrationTool(jetAlgo,_jesFile);
  TH1F *h = new TH1F("","",100,-4.8,4.8);
  h->SetXTitle("Jet #eta_{det}"); h->SetYTitle("Jet mass response");
  h->SetMaximum(1.4); h->SetMinimum(0.6);
  h->Draw();
  
  TLine *l = new TLine(); 
  l->SetLineWidth(2); l->SetLineColor(kGray+2); l->SetLineStyle(4);
  l->DrawLine(-4.65,1,4.65,1);

  for (int i=0;i<massNbins;++i) {
    Graph *g = GetJMSGraph(massEbins[i]);
    FormatGraph(g,i);
    g->Draw("PL");
    double x=0.52, y=0.3-0.05*i;
    if (i>2) { x=0.74; y=0.3-0.05*(i-3); }
    DrawLabel(Form("E = %.0f GeV",massEbins[i]),x,y,i);
  }
  
  tex->SetNDC(); tex->SetTextAlign(12);
  tex->DrawLatex(0.18,0.975,myJES->AbsoluteJES_Description());
  tex->DrawLatex(0.18,0.9,jetAlgo);
}

void readX(){

  TFile *_file0 = TFile::Open("withbdt.root");
  TTree* tree = gFile->Get("SelTree;1");

  TH1F* histo = new TH1F("m", "m", nbin, minMass, maxMass); histo->Sumw2();
  TH1F* bhisto = new TH1F("bm", "bm", nbin, minMass, maxMass); bhisto->Sumw2();
 

  float tzPsi; 
  tree->SetBranchAddress("tzPsi",&tzPsi);
  
  float m;
  tree->SetBranchAddress("mass",&m);
  float bdt;
  tree->SetBranchAddress("bdt",&bdt);

 
  int nevent = tree->GetEntries(); 
  int n = 0;

  for ( int i = 0; i < nevent; ++i){
    tree->GetEntry(i);
    histo->Fill(m);
    //cout << tzPsi << endl;
    if (bdt> -0.6&& tzPsi > 0.3) bhisto->Fill(m);
  }

  histo->SetMinimum(0); 
 histo->Draw();
 //bhisto->Divide(histo); 
 bhisto->Draw("HISTOSAME");
 }   

void seism() {

   TStopwatch sw; sw.Start();
   //set time offset
   TDatime dtime;
   gStyle->SetTimeOffset(dtime.Convert());

   TCanvas *c1 = new TCanvas("c1","Time on axis",10,10,1000,500);
   c1->SetFillColor(42);
   c1->SetFrameFillColor(33);
   c1->SetGrid();

   Float_t bintime = 1; //one bin = 1 second. change it to set the time scale
   TH1F *ht = new TH1F("ht","The ROOT seism",10,0,10*bintime);
   Float_t signal = 1000;
   ht->SetMaximum( signal);
   ht->SetMinimum(-signal);
   ht->SetStats(0);
   ht->SetLineColor(2);
   ht->GetXaxis()->SetTimeDisplay(1);
   ht->GetYaxis()->SetNdivisions(520);
   ht->Draw();

   for (Int_t i=1;i<2300;i++) {
      //======= Build a signal : noisy damped sine ======
      Float_t noise  = gRandom->Gaus(0,120);
      if (i > 700) noise += signal*sin((i-700.)*6.28/30)*exp((700.-i)/300.);
      ht->SetBinContent(i,noise);
      c1->Modified();
      c1->Update();
      gSystem->ProcessEvents(); //canvas can be edited during the loop
   }
   printf("Real Time = %8.3fs, Cpu Time = %8.3fs\n",sw.RealTime(),sw.CpuTime());
}

void PlotMinEtFromSim(Bool_t isPhos = kFALSE){
  gStyle->SetOptTitle(0);
  gStyle->SetOptStat(0);
  gStyle->SetOptFit(0);
  Float_t min = 0;
  float max = 1;
  TString filename, detname;
  if(isPhos){
    min = 0.655;
    max = 0.785;
    detname = "PHOS";
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.PHOS.LHC11a10a_bis.Run139465.root";
  }
  else{
    min = 0.58;
    max = 0.725;
    filename = "rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root";
    detname = "EMCal";
  }
  
  TFile *f = TFile::Open(filename, "READ");
  TList *l = dynamic_cast<TList*>(f->Get("out1"));
  TH1F *fHistSimulatedGammaEnergyAboveThreshold = l->FindObject("fHistSimulatedGammaEnergyAboveThreshold");
  TH1F *fHistSimulatedGammaEnergy = l->FindObject("fHistSimulatedGammaEnergy");
  SetStyles(fHistSimulatedGammaEnergyAboveThreshold,20,TColor::kRed);
  fHistSimulatedGammaEnergyAboveThreshold->Divide(fHistSimulatedGammaEnergy);

    TCanvas *c1 = new TCanvas("c1","Simulation",600,400);
    c1->SetTopMargin(0.02);
    c1->SetRightMargin(0.03);
    c1->SetLeftMargin(0.11745);
    c1->SetBottomMargin(0.11745);
    c1->SetBorderSize(0);
    c1->SetFillColor(0);
    c1->SetFillColor(0);
    c1->SetBorderMode(0);
    c1->SetFrameFillColor(0);
    c1->SetFrameBorderMode(0);
    fHistSimulatedGammaEnergyAboveThreshold->SetMaximum(max +0.1);
    fHistSimulatedGammaEnergyAboveThreshold->SetMinimum(min-0.1);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitle("Centrality bin");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitle("f_{minEt}");
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetLabelSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetYaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->GetXaxis()->SetTitleSize(0.06);
    fHistSimulatedGammaEnergyAboveThreshold->Draw();
    TLine *lineMin = new TLine(-0.5,min,19.5,min);
    lineMin->Draw();
    TLine *lineMax = new TLine(-0.5,max,19.5,max);
    lineMax->Draw();
    lineMin->SetLineColor(TColor::kBlue);
    lineMax->SetLineColor(TColor::kBlue);
    lineMin->SetLineStyle(2);
    lineMax->SetLineStyle(2);

    TString outfile = "/tmp/MinEtFromSim"+detname+".png";
    c1->SaveAs(outfile.Data());
}

void draw_ratio(std::vector<TH1F*> h,
	  TString name, TString xTitle,
	  double xmin, double xmax,
	  TString legHeader = "", bool legRIGHT = true, bool legTOP = true,
	  bool logX = false, bool stat = false, int rebin = -1, int orbin = -1,
	  TString option = "", int nclus = 99) {  //double ymin_ratio, double ymax_ratio,
  TCanvas* can = new TCanvas(name+"_ratio",name+"_ratio",900,450);
  can->cd();

  double legxmin = (legRIGHT ? 0.55 : 0.18);
  double legxmax = legxmin+0.25;
  double legymin = (legTOP ? 0.70 : 0.15);
  double legymax = legymin+0.15;
  TLegend* leg = new TLegend(legxmin,legymin,legxmax,legymax);
  if (legHeader!="") leg->SetHeader(legHeader);
  leg->SetTextSize(0.04);
  leg->SetFillColor(0);
  leg->SetLineColor(0);

  TString options = (option=="" ? "pe" : option);

  if (rebin>0) h[h.size()-1]->Rebin(rebin); //to rebin benchmark before divide
  for (size_t i=0; i<h.size(); i++) {
    //if(h[i]->GetNbinsX() != orbin) cout << "WARNING: orbin for " << h[i]->GetName() << " are " << h[i]->GetNbinsX() << endl; //debug - shift of h[][] wrt clu[][]
    if (rebin>0 && i<(h.size()-1)) h[i]->Rebin(rebin);
    TH1F* ratio = (TH1F*)h[i]->Clone("ratio_"+name);
    ratio->Sumw2();
    ratio->Divide(h[h.size()-1]); //benchmark is at the end.

    if (logX) gPad->SetLogx();
    ratio->SetMarkerStyle(20+i);
    ratio->SetMarkerSize(1.0); //1.2
    ratio->GetXaxis()->SetRangeUser(xmin,xmax);
    ratio->SetMinimum(-0.1);
    ratio->SetMaximum(4);
    if (i==0){    //just for the first one
     ratio->GetXaxis()->SetLabelSize(0.05);
     ratio->GetXaxis()->SetTitle(xTitle);
     ratio->GetXaxis()->SetTitleOffset(1);
     ratio->GetXaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitle("ratio");
     ratio->GetYaxis()->SetTitleSize(0.06);
     ratio->GetYaxis()->SetTitleOffset(0.7);
     ratio->GetYaxis()->SetLabelSize(0.05);
    }
    string nam = "";
    nam = translate(clu[nclus][0].c_str());
    if(i==(h.size()-1)) leg->AddEntry(h[i],nam.c_str(),"l"); //to print only benchmark (first in the list)
    if (i==1) options = options + (stat ? "sames" : "same"); //once is enought
    ratio->Draw(options);
  }
  leg->Draw("same");
  drawPrivate(0.04);
  can->Update();
  can->SaveAs(Outfolder+name+"_ratio.png");
}

void CheckEHF()
{

 
  TH1F * HDat;
  TH1F * HDatnoPUAft;
  TH1F * HDatnoPUBef;

  //  TFile *fzee = new TFile("ZDiffOutputfile_73Weight2010.root");
  TFile *fzee = new TFile("ZDiffOutputfile_Data2011_noPU.root");
  // TFile *fzee = new TFile("ZDiffOutputfile.root");
  HDat  =  (TH1F*)fzee->Get("NVTX1_minEHFZoom_DATA10");
  HDatnoPUAft  =  (TH1F*)fzee->Get("NVTX1noPUAft_minEHFZoom_DATA10");
  HDatnoPUBef  =  (TH1F*)fzee->Get("NVTX1noPUBef_minEHFZoom_DATA10");

  HDat->Scale(1./HDat->GetEntries());
  HDatnoPUAft->Add(HDatnoPUBef);
  HDatnoPUAft->Scale(1./HDatnoPUAft->GetEntries());
  HDatnoPUBef->Scale(1./HDatnoPUBef->GetEntries());
  HDatnoPUAft->SetMinimum(0.);
  NHSetMarker(HDatnoPUAft,2,20,0.8);
  NHSetLine(HDatnoPUAft,2,1,1.);
  HDatnoPUAft->Draw("E1");
  HDatnoPUAft->SetTitle("2011 Z->ee Data, minHF Energy");
  NSetTitle(HDatnoPUAft,"HF Energy [GeV]", "Fraction");


  NHSetMarker(HDatnoPUBef,6,20,0.8);
  NHSetLine(HDatnoPUBef,6,1,1.);
  // HDatnoPUBef->Draw("SAMEP");

  NHSetMarker(HDat,4,20,0.8);
  NHSetLine(HDat,1,1,1.);
  HDat->Draw("SAMEP");

  TLegend *legend = new TLegend(0.6,0.1,0.75,0.4);
  legend->SetTextFont(72);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);

  legend->AddEntry(HDatnoPUBef,"No PU before or after","p");
  //  legend->AddEntry(HDatnoPUBef,"No PU before","p");
  // legend->AddEntry(HDatnoPUAft,"No PU after","p");
  legend->AddEntry(HDat,"PU before && after","p");
  legend->Draw();






  //  NHSetLine(HDatnoPUBef,3,2,1.);
  // HDatnoPUBef->Draw("SAME HIST");

}

//------------------------------------//
// Main
//------------------------------------//
void AnaEField()
{

  // Make Canvas 
  TCanvas* can = makeCanvas("can");
  can->SetLogy();
  can->SetGridx();
  can->SetGridy();

  // Define the function here with parameters:
  // [0] = L
  // [1] = freq
  // [2] = k
  // [3] = n
  // [4] = sqrt(2pi) * mu * mu0
  TString func = "[4]*[0]*[1]*sin(x/57.2957)*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";
  //TString func = "[4]*[1]*sin(x/57.2957)*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";
  //TString func = "[4]*[0]*[1]*exp(-pow([2]*[0],2)*pow(cos(x/57.2957)-1/[3],2)/2)";

  // Now fix some of the variables:
  double freq = 1e9; //600e6; // Hz
  double c    = 3e8;   // m/s
  double n    = 1.3; //1.78;
  double k    = 2*TMath::Pi() * n * freq/c; // 1/m
  double Const= sqrt(2*TMath::Pi()) * 4 * TMath::Pi() * 1e-7 * 1;

  // Make a dummy histogram
  TH1F* h = makeHist("h",100,0,90,"Angle [deg]", "Field Strength",kBlack,0);
  h->Fill(1e10);
  h->SetMinimum(10);
  h->SetMaximum(3000);
  h->Draw();

  // Make two instances of the function
  TF1* f0 = new TF1("f0",func.Data(),0,90);
  setParam(f0,1.2,freq,k,n,Const,kBlack);

  TF1* f1 = new TF1("f1",func.Data(),0,90);
  setParam(f1,0.1,freq,k,n,Const,kBlue);
 
  // Draw
  f0->Draw("same");
  f1->Draw("same");
  
  // Add Legend
  TLegend* leg = makeLegend(0.15,0.3,0.7,0.9);
  leg->SetHeader("f = 600 MHz");
  leg->AddEntry(f0,"L = 1.2m","l");
  leg->AddEntry(f1,"L = 0.1m","l");
  leg->Draw("same");
  
  //can->SaveAs("../plots/quickAnalytic_600MHz.png");

}

void FillGoodRun(int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
     for(int icent=0;icent<ncent;icent++){
 //     for(int ihar=0;ihar<nhar;ihar++){
 //         if(ihar!=1) continue;
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      //std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
        //fin = TFile::Open(Form("/phenix/plhf/xuq/taxi/%s%s/%d/data/%s.root",dataset.Data(),pro.Data(),taxi,GetRun(irun).Data()));
         fin = TFile::Open(Form("/gpfs/mnt/gpfs02/phenix/plhf/plhf1/xuq/phenix/flow/Run16dAu/work/62GeV/output/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_0_0_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
   //   }
     }
}

CheckEnabledChannels(const Char_t *runlist)
{

  ifstream is(runlist);
  Char_t buf[4096];
  Int_t run[1024];
  Int_t nrun = 0;
  while(!is.eof()) {
    is.getline(buf, 4096);
    if (is.eof()) break;
    run[nrun] = atoi(buf);
    printf("added run number %d\n", run[nrun]);
    nrun++;
  }
  printf("%d runs added\n", nrun);
  is.close();

  TH1F *hActive = new TH1F("hActive", "active channels;run;fraction", nrun, 0, nrun);
  TH1F *hReadout = new TH1F("hReadout", "good readout;run;fraction", nrun, 0, nrun);
  for (Int_t irun = 0; irun < nrun; irun++) {
    hr = CheckEnabledChannels(run[irun], kTRUE);
    ha = CheckEnabledChannels(run[irun], kFALSE);
    hReadout->SetBinContent(irun + 1, hr->Integral());
    hActive->SetBinContent(irun + 1, ha->Integral());
    hReadout->GetXaxis()->SetBinLabel(irun + 1, Form("%d", run[irun]));
    delete hr; delete ha;
  }
  
  hReadout->SetMarkerStyle(20);
  hReadout->SetMarkerColor(4);
  hActive->SetMarkerStyle(25);
  hActive->SetMarkerColor(2);
  hReadout->Sumw2();
  hActive->Sumw2();
  hReadout->Divide(hReadout, hActive, 1., 1., "B");
  hActive->Scale(1. / 152928.);
  hReadout->SetMinimum(0.);
  hReadout->SetMaximum(1.);
  hReadout->Draw("E");
  hActive->Draw("E, same");
  TLegend *l = gPad->BuildLegend();
  l->SetFillStyle(0);

}

void FillGoodRun(int icent,int ihar){
    float pi = acos(-1.0);
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    if(nrun<0) exit(1);
       for(int isub=0;isub<nsub;isub++){
            str = choosesub(isub);
            if(str=="ABORT") continue;
            for(int irun=0;irun<nrun;irun++){
      std::cout<<icent<<" "<<ihar<<" "<<isub<<" "<<irun<<std::endl;
         fin = TFile::Open(Form("/store/user/qixu/flow/Run16dAu/62GeV/%s",GetRun(irun).Data()));
        TH1F* hpsi = new TH1F("psi","psi",100,-pi,pi);
        for(int ibbcz=0;ibbcz<nbbcz;ibbcz++){
          TH1F* hpsitemp = (TH1F*)fin->Get(Form("psiFla_%d_%d_%d_%d",icent,ibbcz,ihar,isub));
          hpsi->Add(hpsitemp);
        }
        TF1 *fun = new TF1("fun","pol0",-pi,pi);
      if(hpsi->GetEntries()>1000){
	//hpsi->SetMarkerStyle(20);
	//hpsi->SetMarkerSize(0.6);
	//hpsi->SetMarkerColor(4);
	hpsi->SetMinimum(10);
	hpsi->Fit("fun","QR0");
	//float par=fun->GetParameter(0);
	//hpsi->SetMaximum(1.5*par);
	//hpsi->Draw();
	GoodRunFit[icent][ihar][isub][irun] = fun->GetChisquare()/fun->GetNDF();
        fin->Close();
      }
      else{
        GoodRunFit[icent][ihar][isub][irun] = -9999;
        fin->Close();
     }
     // GoodRunFit[icent][ihar][isub][irun] = 1.;
    }
       }
}