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

//.........这里部分代码省略.........
    hAllFound->Draw("same");  // MC pt distribution
  }
 
 
  /*

  .L ~/ITSupgrade/BuildDetector/DetectorK.cxx+
  
  // All NEW
  DetectorK its("ALICE","ITS");
  its.MakeAliceAllNew(0);
  its.SetMaxRadiusOfSlowDetectors(0.01);
  its.SolveViaBilloir(0);
  TGraph *c = its.GetGraphRecoEfficiency(0,3,2);
  c->Draw("C");


  // Current
  DetectorK its("ALICE","ITS");
  its.MakeAliceCurrent(0,0);
  its.SetMaxRadiusOfSlowDetectors(0.01);
  its.SolveViaBilloir(0);
  TGraph *c = its.GetGraphRecoEfficiency(0,4,2);
  c->Draw("C");

  */

  TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900);
  c3->Divide(2,1); c3->cd(1);
  // Impact parameter

  // Impact parameter resolution ---------------
  h2Ddca->Draw("colz");
  h2Ddca->FitSlicesY() ;
  TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same");
  TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw();
  TGraphErrors *d0 = new TGraphErrors(); 
  for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) {
    d0->SetPoint(     ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns
    d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns
  }
  d0->SetMarkerStyle(21);
  d0->SetMaximum(200);  d0->SetMinimum(0);
  d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
  d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)");
  d0->SetName("dca");  d0->SetTitle("DCAvsPt");

  c3->cd(1);  h2Ddca->Draw("surf2");
  c3->cd(2);  d0->Draw("APE");

  // PT RESOLUTION ------------
  TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900);
  c4->Divide(2,1); c4->cd(1);
  // Impact parameter
  h2Dpt->Draw("colz");
  h2Dpt->FitSlicesY() ;
  TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same");
  TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw("");
  TGraphErrors *gPt = new TGraphErrors(); 
  for (Int_t ibin =1; ibin<=dPtRMS->GetXaxis()->GetNbins(); ibin++) {
    gPt->SetPoint(     ibin-1,dPtRMS->GetBinCenter(ibin),dPtRMS->GetBinContent(ibin)); 
    gPt->SetPointError(ibin-1,0,dPtRMS->GetBinError(ibin)); 
  }
  gPt->SetMarkerStyle(21);
  gPt->SetMaximum(20);  gPt->SetMinimum(0);
  gPt->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
  gPt->GetYaxis()->SetTitle("relative momentum resolution (%)");
  gPt->SetName("dPt");  gPt->SetTitle("DPTvsPt");

  c4->cd(1);  h2Dpt->Draw("surf2");
  c4->cd(2);  gPt->Draw("APE");


  // EXPORT --------

  TFile f("histos.root","RECREATE");

  hMultCount->Write();
  hAllMC->Write();
  hAllFound->Write();
  hImperfect->Write();
  hPerfect->Write();
  hNoMCTrack->Write();

  hPurity->Write();
  hEff->Write();
  hFake->Write();
  hAnna->Write();

  h2Ddca->Write();
  d0->Write();

  h2Dpt->Write();
  gPt->Write();

  f.Close();

  return;

}

void make1DLimitHH(TString combine_dir,bool blind=true){
  //TString combine_dir = "test_runSusyHgg/signalInj_sms_ChiHH_0_175/";

  TGraph obser( (275-125)/25 );
  TGraph graph( (275-125)/25 );
  TGraphAsymmErrors error( (275-125)/25 );

    for(int m=125;m<501;m+=25) {
      int i=(m-125)/25;
      TFile limit_file(Form("%s/higgsCombineChiHH_0_%d.Asymptotic.mH120.root",combine_dir.Data(),m) );
      TTree *limit_tree = (TTree*)limit_file.Get("limit");
      TTreeFormula limit_form("get_limit","limit",limit_tree);

      limit_tree->GetEntry(1);
      float down = limit_form.EvalInstance();
      limit_tree->GetEntry(2);
      float exp = limit_form.EvalInstance();
      limit_tree->GetEntry(3);
      float up = limit_form.EvalInstance();
      limit_tree->GetEntry(5);
      float obs = limit_form.EvalInstance();

      graph.SetPoint(i,float(m),exp);
      error.SetPoint(i,float(m),exp);
      error.SetPointError(i,0,0,exp-down,up-exp);

      obser.SetPoint(i,float(m),obs);
    }

    TGraphErrors* theo = getTheoXSec("/home/amott/HggApp/SusyHgg/xsecs/Higgsino.txt");

    TCanvas cv;
    cv.SetLogy();

    theo->SetMaximum(1e2);
    theo->SetMinimum(1e-2);
    theo->GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    theo->GetXaxis()->SetTitle("m_{chargino}");
    theo->SetFillColor(kBlue);
    theo->SetLineStyle(kDotted);
    theo->SetLineWidth(2.0);
    
    error.SetMaximum(1e2);
    error.SetMinimum(1e-2);
    error.GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    error.GetXaxis()->SetTitle("m_{chargino}");
    error.SetFillColor(kGreen);
    error.Draw("A3");

    theo->Draw("3C");

    graph.SetLineStyle(kDashed);
    graph.SetLineWidth(2);
    graph.Draw("C");

    obser.SetLineStyle(1);
    obser.SetLineWidth(2);
    if(!blind) obser.Draw("C");

    TLegend leg(0.7,0.7,0.85,0.85);
    leg.SetFillColor(0);
    leg.SetBorderSize(0);
    leg.AddEntry(&graph,"expected","l");
    leg.AddEntry(&error,"expected #pm1#sigma","F");
    leg.AddEntry(theo,"theoretical","f");
    if(!blind)     leg.AddEntry(&obser,"observed","l");

    leg.Draw("SAME");

    TLatex prelim(0.65,0.96,"CMS Preliminary");
    prelim.SetNDC();
    prelim.SetTextSize(0.045);
    prelim.Draw();

    TLatex lbl(0.5,0.86,"#sqrt{s} = 8 TeV  #int L dt = 19.78 fb^{-1}");
    lbl.SetNDC();
    lbl.SetTextSize(0.045);
    lbl.Draw();

    cv.SaveAs(combine_dir+"expected_exclusion_HH_1D.png");


}

//.........这里部分代码省略.........
    sampling->SetPoint(ic,ICval[ic]/100.,fit->GetParameter(0));
    sampling->SetPointError(ic,0,fit->GetParError(0));
    cval = sqrt(pow(fit7->GetParameter(1),2)-pow(y0_7,2));
    constantSR7->SetPoint(ic,ICval[ic]/100.,cval);
    constantSR7->SetPointError(ic,0,fit7->GetParameter(1)*fit7->GetParError(1)/cval);    
    //constantSR7->SetPoint(ic,ICval[ic]/100.,fit7->GetParameter(1));
    //constantSR7->SetPointError(ic,0,fit7->GetParError(1));
    samplingSR7->SetPoint(ic,ICval[ic]/100.,fit7->GetParameter(0));
    samplingSR7->SetPointError(ic,0,fit7->GetParError(0));

    if (ic==0) {
      constant->GetPoint(0,x0,y0);
      constantSR7->GetPoint(0,x0_7,y0_7);
      cval = sqrt(pow(fit->GetParameter(1),2)-pow(y0,2));
      constant->SetPoint(ic,ICval[ic]/100.,cval);
      constant->SetPointError(ic,0,fit->GetParameter(1)*fit->GetParError(1)/cval);
      cval = sqrt(pow(fit7->GetParameter(1),2)-pow(y0_7,2));
      constantSR7->SetPoint(ic,ICval[ic]/100.,cval);
      constantSR7->SetPointError(ic,0,fit7->GetParameter(1)*fit7->GetParError(1)/cval);
    }

  }

  mycReso->Update();
  mycReso->Print("PLOTS/ResolutionFitvsIC.pdf");
  
  TLegend *leg = new TLegend(0.6,0.6,0.8,0.8);
  leg->SetFillColor(10);
  leg->AddEntry(sampling,"3#times3 cm^{2}","P");
  leg->AddEntry(samplingSR7,"All detector","P");
  mycR->cd();
  gPad->SetGridy(1);
  sampling->GetYaxis()->SetTitle("Sampling term (GeV^{1/2})");
  sampling->SetMinimum(0.2);
  sampling->SetMaximum(0.3);
  sampling->Draw("APE");
  samplingSR7->Draw("PEsame");
  lat.SetTextSize(0.04);
  sprintf(buf,"Single #gamma, #eta=2.1");
  lat.DrawLatexNDC(0.2,0.87,buf);
  lat.DrawLatexNDC(0.01,0.01,"HGCAL G4 standalone");

  leg->Draw("same");
  mycR->Update();

  mycR->Print("PLOTS/SamplingvsIC.pdf");

  mycC->cd();
  gPad->SetLogx(1);
  gPad->SetGridy(1);
  gStyle->SetOptFit(0);
  //gStyle->SetStatH(0.1);
  //gStyle->SetStatW(0.2);

  constant->GetYaxis()->SetTitle("Constant from intercalib.");
  constant->SetMinimum(0);
  constant->SetMaximum(0.08);
  constant->Draw("APE");
  //constantSR7->Draw("PEsame");
  sprintf(buf,"Single #gamma, #eta=2.1");
  lat.DrawLatexNDC(0.2,0.87,buf);
  lat.DrawLatexNDC(0.01,0.01,"HGCAL G4 standalone");

  TF1 *BE = new TF1("BE","sqrt([0]*[0] + pow(x*1/sqrt([1]),2))",0,1);
  BE->SetParameters(0,30);
  //  BE->SetParLimits(0,1,1);

void plotMerged(Bool_t onlyPlot=0) {

  gStyle->SetPalette(1);
 
  TFile f("histoSum.root","UPDATE");

  TH1F* hAllMC = f.Get("allMC");
  TH1F* hAllFound= f.Get("allFound");
  TH1F* hImperfect= f.Get("imperfect");
  TH1F* hPerfect= f.Get("perfect");
  TH1F* hNoMCTrack= f.Get("noMCtrack");
  
  
  // have to be recalculated
  TH1F* hPurity = f.Get("purity");
  TH1F* hEff= f.Get("efficiency");
  TH1F* hFake= f.Get("fake");
  TH1F* hAnna= f.Get("annaEff");

  TH2D* h2Ddca= f.Get("dca2D");
  TGraphErrors *d0= f.Get("dca");

  TH2D* h2Dpt= f.Get("dPt2D");
  TGraphErrors *gPt= f.Get("dPt");


  if (!onlyPlot) {
    /*    // Get Errors right
    hAllMC->Sumw2();   
    hAllFound->Sumw2();
    hPerfect->Sumw2(); 
    hImperfect->Sumw2(); 
    h2Dpt->Sumw2();
    h2Ddca->Sumw2();
    */

    // Efficiencies - and normalize to 100%
    
    TF1 f1("f1","100+x*0",0.,1.e3);
    
    hPurity->Divide(hPerfect,hAllFound,1,1,"b"); 
    hPurity->Multiply(&f1);
    hPurity->SetMarkerColor(kGreen);
    hPurity->SetMarkerStyle(21);
    hPurity->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hPurity->SetStats(0);
    
    hPurity->GetYaxis()->SetRangeUser(0,100);
    hPurity->SetTitle("Efficiency & Purity");
    
    hEff->Divide(hPerfect,hAllMC,1,1,"b");
    hEff->Multiply(&f1);
    hEff->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hEff->SetMarkerColor(kBlue);
    hEff->SetMarkerStyle(21);
    hEff->SetStats(0);
    
    hFake->Divide(hImperfect,hAllMC,1,1,"b");
    hFake->Multiply(&f1);
    hFake->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hFake->SetMarkerColor(kRed);
    hFake->SetMarkerStyle(21);
    hFake->SetStats(0);
    
    hAnna->Divide(hAllFound,hAllMC,1,1,"b");
    hAnna->Multiply(&f1);
    hAnna->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    hAnna->SetMarkerColor(kBlack);
    hAnna->SetMarkerStyle(21);
    hAnna->SetStats(0);
    
    
    // Impact parameter resolution ---------------
    TCanvas *c3 = new TCanvas("c3","impact");//,200,10,900,900);
    c3->Divide(2,1); c3->cd(1);
    h2Ddca->DrawCopy("colz");
    h2Ddca->FitSlicesY() ;
    TH2D *dcaM = (TH2D*)gDirectory->Get("dca2D_1"); dcaM->Draw("same");
    TH2D *dcaRMS = (TH2D*)gDirectory->Get("dca2D_2"); //dcaRMS->Draw();
    TGraphErrors *d0 = new TGraphErrors(); 
    for (Int_t ibin =1; ibin<=dcaRMS->GetXaxis()->GetNbins(); ibin++) {
      d0->SetPoint(     ibin-1,dcaRMS->GetBinCenter(ibin),dcaRMS->GetBinContent(ibin)*1e4); // microns
      d0->SetPointError(ibin-1,0,dcaRMS->GetBinError(ibin)*1e4); // microns
    }
    d0->SetMarkerStyle(21);
    d0->SetMaximum(200);  d0->SetMinimum(0);
    d0->GetXaxis()->SetTitle("transverse momentum p_{t} (GeV)");
    d0->GetYaxis()->SetTitle("R-#phi Pointing Resolution (#mum)");
    d0->SetName("dca");  d0->SetTitle("DCAvsPt");
    //  c3->cd(1);  h2Ddca->Draw("surf2");
    c3->cd(2);  d0->Draw("APE");
    
    // PT RESOLUTION ------------
    TCanvas *c4 = new TCanvas("c4","pt resolution");//,200,10,900,900);  
    c4->Divide(2,1); c4->cd(1);
    h2Dpt->DrawCopy("colz");
    h2Dpt->FitSlicesY() ;
    TH2D *dPtM = (TH2D*)gDirectory->Get("dPt2D_1"); dPtM->Draw("same");
    TH2D *dPtRMS = (TH2D*)gDirectory->Get("dPt2D_2"); // dPtRMS->Draw("");
    TGraphErrors *gPt = new TGraphErrors(); 
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	  hSmapP->SetBinContent(ii+1,jj/2+1,rSP);
	}
	else {
	  printf("%2d : %2d / %2d / 2 %7.3f %8.4f %8.4f\n",nIter,ieta,iphi,dcorrS,corrS,corrS*drSP);
	  hSmapP->SetBinContent(ii+1,jj/2+1,rSP);
	  hSmapPc->SetBinError(ii+1,jj/2+1,corrS*drSP);
	  break;
	}
	if (nIter==nIterN-1) {
	  printf("%2d : %2d / %2d / 2 %7.3f %8.4f %8.4f\n",nIter,ieta,iphi,dcorrS,corrS,corrS*drSP);
	}
      }
    }
  }
  //fila->Close();

  cout<<endl<<"Rings :  "<<endl;
  cout<<"       E L        "<<"E S        "<<"eta     "<<"delta eta"<<endl;
  double xeta[13], weta[13], reta[13];
  for (int i=29;i<=41;i++) {
    xeta[i-29]=(etaBounds[i-28]+etaBounds[i-29])/2;
    weta[i-29]=(etaBounds[i-28]-etaBounds[i-29]);
    mLEphi[i-29]=mLEphi[i-29]*36/weta[i-29];
    mSEphi[i-29]=mSEphi[i-29]*36/weta[i-29];
    dmLEphi[i-29]=dmLEphi[i-29]*36/weta[i-29];
    dmSEphi[i-29]=dmSEphi[i-29]*36/weta[i-29];
    if (i>39) {  mLEphi[i-29]/=2; mSEphi[i-29]/=2; dmLEphi[i-29]/=2; dmSEphi[i-29]/=2; }
    reta[i-29] = mSEphi[i-29]/mLEphi[i-29];
    cout<<i<<" :  "<<mLEphi[i-29]<<"    "<<mSEphi[i-29]<<"    "<<xeta[i-29]<<"   "<<weta[i-29]<<endl;
  }
  TCanvas *cgL = new TCanvas("cgL","cgL",300,300,600,600);
  TGraphErrors *grL = new TGraphErrors(13,xeta,mLEphi,0,dmLEphi);
  grL->SetTitle("HFP L;#eta;E_{Ring} / #Delta#eta_{Ring} ,  GeV");
  grL->SetMinimum(0);
  grL->SetMarkerStyle(20);
  grL->Draw("1+PAl");
  cgL->Print("pictHFplot/etaProfHFPL.gif");
  mSEphi[12]/=2; mSEphi[11]/=2;
  TCanvas *cgS = new TCanvas("cgS","cgS",300,300,600,600);
  TGraphErrors *grS = new TGraphErrors(13,xeta,mSEphi,0,dmSEphi);
  grS->SetTitle("HFP S;#eta;E_{Ring} / #Delta#eta_{Ring} ,  GeV");
  grS->SetMinimum(0);
  grS->SetMarkerStyle(20);
  grS->Draw("1+PAl");
  cgS->Print("pictHFplot/etaProfHFPS.gif");
  TCanvas *crg = new TCanvas("crg","crg",300,300,600,600);
  TGraphErrors *rg = new TGraphErrors(13,xeta,reta,0,0);
  rg->SetTitle("HFP;#eta;E(S) / E(L)");
  rg->SetMinimum(0);
  rg->Draw("1+PAl");
  crg->Print("pictHFplot/SoverLetaHFP.gif");

  TCanvas *cL0 = new TCanvas("cL0","cL0",0,0,650,600);
  hLmapP0->Draw("colz");
  cL0->Update();
  TCanvas *cS = new TCanvas("cS0","cS0",1000,0,650,600);
  hSmapP0->Draw("colz");
  cS0->Update();

  //TFile *histf = new TFile("HFPmc.root","RECREATE");

  FILE *ft1;
  //sprintf(ctit,"corrHFPmc_%d_%d.txt",((int) Ethr1),((int) Ethr2));
  sprintf(ctit,"corrHFP_%s_%d_%d.txt",ftit,((int) Ethr1),((int) Ethr2));
  if ((ft1 = fopen(ctit,"w"))==NULL){               // Open new file
    printf("\nNo file %s open => EXIT\n\n",file);

//.........这里部分代码省略.........
  double wpDef;
  double clsWpDef;
  ifstream wpChannel;
  
  rolls.open("data/Barrel.txt");
  rolls>>rollName;      
  while (!rolls.eof()){
    //Being charged the working point per channel associated with the current roll
    /*    wpChannel.open(("results/"+rollName+"/wpChannel.txt").c_str());
    wpChannel>>wpDef;
    cout<<"wpDef= "<<wpDef<<endl;
    wpChannel.close();
    */
    //Being charged the run data for the current roll
    runsData.open(("results/"+rollName+"/runsData.txt").c_str());
    for (int n=0;n<runsAmmount;n++){
      runsData>>HV[n]>>EFF[n]>>ERR[n]>>EXP[n]>>CLS[n];
      //The first three runs (lower high voltage) are not plotted
      if (n>2){
	hv[n-3]=HV[n];
	cls[n-3]=CLS[n];
      }
    }
    runsData.close();

    ///////////////////////////////// The error in this cases is defined in this way,  this is something that must be improved    
    for(int n=0;n<runsAmmount;n++){
      if (n>2){
	exc[n-3]=0.0001;
	ecls[n-3] = 0.0001;
      }
    }
    
    //////////////////////////////////////////////////////
    //Being charged the fit data obtained for the current roll
    hvcls = new TGraphErrors(runsHavingInAccount, hv, cls, exc, ecls);
    cout<<rollName;
    fitData.open(("results/"+rollName+"/fitDataCLS.txt").c_str());
    fitData>>ca>>cb>>chi2>>wp>>clswp>>wpDef>>clsWpDef;
    cout<<ca<<" "<<cb<<" "<<wp<<" "<<clswp<<endl;
    fitData.close();
    f1 = new TF1("f1",expFunc,8.5,9.9,2); //3 es el numero de parametros del fit//
    f1->SetParameter(0, ca);
    f1->SetParameter(1, cb);//
    
    //Being made the plot of Working Point and Knee  
    knee = wp - 0.100;
    clsknee=expFuncAMano(knee,ca,cb);
    hvc[0]=wp - 0.100;
    hvc[1]=wp;
    CLSc[0]=clsknee;
    CLSc[1]=clswp;
    TGraph *gr2 = new TGraph(2,hvc,CLSc);
    
    //Being made the plot of Working Point per Channel
    wpCh[0]=wpDef;
    clsWpCh[0]=clsWpDef;
    TGraph *gr4 = new TGraph(1,wpCh,clsWpCh);
    
    //Being made the plot of the sigmoid
    for (int k=0; k<P; k++){
      x[k]=xmin+k*(xmax-xmin)/(P-1);
      y[k]=expFuncAMano(x[k],ca,cb);
      x[k]=x[k];
    }
    TGraph *gr3 = new TGraph(200,x,y);
   
    
    //Being set plot parameters 
    hvcls->SetLineColor(kRed);
    hvcls->SetMarkerStyle(20);
    hvcls->SetMarkerSize(2.0);
    hvcls->SetMinimum(-0.01);
    hvcls->SetMaximum(6);
    TAxis *axis = hvcls->GetXaxis();
    axis->SetLimits(8.5,9.9); 
    hvcls->SetTitle(("CLS vs HV_Eff " + rollName).c_str());
    hvcls->GetXaxis()->SetTitle("HV_Eff(kV)");
    hvcls->GetYaxis()->SetTitle("CLS");
    hvcls->Draw("AP");   
    gr3->SetLineColor(kBlue);
    gr3->Draw("C");// superimpose the second graph by leaving out the axis option "A"
    gr2->SetMarkerStyle(28);
    gr2->SetMarkerSize(3);
    gr2->SetLineColor(kBlue);
    gr2->Draw("P");
    gr4->SetMarkerStyle(24);
    gr4->SetMarkerSize(5);
    gr4->SetLineColor(kRed);
    gr4->Draw("P");
    
    //Being stored the plot as png file
    gSystem->mkdir(("results/"+rollName).c_str());
    c1->SaveAs(("results/"+rollName+"/CLSvsHV.png").c_str());
    c1->Clear();     
    //}
    rolls>>rollName;
  }
  exit(0);
}

//.........这里部分代码省略.........
  //l1->Scale(PDRate/nentries);

  for(unsigned int ii=0; ii<l1Names->size(); ++ii)
    l1->GetXaxis()->SetBinLabel(ii+1,l1Names->at(ii).c_str());
  //l1->GetYaxis()->SetTitle("L1 Rate @4E33 [Hz]");
  l1->SetMaximum(l1->GetMaximum()+200);
  l2->SetLineColor(kRed);
  
  l1->Draw();
  l2->Draw("same");
  c3->Update();

  // TGaxis *l1axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),gPad->GetUxmax(), gPad->GetUymax(),
  // 			      l1->GetMinimum()*lumiScaleFactor,
  // 			      l1->GetMaximum()*lumiScaleFactor,510,"+L");

  // c2->SetTicky(0);
  // l1axis->SetLineColor(kRed);
  // l1axis->SetLabelColor(kRed);
  // l1axis->SetTextColor(kRed);
  // l1axis->SetTitleOffset(1.3);
  // l1axis->SetLabelSize(0.03);
  // l1axis->SetTitle("L1 Rate @1E34 [Hz]");
  // l1axis->Draw();
  
  
  //return 0;

  //##############################################
  //##############################################

  //book graphs and plots
  TGraphErrors* totRateVsCut = new TGraphErrors();
  totRateVsCut->SetMinimum(0);
  TGraphErrors* pureRateVsCut450 = new TGraphErrors();
  TGraphErrors* pureRateVsCut280 = new TGraphErrors();

  //loops
  int bin = 0;
  for (int cut = 350; cut < 500; cut=cut+10)
    {
      int mjjPassed = 0;
      int HT250Calo_Passed = 0;
      int excl410_passed = 0;
      int excl250_passed = 0;
      for (Long64_t jentry=0; jentry<nentries;++jentry) 
	{
	  tt->GetEntry(jentry);

	  if (hltAccept->at(HT250Calo) == 1)
	    ++HT250Calo_Passed;

	  //if (caloMjj > cut && !hltAccept->at(HT410PF))
	  if (caloMjj > cut && l1Accept->at(L1scenario) == 1 && hltAccept->at(HT410PF)==0)
	    ++excl410_passed;
	  if (caloMjj > cut && l1Accept->at(L1scenario)==1 && hltAccept->at(HT250Calo)==0)
	    ++excl250_passed;
	  if (caloMjj > cut && l1Accept->at(L1scenario)==1)
	    ++mjjPassed;

	  // if (hltAccept->at(HT250Calo) == 0 && mjj > cut)
	  //   std::cout << "ref trigger doesn't completely cover cut at " << cut << std::endl;
	}
      // float mjjTotalRate = (float)mjjPassed/(float)HT250Calo_Passed*HT250Calo_rate;
      // float mjjPureRate = (float)exclPassed/(float)HT250Calo_Passed*HT250Calo_rate;

TGraphErrors *  drawAll()
{
  TFile *_file1 = TFile::Open("here.root");

  // use also an old style file??

  bool old = false;
  TFile *_file2;
  if (old == true){
  _file2= TFile::Open("jetTagAnalysis.root");
  }
  vector<TFile *> files;
  vector<string> algos;
  algos.push_back("trackCountingHighPurBJetTags");
  algos.push_back("trackCountingHighEffBJetTags");
  algos.push_back("jetProbabilityBJetTags");
  algos.push_back("jetBProbabilityBJetTags");
  algos.push_back("simpleSecondaryVertexBJetTags");
  algos.push_back("combinedSecondaryVertexBJetTags");
  algos.push_back("softMuonBJetTags");
  algos.push_back("softMuonByIP3dBJetTags");
  algos.push_back("softMuonByPtBJetTags");
  algos.push_back("softElectronBJetTags");

  TLegend * leg = new TLegend(0.2,0.65,0.4,0.9);
  TCanvas * c1 = new TCanvas();
  c1->SetLogy();  
  c1->SetGridy();  
  c1->SetGridx();  
  for(int i = 0 ; i < algos.size() ; i++)
   {
      cout << algos[i] << endl;
      // get the eff vs dicriminant

      TH1F* effdiscrb =  (TH1F *) getHistogram(_file1,algos[i],"effVsDiscrCut_discr","GLOBAL", "B");
      TH1F* effdiscruds = (TH1F *) getHistogram(_file1,algos[i],"effVsDiscrCut_discr","GLOBAL", "DUS");
      
      cout <<" HISTOS "<<effdiscrb <<" " <<effdiscruds<<endl;

      TGraphErrors * h = computeGraph(effdiscrb,effdiscruds);
      //      h->SetMaximum(1.);
            h->SetMinimum(1e-5);
      //      h->GetXaxis()->SetLimits(0,.2);
      //      h->GetYaxis()->SetLimits(1e-4,1);
     cout << h << endl;
     //     return h;
     setStyle(i,h);
     if(i==0) h->Draw("ALP"); else h->Draw("lpsame"); 
     leg->AddEntry(h,algos[i].c_str(),"p");
   }

  if (old == true){
    
    for(int i = 0 ; i < algos.size() ; i++)
      {
	cout << algos[i] << endl;
	// get the eff vs dicriminant
	
	TH1F* effdiscrb =  (TH1F *) getHistogram2(_file2,algos[i],"effVsDiscrCut_discr","GLOBAL", "B");
	TH1F* effdiscruds = (TH1F *) getHistogram2(_file2,algos[i],"effVsDiscrCut_discr","GLOBAL", "DUS");
      
	cout <<" HISTOS "<<effdiscrb <<" " <<effdiscruds<<endl;
	
	TGraphErrors * h = computeGraph(effdiscrb,effdiscruds);
	//      h->SetMaximum(1.);
	h->SetMinimum(1e-5);
	//      h->GetXaxis()->SetLimits(0,.2);
	//      h->GetYaxis()->SetLimits(1e-4,1);
	cout << h << endl;
	//     return h;
     setStyle(i+10,h);
     //     if(i==0) h->Draw("ALP"); else 
     h->Draw("lpsame"); 
     leg->AddEntry(h,algos[i].c_str(),"p");
      }
    
  }



  leg->Draw("same");

}

//.........这里部分代码省略.........
  tt->SetBranchAddress("l1DeltaEta", &l1DeltaEta_, &b_l1DeltaEta);

  tt->SetBranchAddress("l1JetPt", &l1JetPt_, &b_l1JetPt);
  tt->SetBranchAddress("l1JetEta", &l1JetEta_, &b_l1JetEta);
  tt->SetBranchAddress("l1JetPhi", &l1JetPhi_, &b_l1JetPhi);

  tt->SetBranchAddress("hltAccept", &hltAccept, &b_hltAccept);
  tt->SetBranchAddress("l1Accept", &l1Accept, &b_l1Accept);
  tt->SetBranchAddress("l1Names", &l1Names, &b_l1Names);

  int nentries = tt->GetEntries();
  std::cout << "Number of entries: " << nentries << std::endl;
  
  //book graphs and plots
  float min = 0.;
  float max = 1000.;
  int nBins = 20;

  TF1* f1 = new TF1("f1","[0]*TMath::Erf((x-[1])/[2])-[0]*TMath::Erf((-x-[1])/[2])",min,max);
  f1->SetParameters(0.5,350,40);  
  f1->FixParameter(0,0.5);
  f1->SetLineWidth(2.);
  f1->SetLineColor(kRed);

  TF1* f2 = (TF1*)f1->Clone("f2");
  f2->SetParameters(0.5,150,10);
  f2->SetLineColor(kBlack);

  //##############################################
  //##############################################

  //book graphs and plots
  TGraphErrors* totL1RateVsCut = new TGraphErrors();
  totL1RateVsCut->SetMinimum(0);
  TGraphErrors* pureL1RateVsCut = new TGraphErrors();

  TGraphErrors* totHltRateVsCut = new TGraphErrors();
  totHltRateVsCut->SetMinimum(0);
  TGraphErrors* pureHltRateVsCut = new TGraphErrors();

  //loops
  int bin = 0;
  for (int cut = 200; cut < 350; cut=cut+10)
    {
      std::cout << "analyzing point at " << cut << " GeV" << std::endl;

      int mjjHltPassed = 0;
      int excl_mjjHltPassed = 0;

      int mjjL1Passed = 0;
      int excl_mjjL1Passed = 0;

      int HTT240Passed = 0;

      for (Long64_t jentry=0; jentry<nentries;++jentry) 
	{
	  tt->GetEntry(jentry);

	  //### Sanity checks ###
	  if (l1Accept->at(HTT240) == 1)
	    ++HTT240Passed;
	  
	  //### L1 ###
	  bool l1Pass = (l1Mjj>cut &&
			 l1Jet1Pt_ > 15. &&
			 l1Jet2Pt_ > 15. &&

//.........这里部分代码省略.........
    ex[1] = 0;
    ex[2] = 0;
    ex[3] = 0;
      
    bool good=true;
      
    for(i=0;i<3;i++){
      if(avEFF[i]-avEFF[i+1]>5.)good = false;
      if(avEFF[i]==0. || avEFF[i+1]==0.) good = false;
    }
      
    if(!good){
      //cout<<Roll<<" bad scan"<<endl;
    }else{
      //cout<<Roll<<" good scan"<<endl;
    }

    //cout << x[0] <<" : " << avEFF[0] << "------" << avERR[0] << endl;
    //cout << x[1] <<" : " << avEFF[1] << "------" << avERR[1] << endl;
    //cout << x[2] <<" : " << avEFF[2] << "------" << avERR[2] << endl;
    //cout << x[3] <<" : " << avEFF[3] << "------" << avERR[3] << endl;
      
    string titlename = "Efficiency vs HV " + Roll; 
    namesave = "GlobalEff.png";
    string foldername = Roll;
    string fullname = foldername + string("/") + namesave;
    
    grp = new TGraphErrors(4, x, avEFF, ex, avERR);
    
    grp->SetTitle(titlename.c_str());
    grp->SetLineColor(kRed);
    grp->SetMarkerStyle(20);
    grp->SetMarkerSize(1.0);
    grp->SetMinimum(-0.01);
    grp->SetMaximum(100);
    grp->GetXaxis()->SetTitle("HV(kV)");
    grp->GetYaxis()->SetTitle("Efficiency(%)");
    grp->Fit(f1);

    grp->Draw("AP");
    
    double emax=f1->GetParameter(0);
    
    double parameters[3];
    parameters[0]=emax;
    parameters[1]=f1->GetParameter(1);
    parameters[2]=f1->GetParameter(2);

    for(i=0;i<3;i++){
      cout<<"i="<<parameters[i]<<" "<<endl;
    }
    
    double hv[1]={0};
    double wp=0;
    bool found=false;

    for(hv[0]=8.5;hv[0]<=9.7;hv[0]=hv[0]+0.01){
      if(emax-amano(hv[0],f1->GetParameter(1),f1->GetParameter(0),f1->GetParameter(2))<=0.05 && found==false){
	cout<<hv[0]<< " " <<amano(hv[0],f1->GetParameter(1),f1->GetParameter(0),f1->GetParameter(2))<<" "<<emax<<endl;
    	wp=hv[0];
	found=true;
      }
    }
    
    string cond = "bad";
    if(good){

//.........这里部分代码省略.........
    if(Buffer[n]=='p') Buffer[n] = '+';
    if(Buffer[n]=='m') Buffer[n] = '-';
  }
  if(!strcmp(Iso, "32"))  sprintf(Buffer, "%s^{3/2}",  Buffer, Iso);
  if(!strcmp(Iso, "p12")) sprintf(Buffer, "%s^{p1/2}", Buffer, Iso);
  PlotsRe->GetYaxis()->SetTitle(Buffer);

  //Open text file with model values for given p pi0 multipole
  sprintf(Buffer, "model/ppi0/%s.txt", Mlp);
  InModel_0 = fopen(Buffer, "r");
  //Skip two lines with table header
  fgets(Buffer, sizeof(Buffer), InModel_0);
  fgets(Buffer, sizeof(Buffer), InModel_0);
  //Read multipole model values from file
  MoPts_0 = 0;
  while(!feof(InModel_0))
  {
    if(fscanf(InModel_0, "%lf %lf %lf", &W, &Re, &Im)==3)
    {
      MoW_0[MoPts_0]  = W;
      MoRe_0[MoPts_0] = Re;
      MoIm_0[MoPts_0] = Im;
      MoPts_0++;
    }
  }
  //Close file with model values
  fclose(InModel_0);

  //Open text file with model values for given n pi+ multipole
  sprintf(Buffer, "model/npip/%s.txt", Mlp);
  InModel_p = fopen(Buffer, "r");
  //Skip two lines with table header
  fgets(Buffer, sizeof(Buffer), InModel_p);
  fgets(Buffer, sizeof(Buffer), InModel_p);
  //Read multipole model values from file
  MoPts_p = 0;
  while(!feof(InModel_p))
  {
    if(fscanf(InModel_p, "%lf %lf %lf", &W, &Re, &Im)==3)
    {
      MoW_p[MoPts_p]  = W;
      MoRe_p[MoPts_p] = Re;
      MoIm_p[MoPts_p] = Im;
      MoPts_p++;
    }
  }
  //Close file with model values
  fclose(InModel_p);

  //Create selected isospin multipole from p pi0 and n pi+ multipoles
  for(Int_t wp=0; wp<MoPts_p; wp++)
  {
    //Find corresponding energy bin between n pi+ and p pi0 multipoles
    Int_t w0;
    for(Int_t w0=0; w0<MoPts_0; w0++)
      if(MoW_p[wp]==MoW_0[w0]) break;
    //Create isospin multipoles
    if(!strcmp(Iso, "32"))  { MoRe[wp] = A_32(MoRe_0[w0], MoRe_p[wp]); MoIm[wp] = A_32(MoIm_0[w0], MoIm_p[wp]); }
    if(!strcmp(Iso, "p12")) { MoRe[wp] = A_12(MoRe_0[w0], MoRe_p[wp]); MoIm[wp] = A_12(MoIm_0[w0], MoIm_p[wp]); }
  }

  //Create graphs for real and imaginary parts of model multipole
  ModelRe = new TGraph(MoPts_p, MoW_p, MoRe);
  ModelIm = new TGraph(MoPts_p, MoW_p, MoIm);

  //Color, line size adjustments
  ModelRe->SetLineColor(kRed);
  ModelIm->SetLineColor(kBlue);
  ModelRe->SetLineWidth(2);
  ModelIm->SetLineWidth(2);

  //Set plot titles and object names
  sprintf(Buffer, "%s_%s", Mlp, Iso);
  ModelRe->SetTitle(Buffer);
  ModelIm->SetTitle(Buffer);
  sprintf(Buffer, "Model_Re%s_%s", Mlp, Iso);
  ModelRe->SetName(Buffer);
  sprintf(Buffer, "Model_Im%s_%s", Mlp, Iso);
  ModelIm->SetName(Buffer);

  //Plot graphs
  ModelRe->Draw("L"); //Plot as line, into same frame
  ModelIm->Draw("L"); //Plot as line, into same frame

  //Adjust drawing ranges for y-axis
  if((Lo==0.0) && (Hi==0.0))
  {
    PlotsRe->SetMinimum(Min*1.05);
    PlotsRe->SetMaximum(Max*1.05);
  }
  else
  {
    PlotsRe->SetMinimum(Lo);
    PlotsRe->SetMaximum(Hi);
  }

  if(SAVE) PlotsRe->SetTitle("");
  if(SAVE) sprintf(Buffer, "isospin/%s_%s.eps", Mlp, Iso);
  if(SAVE) Canvas->SaveAs(Buffer);
}

void PlotPtTresh2() {
	//  gROOT->ProcessLine(".L ./tdrstyle.C");
	//  setTDRStyle();
	Int_t r=1; 
	

	TFile *f1 = new TFile("histo_file_QCDPU50GEM2019.root");
	
	char  histo[40], histo1[40], histo2[40], histoDC[40];
	
	TH1F *hInitialBarrel = (TH1F*)f1->Get("hInitialBarrel"); 
	TH1F *hInitialEndcap = (TH1F*)f1->Get("hInitialEndcap"); 	
	
	TH1F *hPtBarrel = (TH1F*)f1->Get("hMuonPtFinal_Barrel"); 
	TH1F *hPtHE = (TH1F*)f1->Get("hMuonPtFinal_HE"); 
	TH1F *hPtFull = (TH1F*)f1->Get("hMuonPtFinal_Full"); 
	TH1F *hInitial = (TH1F*)f1->Get("entries_after_vertex"); 
	TH1F *hWP2012Pt = (TH1F*)f1->Get("hMuonPtWP2012"); 
	TH1F *h_HTot = (TH1F*)f1->Get("H_tot"); 
	TH1F *initial_entries = (TH1F*)f1->Get("initial_entries"); 
	
	cout<<"Htot"<<h_HTot->Integral()<<"  initial entries  "<<initial_entries->Integral()<<endl;
	
	int  Ninitial=hInitial->Integral();    
	cout<<" initial  "<<Ninitial<<endl;
	cout<<" initial  barrel "<<hInitialBarrel->Integral()<<endl;
	cout<<" initial endcap"<<hInitialEndcap->Integral()<<endl;
	
	cout<<" after muon barrel  "<<hPtBarrel->Integral()<<endl;
	cout<<" after muon HE  "<<hPtHE->Integral()<<endl;
	
	float MuonEff2012[1] =0;
	MuonEff2012[0]= hWP2012Pt->Integral()/hInitial->Integral();
	float Err_MuonEff2012[1];
	Err_MuonEff2012[0]=(TMath::Sqrt(Err_MuonEff2012[0]*(1-Err_MuonEff2012[0])))/Ninitial;
	//cout<<"Muon Eff WP 2012 "<<hWP2012Pt->Integral()/hInitial->Integral()<<endl;
	
	/////////////////////////compute selection eff as a function of muon pt threshold/////////////////////////////
	int EvAfterPtCut_HE[100]; float EffSel[100], xgr[100], ex[100],ey[100], EffSelFull[100], eyfull[100], Eff_Barrel[100], Eff_Endcap[100], err_Eff_barrel[100],err_Eff_endcap[100];
	int EvAfterPtCut_Barrel[100]; 
	int EvAfterPtCut[100];
	int EvAfterPtCut_BB[100];	
	int EvAfterPtCut_EE[100];	
	
	//  std::cout<<" Muon Pt in |eta|<2.4 "<<hMuonTot->Integral(0,200)<<endl;
	
	for(int i=0;i<40;i++){
		EvAfterPtCut_Barrel[i]=hPtBarrel->Integral(10,100);
		EvAfterPtCut_HE[i]=hPtHE->Integral(i+1,100);
		
		EvAfterPtCut[i]=hPtFull->Integral(i+1,100); 
		EvAfterPtCut_BB[i] =hPtBarrel->Integral(i+1,100);
		EvAfterPtCut_EE[i] =hPtHE->Integral(i+1,100);
		
		EffSelFull[i] = EvAfterPtCut[i] / (hInitialBarrel->Integral() + hInitialEndcap->Integral() );		
		Eff_Barrel[i]= EvAfterPtCut_BB[i]/hInitialBarrel->Integral();  
		Eff_Endcap[i]=EvAfterPtCut_EE[i]/hInitialEndcap->Integral();  
		
		std::cout<<i*2<<" nev after pt cut barrel "<<EvAfterPtCut_BB[i]<<"  nev after pt cut endcap  "<<EvAfterPtCut_HE[i]<<"  eff sel  "<<endl;
		
		xgr[i]=i*2;
		ex[i]=0; 
		EffSel[i]=(EvAfterPtCut_Barrel[i]+EvAfterPtCut_HE[i])/(hInitial->Integral());
		
		ey[i]=TMath::Sqrt((EffSel[i]*(1-EffSel[i])))/hInitial->Integral();
		eyfull[i]=TMath::Sqrt((EffSelFull[i]*(1-EffSelFull[i])))/hInitial->Integral();
		//  std::cout<<i*2<<" nev after pt cut barrel || endcap "<<EvAfterPtCut_Barrel[i]+EvAfterPtCut_HE[i]<<" eff sel  "<<EffSel[i]<<" errore "<<ey[i]<<endl;
	}
	
	double xgr12[1]; xgr12[0]=13.0; 
	double ex12[1]; ex12[0]=0.; 
	TGraphErrors* gr = new TGraphErrors(26,xgr,EffSel,ex,ey);
	TGraphErrors* gr_full = new TGraphErrors(26,xgr,EffSelFull,ex,eyfull);
	TGraph* gr12 = new TGraph(1);
	gr12->SetPoint(0,20,MuonEff2012[0]);
	
	TGraphErrors* gr_BB = new TGraphErrors(26,xgr,Eff_Barrel,ex,ey);	
	TGraphErrors* gr_EE = new TGraphErrors(26,xgr,Eff_Endcap,ex,ey);	
	
	///////////////////////////////////////////************************************////////////////////////////////////////////
	TCanvas *Canv5 = new TCanvas("Canv5","Palle",800,700);
	gr->SetMarkerStyle(24);  
	gr12->SetMarkerStyle(23); 
	gr->SetMarkerColor(kRed);  
	gr12->SetMarkerColor(kBlue);  
	gr->SetTitle(" Selection Efficiency as a function of muon pt threshold ");
	gr->GetYaxis()->SetTitle("Selection Efficiency");
	gr->GetYaxis()->SetTitleOffset(1.4);
	gr->GetXaxis()->SetTitle("P_{T}^{#mu Endcap} Cut (GeV/c)");
	
	gr->SetMarkerSize(1.2);  
	gr12->SetMarkerSize(1.4); 
	
	gr->Draw("ACP");
	gr12->Draw("CP");
	gr->SetMinimum(0.5);
	TLegend* legend = new TLegend( 0.15, 0.5, 0.75, 0.7);
	legend->SetTextSize(0.025);
	legend->SetFillColor(0);
	legend->AddEntry(gr, "PU=50 wGEM (P_{T}^{#mu, barrel}>20, |#eta|<2.1)","p");
//.........这里部分代码省略.........

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

   //
   // prediction w/ noise-free observables
   //

   int m = 6;
   double t [m];
   t[0] = 0.05;
   t[1] = 0.29;
   t[2] = 0.31;
   t[3] = 0.50;
   t[4] = 0.75;
   t[5] = 0.95;

   TVectorD ft(m);
   ft[0] = 1.0;
   ft[1] = 1.0;
   ft[2] = 2.0;
   ft[3] = 1.0;
   ft[4] = 1.0;
   ft[5] = 1.0;

   TMatrixD K11(m,m);
   TMatrixD K12(m,n);
   TMatrixD K21(n,m);
   TMatrixD K22(n,n);

   for(int i=0; i < n; i++){
      for(int j=0; j < n; j++){

         if( i < m and j < m ){
            K11[i][j] = kern(t[i],t[j],1);
            if( i==j ) K11[i][j] += 0.00001;
         }
         if( i < m and j < n ){
            K12[i][j] = kern(t[i],x[j],1);
            if( i==j ) K12[i][j] += 0.00001;
         }
         if( i < n and j < m ){
            K21[i][j] = kern(x[i],t[j],1);
            if( i==j ) K21[i][j] += 0.00001;
         }
         if( i < n and j < n ){
            K22[i][j] = kern(x[i],x[j],1);
            if( i==j ) K22[i][j] += 0.00001;
         }

      }
   }

   K11[2][2] += 0.1;

   // use Cholesky decomposition to obtain vector distributed as multivariate gaussian
   
   // compute covariance matrix
   TMatrixD fC(n,n);
   TMatrixD K11inv(m,m);
   K11inv = K11;
   K11inv.Invert();
   fC = K22 - K21*K11inv*K12;
   TDecompChol fD(fC);
   TMatrixD fL(n,n);
   fL = fD.GetU();

   //
   // result
   //
   
   TVectorD v(n);
   TVectorD fy(n);
   TVectorD fy_err(n);
   // mean value and error band
   for(int i=0; i < n; i++) v[i] = 1.0;
   fy = K21*K11inv*ft;
   fy_err = fL*v;

   TGraphErrors *gtest = new TGraphErrors();
   for(int i=0; i < n; i++){
      gtest->SetPoint(i, x[i], fy[i]);
      gtest->SetPointError(i, 0.0, fy_err[i]);
   }

   TCanvas * canvas = new TCanvas("canvas","canvas",800,800);
   canvas->cd();

   TGraph *gtrain = new TGraph();
   for(int i=0; i < m; i++) gtrain->SetPoint(i, t[i], ft[i]);

   gtest->SetFillStyle(3002);
   gtest->SetMinimum(-3.5);
   gtest->SetMaximum(4.5);
   gtest->Draw("AC3");
   gtrain->SetMarkerStyle(20);
   gtrain->SetMarkerColor(2);
   gtrain->Draw("P");

   canvas->Draw();

   return;
}

void make1DLimit(TString combine_dir,TString type= "WH",bool blind=true){
  //TString combine_dir = "test_runSusyHgg/signalInj_sms_ChiWH_0_175/";
  //WH
  wh_limits.push_back(wh_125);
  wh_limits.push_back(wh_150);
  wh_limits.push_back(wh_175);
  wh_limits.push_back(wh_200);
  //HH
  hh_limits.push_back(hh_125);
  hh_limits.push_back(hh_150);
  hh_limits.push_back(hh_175);
  hh_limits.push_back(hh_200);

  TGraph obser( (200-125)/25 );
  TGraph graph( (200-125)/25 );
  TGraphAsymmErrors error( (200-125)/25 );
  TGraphAsymmErrors error2S( (200-125)/25 );

  TGraph obser_r( (200-125)/25 );
  TGraph graph_r( (200-125)/25 );
  TGraphAsymmErrors error_r( (200-125)/25 );
  TGraphAsymmErrors error_r2S( (200-125)/25 );
  
  TGraphErrors* theo = 0;
  if(type=="WH") theo = getTheoXSec("xsecs/CharginoNeutralino.txt");
  else theo = getTheoXSec("xsecs/Higgsino_ElectroHiggs.txt");
  //else theo = getTheoXSec("/home/amott/HggApp/SusyHgg/xsecs/Higgsino.txt");

  for(int m=125;m<=200;m+=25) {
    int i=(m-125)/25;
    TFile limit_file(Form("%s/higgsCombineChi%s_0_%d.Asymptotic.mH120.root",combine_dir.Data(),type.Data(),m) );
    TTree *limit_tree = (TTree*)limit_file.Get("limit");
    TTreeFormula limit_form("get_limit","limit",limit_tree);
    
    float down_2s = -1;
    float down = -1;
    float exp = -1;
    float up = -1;
    float up_2s = -1;
    float obs = -1;

    if( type == "WH" )
      {
	down_2s = wh_limits.at(i)[0];
	down = wh_limits.at(i)[1];
	exp = wh_limits.at(i)[2];
	up = wh_limits.at(i)[3];
	up_2s = wh_limits.at(i)[4];
	obs = wh_limits.at(i)[5];
      }
    else if ( type == "HH")
      {
	down_2s = hh_limits.at(i)[0];
	down = hh_limits.at(i)[1];
	exp = hh_limits.at(i)[2];
	up = hh_limits.at(i)[3];
	up_2s = hh_limits.at(i)[4];
	obs = hh_limits.at(i)[5];
      }
    else
      {
	std::cerr << "UNRECOGNIZED OPTION!!! QUITTING" << std::endl;
      }
    
    if(i==0) m+=5; //first point is actually at m=130
    graph.SetPoint(i,float(m), exp);
    error.SetPoint(i,float(m), exp);
    error2S.SetPoint(i, float(m), exp);
    error.SetPointError(i, 0, 0, exp-down, up-exp);
    error2S.SetPointError(i, 0 , 0 , exp-down_2s, up_2s-exp);
  
    graph_r.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r2S.SetPoint(i,float(m),exp/theo->Eval(m));
    error_r.SetPointError(i,0,0,(exp-down)/theo->Eval(m),(up-exp)/theo->Eval(m));
    error_r2S.SetPointError(i, 0, 0, (exp-down_2s)/theo->Eval(m), (up_2s-exp)/theo->Eval(m) );
    
  
    obser.SetPoint(i,float(m),obs);
    obser_r.SetPoint(i,float(m),obs/theo->Eval(m));
    if(i==0) m-=5;
  }


    TCanvas cv;
    cv.SetLogy();
    cv.SetGrid(1,1);
    theo->SetMaximum(1e2);
    theo->SetMinimum(1e-2);
    theo->GetYaxis()->SetLabelSize(0.05);
    theo->GetYaxis()->SetTitleSize(0.06);
    theo->GetYaxis()->SetTitleOffset(0.8);
    theo->GetYaxis()->SetTitle("95% CL #sigma upper limit (pb)");
    theo->GetXaxis()->SetTitle("m_{chargino} (GeV)");
    if(type=="HH") theo->GetXaxis()->SetTitle("m_{neutralino} (GeV)");
    theo->SetFillColor(kBlue);
    theo->SetLineStyle(kDotted);
    theo->SetLineWidth(2);
    
    error.SetMaximum(1e2);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
					}
					std::string fileName = filepath + material[mat] + interaction[inter] + std::string("_QE_") + std::string(pdeChar) + std::string("_SPTR_80_ASIC_30_DT300_histos.root"); std::string(Form("%d", 1));
					std::cout << fileName << std::endl;
					TFile *fIn = new TFile(fileName.c_str(), "read");

					std::string baseName = "AVG.";
					std::string smearName = "";
					if(smear==1){
						smearName = "Smear";
					}

					//NPE
					for(int npe=0;npe<10;npe++){
						std::string histName = baseName + smearName + std::string("DTOF2") + std::string("N") + std::string(Form("%d", npes[npe]));
						std::cout << histName << std::endl;

						TH1F *h1 = (TH1F*) fIn->Get(histName.c_str());
						TF1* gauF1 = new TF1("gauF1","gaus",-200,200);
						h1->Fit("gauF1","","e",-200,200);
						h1->Draw();

						std::string path = "/home/jmbenlloch/next/petalo/work/histo/phys/5ps/jitter_avg/npe/";
						std::string histo = path + material[mat] + interaction[inter] + std::string("_QE_") + std::string(pdeChar) + std::string("_") + histName + std::string(".png");
						std::cout << histo << std::endl;

						c1->Print(histo.c_str());

						// FWHM
						sigmasNPE[mat][inter][smear][pde][npe] = 2.35 * h1->GetFunction("gauF1")->GetParameter(2);
						errorsNPE[mat][inter][smear][pde][npe] = 2.35 * h1->GetFunction("gauF1")->GetParError(2);
					}

					fIn->Close();
				}
			}
		}
	}
	c1->Close();

	double npesD[10] = {1,2,3,4,5,6,7,8,9,10};
	//EColor colors[5] = {kRed,kBlue,kGreen+3,kCyan};
	EColor colors[5] = {kRed,kBlack,kGreen+3,kBlue};

	for(int mat=0; mat<2;mat++){
		if (mat==0){
			int pdes[5] = {10,15,20,100};
		}else{
			int pdes[5] = {30,50,70,100};
		}

		for(int inter=1; inter<2;inter++){
			for(int smear=0; smear<1;smear++){
				std::string smearName = "";
				if(smear==1){
					smearName = "_smear";
				}
				//NPE
				//TCanvas *c2 = new TCanvas("c2","multipads",900,700);
				TCanvas *c2 = tcanvset("c");

				//TLegend *leg = new TLegend(0.7, 0.7, 0.9, 0.9);
				TLegend *leg = tlegset();
				leg->SetFillColor(0);

				for(int pde=0; pde<4;pde++){
					TGraphErrors *gLxe = new TGraphErrors(10, npesD, sigmasNPE[mat][inter][smear][pde], 0, errorsNPE[mat][inter][smear][pde]);
					gLxe->SetLineColor(colors[pde]);
					gLxe->SetLineWidth(2);
					gLxe->SetLineStyle(4-pde);
					grerrplot(*gLxe);

					leg->AddEntry(gLxe,Form("PDE %d%%", pdes[pde]), "lp");
					std::cout << "pde: " << pde << std::endl;

					if(pde==0){
						gLxe->SetTitle("");
						gLxe->GetXaxis()->SetTitle("Number of photoelectrons");
						gLxe->GetXaxis()->SetLimits(0,10);
						gLxe->GetYaxis()->SetTitle("CRT (ps)");
						gLxe->SetMinimum(0.);
						gLxe->SetMaximum(250.);
						gLxe->Draw();

					}else{
						gLxe->Draw("same");
					}
					if(pde==3){
						leg->Draw("same");
					}
				}
				std::string path = "/home/jmbenlloch/next/petalo/PESOA/PETALOS/root/paper/plots/";
				std::string histo = path + material[mat] + interaction[inter] + std::string("_avg_npe") + smearName + std::string(".png");

				c2->Print(histo.c_str());
				c2->Close();

			}
		}
	}
}