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

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

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

}

//.........这里部分代码省略.........
TH1F *hZPApmcUncalib = new TH1F("hZPApmcUncalib","hZPApmcUncalib",200.,0.,2000.);
TH1F *hZEM1 = new TH1F("hZEM1","hZEM1",200.,0.,2000.);
TH1F *hZEM2 = new TH1F("hZEM2","hZEM2",200.,0.,2000.);

ttree->SetBranchAddress("run",&runNumber);
ttree->SetBranchAddress("ZNC_mean_value",&ZNC_mean);
ttree->SetBranchAddress("ZNA_mean_value",&ZNA_mean);
ttree->SetBranchAddress("ZPC_mean_value",&ZPC_mean);
ttree->SetBranchAddress("ZPA_mean_value",&ZPA_mean);
ttree->SetBranchAddress("ZNCuncalib_mean",&ZNCuncalib_mean);
ttree->SetBranchAddress("ZNAuncalib_mean",&ZNAuncalib_mean);
ttree->SetBranchAddress("ZPCuncalib_mean",&ZPCuncalib_mean);
ttree->SetBranchAddress("ZPAuncalib_mean",&ZPAuncalib_mean);
ttree->SetBranchAddress("ZEM1_mean_value",&ZEM1_mean);
ttree->SetBranchAddress("ZEM2_mean_value",&ZEM2_mean);
ttree->SetBranchAddress("ZNC_X_Centroid",&ZNC_XCent);
ttree->SetBranchAddress("ZNC_Y_Centroid",&ZNC_YCent);
ttree->SetBranchAddress("ZNA_X_Centroid",&ZNA_XCent);
ttree->SetBranchAddress("ZNA_Y_Centroid",&ZNA_YCent);
ttree->SetBranchAddress("ZNC_X_Centroid_Err",&ZNC_XCent_err);
ttree->SetBranchAddress("ZNC_Y_Centroid_Err",&ZNC_YCent_err);
ttree->SetBranchAddress("ZNA_X_Centroid_Err",&ZNA_XCent_err);
ttree->SetBranchAddress("ZNA_Y_Centroid_Err",&ZNA_YCent_err);
ttree->SetBranchAddress("ZN_TDC_Sum",&ZN_TDC_Sum);
ttree->SetBranchAddress("ZN_TDC_Diff",&ZN_TDC_Diff);
ttree->SetBranchAddress("ZN_TDC_Sum_Err",&ZN_TDC_Sum_err);
ttree->SetBranchAddress("ZN_TDC_Diff_Err",&ZN_TDC_Diff_err);
//ttree->SetBranchAddress("ZNC_TDC",&ZNC_TDC);
//ttree->SetBranchAddress("ZNA_TDC",&ZNA_TDC);

printf(" branch addresses set\n");

TH1F *hznc = new TH1F("hznc","ZNC average signal",3,-1,1);
hznc->GetXaxis()->SetRangeUser(-1.,1.);
hznc->SetDrawOption("EP");
hznc->SetMarkerStyle(20);
hznc->SetMarkerColor(kRed);
hznc->SetLineColor(kRed);

TH1F *hzna = new TH1F("hzna","ZNA average signal",3,-1,1);
hzna->GetXaxis()->SetRangeUser(-1.,1.);
hzna->SetDrawOption("EP");
hzna->SetMarkerStyle(20);
hzna->SetMarkerColor(kRed);
hzna->SetLineColor(kRed);

TH1F *hzpc = new TH1F("hzpc","ZPC average signal",3,-1,1);
hzpc->GetXaxis()->SetRangeUser(-1.,1.);
hzpc->SetDrawOption("EP");
hzpc->SetMarkerStyle(20);
hzpc->SetMarkerColor(kRed);
hzpc->SetLineColor(kRed);

TH1F *hzpa = new TH1F("hzpa","ZPA average signal",3,-1,1);
hzpa->GetXaxis()->SetRangeUser(-1.,1.);
hzpa->SetDrawOption("EP");
hzpa->SetMarkerStyle(20);
hzpa->SetMarkerColor(kRed);
hzpa->SetLineColor(kRed);

TH1F *hzncUncalib = new TH1F("hzncUncalib","ZNC uncalibrated average signal",3,-1,1);
hzncUncalib->GetXaxis()->SetRangeUser(-1.,1.);
hzncUncalib->SetDrawOption("EP");
hzncUncalib->SetMarkerStyle(20);
hzncUncalib->SetMarkerColor(kAzure+10);
hzncUncalib->SetLineColor(kAzure+10);

void Getvn2D(bool usingCNTEP=0){
    TString str;
    TFile *fin;
    int nrun = GetTotalRun();
    std::cout<<"Totally we have "<<nrun<<" runs/segments!"<<std::endl;
    FillGoodRun();
    std::cout<<"Filling Good run finished!"<<std::endl;

    if(nrun<0) exit(1);

     ofstream fout, foutraw, fout1, fout2;
     int iharE=0;
     if(nhar==1||nhar==2) iharE=1;
     for(int icent=0;icent<ncent;icent++){
      for(int ihar=0;ihar<nhar;ihar++){
          if(icent!=0) continue;
       for(int isub=0;isub<nsub;isub++){
        int n = ihar+1.0+iharE;
        str = choosesub(isub);
        TString UseCNTEP;
        if(str=="ABORT") continue;
        if(usingCNTEP)
         UseCNTEP = "UseCNTEP";
        else
         UseCNTEP = "NoUseCNTEP";
        std::cout<<UseCNTEP<<std::endl;
        std::cout<<"starting doing "<<str<<" v"<<n<<" analysis!"<<std::endl;
        cout<<Form("Result/%s/res%d_%d_%s.dat",UseCNTEP.Data(),n,icent,str.Data())<<endl;
         fout1.open(Form("Result/%s/res%d_%d_%s.dat",UseCNTEP.Data(),n,icent,str.Data())); //using str as event plane detector
         fout2.open(Form("Result/%s/psi%d_%d_%s.dat",UseCNTEP.Data(),n,icent,str.Data())); //using str as event plane detector
         float reso = GetReso(icent,ihar,isub,usingCNTEP);
         fout1<<reso<<std::endl;
         if(reso<0) {std::cout<<"resolution is wrong!"<<std::endl; reso = -9999;}
        for(int irun=0;irun<nrun;irun++){
         fout2<<GetRun(irun)<<" "<<GoodRunFit[icent][ihar][isub][irun]<<std::endl;
        }
         TH2F* hvobsall = new TH2F(Form("hvobs_%d_%d_%d",icent,ihar,isub),Form("hvobs_%d_%d_%d",icent,ihar,isub),60,0,6,220,-1.1,1.1);
        // TH2F* hvobsallsq = new TH2F(Form("hvobssq_%d_%d_%d",icent,ihar,isub),Form("hvobssq_%d_%d_%d",icent,ihar,isub),60,0,6,220,-1.1,1.1);
         
        for(int iphi=0;iphi<nphi+1;iphi++){
         TH2F* hvobs = new TH2F(Form("hvobs_%d_%d_%d_%d",icent,ihar,isub,iphi),Form("hvobs_%d_%d_%d_%d",icent,ihar,isub,iphi),60,0,6,220,-1.1,1.1);
        // TH2F* hvobssq = new TH2F(Form("hvobssq_%d_%d_%d_%d",icent,ihar,isub,iphi),Form("hvobssq_%d_%d_%d_%d",icent,ihar,isub,iphi),60,0,6,220,-1.1,1.1);
         string phistr = (iphi==0)?"_east":"_west";
         if(iphi==nphi) phistr = "";
         cout<<"open v2 file"<<endl;
         fout.open(Form("Result/%s/v%d_%d%s_%s.dat",UseCNTEP.Data(),n,icent,phistr.c_str(),str.Data())); //using str as event plane detector
         cout<<"open v2raw file"<<endl;
         foutraw.open(Form("Result/%s/v%draw_%d%s_%s.dat",UseCNTEP.Data(),n,icent,phistr.c_str(),str.Data())); //using str as event plane detector
         if(iphi<nphi){
        for(int irun=0;irun<nrun;irun++){
         //std::cout<<"cent = "<<icent<<"; n = "<<n<<" ;isub = "<<str<<" ;run = "<<irun<<" "<<phistr<<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/39GeV/output/%s",GetRun(irun).Data()));
         fin = TFile::Open(Form("/gpfs/mnt/gpfs02/phenix/plhf/plhf1/xuq/phenix/flow/Run16dAu/work/39GeV/treeout/%s",GetRun(irun).Data()));
         if(!(GoodRunFit[icent][ihar][isub][irun]>0.2 && GoodRunFit[icent][ihar][isub][irun]<3.0)){
         std::cout<<"cent = "<<icent<<"; n = "<<n<<" ;isub = "<<str<<" ;run = "<<GetRun(irun)<<" is bad run!"<<std::endl;
         fin->Close();
        continue;
         }
         TH2F* hvobstemp = (TH2F*)fin->Get(Form("vobs%s_0_0_%d_%d_%d",str.Data(),icent,ihar,iphi));
         //TH2F* hvobssqtemp = (TH2F*)fin->Get(Form("vobs%ssq_%d_%d_%d",str.Data(),icent,ihar,iphi));
         hvobs->Add(hvobstemp);
         //hvobssq->Add(hvobssqtemp);
         fin->Close();
        }
         }
        hvobsall->Add(hvobs);
        //hvobsallsq->Add(hvobssq);
        if(iphi==nphi){
        hvobs = hvobsall;
        //hvobssq = hvobsallsq;
        }
            TH1F* ptProj = (TH1F*)hvobs->ProjectionX(Form("hptProj"),0,-1);
         for(int ipt=0;ipt<npt-1;ipt++){
             int xbinmin = hvobs->GetXaxis()->FindBin(ptbin[ipt]+eps);
             int xbinmax = hvobs->GetXaxis()->FindBin(ptbin[ipt+1]-eps);
           //  std::cout<<xbinmin<<" "<<xbinmax<<std::endl;
           //  std::cout<<ptbin[ipt]<<" "<<ptbin[ipt+1]<<std::endl;
            TH1F* hvobsProj = (TH1F*)hvobs->ProjectionY(Form("hvobsProj_%d",ipt),xbinmin,xbinmax);
           // TH1F* hvobssqProj = (TH1F*)hvobssq->ProjectionY(Form("hvobssqProj_%d",ipt),xbinmin,xbinmax);
            float vobs = hvobsProj->GetMean();
            float Ntracks = hvobsProj->Integral();
           // float vobssq = hvobssqProj->GetMean();
            float v = vobs/reso;
            if(Ntracks>0)
            float verr = hvobsProj->GetRMS()/reso/sqrt(Ntracks);
            else verr = -9999;
           // float verr = sqrt(vobssq/reso/reso-(v*v))/sqrt(Ntracks);
            ptProj->GetXaxis()->SetRangeUser(ptbin[ipt]+eps,ptbin[ipt+1]-eps);
            float pt = ptProj->GetMean();
            fout<<pt<<" "<<v<<" "<<verr<<" "<<std::endl;
            foutraw<<pt<<" "<<vobs<<" "<<verr*reso<<" "<<std::endl;
         }
        fout.close();
        foutraw.close();
         }
        fout1.close();
        fout2.close();
        }
        }
//.........这里部分代码省略.........

void readTree_background() {
  Char_t *filename = "Background.root";
  // Retrieve the TTree
  TFile* myFile = TFile::Open(filename);
  TTree* tree = (TTree*)(myFile->Get("tree"));
  Double_t Et1, eta1, phi1, Et2, eta2, phi2;
  tree->SetBranchAddress("Et1" ,&Et1 );
  tree->SetBranchAddress("eta1",&eta1);
  tree->SetBranchAddress("phi1",&phi1);
  tree->SetBranchAddress("Et2" ,&Et2 );
  tree->SetBranchAddress("eta2",&eta2);
  tree->SetBranchAddress("phi2",&phi2);
  // Book histograms
  TH1F* hmass = new TH1F("hmass","m_{#gamma#gamma}",60,100.,160.);
  hmass->GetXaxis()->SetTitle("Invariant mass [GeV]");
  hmass->GetYaxis()->SetTitle("Events");
  // Loop over the events
  Long64_t events = tree->GetEntries();  
  for (int i=0; i<events; i++) {
    tree->GetEntry(i);
    TLorentzVector g1,g2;
    g1.SetPtEtaPhiM(Et1,eta1,phi1,0.);
    g2.SetPtEtaPhiM(Et2,eta2,phi2,0.);
    TLorentzVector gg=g1+g2;
    hmass->Fill( gg.M() );
  }
  // Test of different background options
  hmass->DrawClone("e");
  TCanvas* myCanvas = new TCanvas("myCanvas","Background fits",800,800);
  myCanvas->Divide(2,2);
  // Linear background
  TF1* myBack1 = new TF1("myBack1","[0]+[1]*x",100.,160.);
  myBack1->SetParameter(0,events);
  myBack1->SetParameter(1,-100.);
  myCanvas->cd(1);
  hmass->Fit(myBack1);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack1);
  // Quadratic background
  TF1* myBack2 = new TF1("myBack2","[0]+[1]*x+[2]*x**2",100.,160.);
  myBack2->SetParameter(0,events);
  myBack2->SetParameter(1,-100.);
  myBack2->SetParameter(1,0.);
  myCanvas->cd(2);
  hmass->Fit(myBack2);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack2);
  // Exponential background
  TF1* myBack3 = new TF1("myBack3","[0]*exp(-x/[1])",100.,160.);
  myBack3->SetParameter(0,events);
  myBack3->SetParameter(1,100.);
  myCanvas->cd(3);
  hmass->Fit(myBack3);
  hmass->DrawClone("e");
  EvaluatePvalue(myBack3);
  // Cubic background
  TF1* myBack4 = new TF1("myBack4","[0]+[1]*x+[2]*x**2+[3]*x**3",100.,160.);
  myBack4->SetParameter(0,0.);
  myBack4->SetParameter(1,0.);
  myBack4->SetParameter(2,0.);
  myBack4->SetParameter(3,0.);
  myCanvas->cd(4);
  TFitResultPtr fit4 = hmass->Fit(myBack4,"S");
  hmass->DrawClone("e");
  EvaluatePvalue(myBack4);
}

//.........这里部分代码省略.........
   TFile *f;

   TH1D *hFake;
   TH1D *hReal;
   
   f = TFile::Open(fpath.c_str());
   
     {	
	TH1D *hFake_c = (TH1D*)f->Get(histnameFake.c_str());
	hFake_c->SetMarkerSize(0.0);

	TH1D *hReal_c = (TH1D*)f->Get(histnameReal.c_str());
	hReal_c->SetMarkerSize(0.0);

	hFake = (TH1D*)hFake_c->Clone("hFake");
	hReal = (TH1D*)hReal_c->Clone("hReal");
     }
   
   addbin(hFake);
   addbin(hReal);
   
   std::cout << hReal->Integral() << std::endl;
   std::cout << hFake->Integral() << std::endl;
   
   double iFake = hFake->Integral();
   double iReal = hReal->Integral();
//   double iAll = iFake+iReal;
   
//   hFake->Scale(1./iAll);
//   hReal->Scale(1./iAll);
   
   hFake->SetLineColor(9);
   hFake->SetFillColor(9);
   hReal->SetLineColor(46);
   hReal->SetFillColor(46);
   
   THStack *hst = new THStack();
   hst->Add(hReal);
   hst->Add(hFake);
   
   hst->Draw("hist e1");

   float max = hst->GetMaximum();
   
   hst->SetMaximum(1.3*max);
   hst->SetMinimum(0.);

   if( histnameFake == "h_muFake_pt" ) hst->GetXaxis()->SetTitle("p_{T} [GeV]");
   
   hst->GetYaxis()->SetTitle("Number of muons");
   
   leg->AddEntry(hFake,"Fake","f");
   leg->AddEntry(hReal,"Real","f");
   
   leg->Draw();
   
   c1->Print("pics/fake.eps");
   c1->Clear();

   TH1F *hRealScaled = (TH1F*)hReal->Clone("hRealScaled");
   hRealScaled->Scale(1./iReal);

   TH1F *hFakeScaled = (TH1F*)hFake->Clone("hFakeScaled");
   hFakeScaled->Scale(1./iFake);

   hFakeScaled->SetLineColor(9);
   hFakeScaled->SetFillColor(0);
   hRealScaled->SetLineColor(46);
   hRealScaled->SetFillColor(0);
   
   hRealScaled->Draw("hist e1");
   hFakeScaled->Draw("hist e1 same");

   float max1 = hRealScaled->GetMaximum();
   float max2 = hFakeScaled->GetMaximum();
   max = (max1 > max2) ? max1 : max2;
   
   hRealScaled->SetMaximum(1.3*max);
   hRealScaled->SetMinimum(0.);

   if( histnameFake == "h_muFake_pt" ) hRealScaled->GetXaxis()->SetTitle("p_{T} [GeV]");
   if( histnameFake == "h_muFake_trackerLayersWithMeasurement" ) hRealScaled->GetXaxis()->SetTitle("trackerLayersWithMeasurement");
   if( histnameFake == "h_muFake_numberOfValidMuonHits" ) hRealScaled->GetXaxis()->SetTitle("numberOfValidMuonHits");
   if( histnameFake == "h_muFake_numberOfMatches" ) hRealScaled->GetXaxis()->SetTitle("numberOfMatches");
   if( histnameFake == "h_muFake_numberOfMatchedStations" ) hRealScaled->GetXaxis()->SetTitle("numberOfMatchedStations");
   if( histnameFake == "h_muFake_numberOfValidHits" ) hRealScaled->GetXaxis()->SetTitle("numberOfValidHits");
   if( histnameFake == "h_muFake_numberOfValidPixelHits" ) hRealScaled->GetXaxis()->SetTitle("numberOfValidPixelHits");
   if( histnameFake == "h_muFake_numberOfHits" ) hRealScaled->GetXaxis()->SetTitle("numberOfHits");
   if( histnameFake == "h_muFake_normalizedChi2GlobalTrack" ) hRealScaled->GetXaxis()->SetTitle("normalizedChi2GlobalTrack");
   if( histnameFake == "h_muFake_normalizedChi2InnerTrack" ) hRealScaled->GetXaxis()->SetTitle("normalizedChi2InnerTrack");
   
   hRealScaled->GetYaxis()->SetTitle("Normalized to unity");
   
   leg->Draw();
   
   c1->Print("pics/fakeComp.eps");
   c1->Clear();
   
   gApplication->Terminate();
}

//.........这里部分代码省略.........
            if(fits2b == s2bValues[s2bStep]) {

//             tphArray[tphStep]  = fittph;
               tphArray[tphStep]  = sqrt(1.0/(1.0+fittph));

               if (fittph < tphMin) tphMin = fittph;

               if (fittph > tphMax) tphMax = fittph;

               fitxArray[tphStep] = fitx;

               if (fitx < fitxMin) fitxMin = fitx;

               if (fitx > fitxMax) fitxMax = fitx;

               tphStep++;


              TMarker *m = new TMarker(fitxArray[tphStep],tphArray[tphStep],20);
              m->SetMarkerSize(2);
              m->SetMarkerColor(31+tphStep);
              m->Draw();



            }
         }

         if (s2bStep == 0) {

            TH1F* frame = MyC->DrawFrame(0.0,0.0,1.1*fitxMax,1.1);
//          TH1F* frame = MyC->DrawFrame(0.7*fitxMin,0.7*tphMin,1.1*fitxMax,1.1*tphMax);

            TAxis *xaxis = frame->GetXaxis();
            TAxis *yaxis = frame->GetYaxis();
   
            xaxis->SetTitle("x = u^{2}/v^{2}");
            xaxis->CenterTitle();
            xaxis->SetTitleOffset(1.);
            xaxis->SetDecimals();
            xaxis->SetLabelSize(0.03);
            xaxis->SetLabelOffset(0.01);
        
            yaxis->SetTitle("tan^{2}(#phi)");
            yaxis->CenterTitle();
            yaxis->SetTitleOffset(1.);
            yaxis->SetDecimals();
            yaxis->SetLabelSize(0.03);
            yaxis->SetLabelOffset(0.01);

            frame->SetTitle(plottitle.c_str());

         }

         tphplots[s2bStep] = new TGraph(tphSteps,fitxArray,tphArray);

         tphplots[s2bStep]->SetMarkerStyle(20);

         tphplots[s2bStep]->SetMarkerSize(0.4);
   
         tphplots[s2bStep]->Draw("CP");

      }
   }

   else if (strpltmd.compare("s2b") == 0) {

void DoEvolutions( const TString &sim, Int_t time, Int_t Nbins=1, const TString &options="") { 
  
#ifdef __CINT__  
  gSystem->Load("libptools.so");
#endif

  PGlobals::Initialize();
    
  // Palettes!
  gROOT->Macro("PPalettes.C");

  TString opt = options;
  // cout << "options = " << opt << endl;

  // Load PData
  PData *pData = PData::Get(sim.Data());
  pData->LoadFileNames(time);
  if(!pData->IsInit()) return;

  Bool_t CYL = kFALSE;
  if(sim.Contains("cyl")) { CYL = kTRUE; opt += "cyl"; } 
    
  Bool_t ThreeD = kFALSE;
  if(sim.Contains("3D")) ThreeD = kTRUE; 
  
  // Some plasma constants
  Double_t n0 = pData->GetPlasmaDensity();
  Double_t kp = pData->GetPlasmaK();
  Double_t skindepth = 1.0;
  if(kp!=0.0) skindepth = 1/kp;
  Double_t E0 = pData->GetPlasmaE0();

  // Some initial beam properties:
  Float_t Ebeam = pData->GetBeamEnergy() * PUnits::MeV;
  Float_t gamma = pData->GetBeamGamma();
  Float_t vbeam = pData->GetBeamVelocity();
  
  Double_t rms0 = pData->GetBeamRmsY() * kp;
  if(CYL)  rms0 = pData->GetBeamRmsR() * kp;
  
  // Time in OU
  Float_t Time = pData->GetRealTime();
  // z start of the plasma in normalized units.
  Float_t zStartPlasma = pData->GetPlasmaStart() * kp;
  // z start of the beam in normalized units.
  Float_t zStartBeam = pData->GetBeamStart() * kp;
  
  if(opt.Contains("center")) {
    Time -= zStartPlasma;
    if(opt.Contains("comov"))      // Centers on the head of the beam.
      Time += zStartBeam;
  }

  // Beam charge 2D and 1D histogram (on-axis)
  // ------------------------------------------------------------------
  TH2F *hDen2D = NULL;
  if(pData->GetChargeFileName(1)) {
    char hName[24];
    sprintf(hName,"hDen2D");
    hDen2D = (TH2F*) gROOT->FindObject(hName);
    if(hDen2D) { delete hDen2D; hDen2D = NULL; }

    if(!ThreeD)
      hDen2D = pData->GetCharge(1,opt);
    else
      hDen2D = pData->GetCharge2DSliceZY(1,-1,1,opt+"avg");

    hDen2D->SetName(hName);
    hDen2D->GetXaxis()->CenterTitle();
    hDen2D->GetYaxis()->CenterTitle();
    hDen2D->GetZaxis()->CenterTitle();
    
    if(opt.Contains("comov"))
      hDen2D->GetXaxis()->SetTitle("k_{p}#zeta");
    else
      hDen2D->GetXaxis()->SetTitle("k_{p}z");
    
    if(CYL) 
      hDen2D->GetYaxis()->SetTitle("k_{p}r");
    else
      hDen2D->GetYaxis()->SetTitle("k_{p}y");

    hDen2D->GetZaxis()->SetTitle("n_{b}/n_{0}");

 
  }
  
  // Define ranges from the charge 2D histogram:
  // Binning for 2D histograms:
  // We get this values from the 2D density histogram.
  Int_t   x1Nbin    = hDen2D->GetNbinsX();
  Float_t x1Range   = (hDen2D->GetXaxis()->GetXmax() - hDen2D->GetXaxis()->GetXmin());
  Float_t x1Mid     = (hDen2D->GetXaxis()->GetXmax() + hDen2D->GetXaxis()->GetXmin())/2.;
  Float_t x1Min     = hDen2D->GetXaxis()->GetXmin();
  Float_t x1Max     = hDen2D->GetXaxis()->GetXmax();
  
  Int_t   x2Nbin    = hDen2D->GetNbinsY();      
  Float_t x2Range   = (hDen2D->GetYaxis()->GetXmax() - hDen2D->GetYaxis()->GetXmin());
  Float_t x2Mid     = (hDen2D->GetYaxis()->GetXmax() + hDen2D->GetYaxis()->GetXmin())/2.;
  Float_t x2Min     = x2Mid - x2Range/2;
//.........这里部分代码省略.........

void PlotShapeSystematics() {

  TFile *file = 0;
  TH1F *DefaultShape = 0;
  TH1F *UpShape = 0;
  TH1F *DownShape = 0;
  TCanvas *cv = 0;
  TLegend *legend = 0;

  //*********************************************************
  //0 Jet Bin - OF
  //*********************************************************

  file = new TFile("/data/smurf/sixie/data/Thesis/cards/130/hwwof_0j.input.root","READ");

  DefaultShape = (TH1F*)file->Get("histo_qqWW");
  UpShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwUp");
  DownShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwDown");

  assert(DefaultShape);
  assert(UpShape);
  assert(DownShape);

  NormalizeHist(DefaultShape);
  NormalizeHist(UpShape);
  NormalizeHist(DownShape);


  cv = new TCanvas("cv","cv",800,600);

  legend = new TLegend(0.2, 0.6, 0.5, 0.9);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);
  legend->SetFillColor(kWhite);
  legend->AddEntry(DefaultShape, "Default Shape", "L");
  legend->AddEntry(UpShape, "Bounding Shape (Up)", "L");
  legend->AddEntry(DownShape, "Bounding Shape (Down)", "L");

  DefaultShape->SetLineColor(kBlack);
  UpShape->SetLineColor(kBlue);
  DownShape->SetLineColor(kRed);
  DefaultShape->SetLineWidth(2);
  UpShape->SetLineWidth(2);
  DownShape->SetLineWidth(2);
  DefaultShape->SetTitle("");
  DefaultShape->GetXaxis()->SetTitle("MVA discriminator");
  DefaultShape->GetYaxis()->SetTitle("Fraction of Events");
  DefaultShape->GetYaxis()->SetTitleOffset(1.4);
  DefaultShape->GetXaxis()->SetRangeUser(-1.0,1.0);
  DefaultShape->GetYaxis()->SetRangeUser(0.0,0.25);


  DefaultShape->Draw("hist");
  UpShape->Draw("same,hist");
  DownShape->Draw("same,hist");
  legend->Draw();

  cv->SaveAs("WWBkgShapeVariation_MCAtNLOScaleVariation_OF0Jet.png");
  cv->SaveAs("WWBkgShapeVariation_MCAtNLOScaleVariation_OF0Jet.eps");



  //*********************************************************
  //0 Jet Bin - SF
  //*********************************************************

  file = new TFile("/data/smurf/sixie/data/Thesis/cards/130/hwwsf_0j.input.root","READ");

  DefaultShape = (TH1F*)file->Get("histo_qqWW");
  UpShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwUp");
  DownShape = (TH1F*)file->Get("histo_qqWW_CMS_MVAWWNLOBounding_hwwDown");

  assert(DefaultShape);
  assert(UpShape);
  assert(DownShape);

  NormalizeHist(DefaultShape);
  NormalizeHist(UpShape);
  NormalizeHist(DownShape);


  cv = new TCanvas("cv","cv",800,600);

  legend = new TLegend(0.2, 0.6, 0.5, 0.9);
  legend->SetTextSize(0.04);
  legend->SetBorderSize(0);
  legend->SetFillColor(kWhite);
  legend->AddEntry(DefaultShape, "Default Shape", "L");
  legend->AddEntry(UpShape, "Bounding Shape (Up)", "L");
  legend->AddEntry(DownShape, "Bounding Shape (Down)", "L");

  DefaultShape->SetLineColor(kBlack);
  UpShape->SetLineColor(kBlue);
  DownShape->SetLineColor(kRed);
  DefaultShape->SetLineWidth(2);
  UpShape->SetLineWidth(2);
  DownShape->SetLineWidth(2);
  DefaultShape->SetTitle("");
  DefaultShape->GetXaxis()->SetTitle("MVA discriminator");
  DefaultShape->GetYaxis()->SetTitle("Fraction of Events");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
   Double_t Graph0_fehy3021[10] = {
   5.988963,
   6.218694,
   6.307709,
   6.368316,
   5.847823,
   5.655323,
   5.083818,
   4.561256,
   3.722445,
   2.542957};
   TGraphAsymmErrors *grae = new TGraphAsymmErrors(10,Graph0_fx3021,Graph0_fy3021,Graph0_felx3021,Graph0_fehx3021,Graph0_fely3021,Graph0_fehy3021);
   grae->SetName("Graph0");
   grae->SetTitle("Graph");

   Int_t ci;      // for color index setting
   TColor *color; // for color definition with alpha
   ci = TColor::GetColor("#ffff00");
   grae->SetFillColor(ci);

   ci = TColor::GetColor("#ff0000");
   grae->SetLineColor(ci);
   grae->SetLineWidth(4);
   grae->SetMarkerStyle(20);
   grae->SetMarkerSize(0);
   
   TH1F *Graph_Graph3021 = new TH1F("Graph_Graph3021","Graph",100,-2.5,2.5);
   Graph_Graph3021->SetMinimum(0);
   Graph_Graph3021->SetMaximum(140);
   Graph_Graph3021->SetDirectory(0);
   Graph_Graph3021->SetStats(0);
   Graph_Graph3021->SetLineStyle(0);
   Graph_Graph3021->SetMarkerStyle(20);
   Graph_Graph3021->GetXaxis()->SetNdivisions(505);
   Graph_Graph3021->GetXaxis()->SetLabelFont(42);
   Graph_Graph3021->GetXaxis()->SetLabelOffset(0.007);
   Graph_Graph3021->GetXaxis()->SetLabelSize(0);
   Graph_Graph3021->GetXaxis()->SetTitleSize(0.07142857);
   Graph_Graph3021->GetXaxis()->SetTitleOffset(1.1);
   Graph_Graph3021->GetXaxis()->SetTitleFont(42);
   Graph_Graph3021->GetYaxis()->SetTitle("d#sigma (W^{+}#rightarrow#font[12]{l}^{+}#nu) / d#eta_{lab} [nb]");
   Graph_Graph3021->GetYaxis()->SetLabelFont(42);
   Graph_Graph3021->GetYaxis()->SetLabelOffset(0.007);
   Graph_Graph3021->GetYaxis()->SetLabelSize(0.07142857);
   Graph_Graph3021->GetYaxis()->SetTitleSize(0.07142857);
   Graph_Graph3021->GetYaxis()->SetTitleOffset(1.05);
   Graph_Graph3021->GetYaxis()->SetTitleFont(42);
   Graph_Graph3021->GetZaxis()->SetLabelFont(42);
   Graph_Graph3021->GetZaxis()->SetLabelOffset(0.007);
   Graph_Graph3021->GetZaxis()->SetLabelSize(0.05);
   Graph_Graph3021->GetZaxis()->SetTitleSize(0.06);
   Graph_Graph3021->GetZaxis()->SetTitleFont(42);
   grae->SetHistogram(Graph_Graph3021);
   
   grae->Draw("a2");
   
   Double_t Graph1_fx3022[10] = {
   2.2,
   1.75,
   1.25,
   0.75,
   0.25,
   -0.25,
   -0.75,
   -1.25,
   -1.75,

//------------------------------------------------------------------------------
//Subtraction 
//------------------------------------------------------------------------------
void Subtraction(TString  hname,
		   TString  xtitle,
		   Int_t    ngroup       = -1,
		   Int_t    precision    = 1,
		   TString  units        = "NULL",
		   Double_t xmin         = -999,
		   Double_t xmax         =  999,
		   Bool_t   moveOverflow = true)
{
  TCanvas* canvas = new TCanvas(hname, hname, 800, 800);

  TPad* pad1 = new TPad("pad1", "pad1", 0, 0.0, 1, 1.0);

  pad1->SetTopMargin   (0.08);
  //pad1->SetBottomMargin(0.02);
  pad1->Draw();
      

  //----------------------------------------------------------------------------
  // pad1
  //----------------------------------------------------------------------------
  pad1->cd();

  pad1->SetLogy(_setLogy);

  TH1F* hist[nProcesses];

  for (UInt_t ip=0; ip<nProcesses; ip++) {

    hist[ip] = (TH1F*)input[ip]->Get(hname);
    hist[ip]->SetName(hname + process[ip]);

    if (moveOverflow) MoveOverflowBins  (hist[ip], xmin, xmax);
    else              ZeroOutOfRangeBins(hist[ip], xmin, xmax);

	if (ngroup > 0) hist[ip]->Rebin(ngroup);

	if (_dataDriven && ip == iWW)    hist[ip]->Scale(WWScale[_njet]);
	if (_dataDriven && ip == iDY)    hist[ip]->Scale(ZjScale[_njet]);
	if (_dataDriven && ip == iDYtau) hist[ip]->Scale(ZjScale[_njet]);

  }

  // Data subtraction for Top background estimation
  //----------------------------------------------------------------------------
  TH1F* subData = (TH1F*)hist[iData]->Clone("subData");
  for (UInt_t ip=0; ip<nProcesses; ip++) {
	  if (ip == itt) continue;
	  if (ip == itW) continue;
	  if (ip == iData ) continue;
	  subData->Add(hist[ip],-1);
  }
  subData->SetLineColor(kRed+1);
  Double_t subData_Yield = subData->Integral();
  //subData->SetLineColor();

  // Top background
  //----------------------------------------------------------------------------
  TH1F* Top = (TH1F*)hist[itt]->Clone("Top");
  Top->Add(hist[itW]);
  Top->SetLineColor(kBlue+1);
  Double_t Top_Yield = Top->Integral();

  // Axis labels
  //----------------------------------------------------------------------------
  TAxis* xaxis = subData->GetXaxis();
  TAxis* yaxis = subData->GetYaxis();

  TString ytitle = Form("entries / %s.%df", "%", precision);

  xaxis->SetTitle(xtitle);
  yaxis->SetTitle(Form(ytitle.Data(), subData->GetBinWidth(0)));
  yaxis->SetTitleOffset(1.6);

  if (!units.Contains("NULL")) {
    
    xaxis->SetTitle(Form("%s [%s]", xaxis->GetTitle(), units.Data()));
    yaxis->SetTitle(Form("%s %s",   yaxis->GetTitle(), units.Data()));
  }


  // Draw
  //----------------------------------------------------------------------------
  xaxis->SetRangeUser(xmin, xmax);

  subData->Draw("hist");
  Top->Draw("hist same");

  // Adjust scale
  //----------------------------------------------------------------------------
  subData->SetMinimum(0.0);
  
  Float_t theMax   = GetMaximumIncludingErrors(subData, xmin, xmax);
  Float_t theMaxMC = GetMaximumIncludingErrors(Top,       xmin, xmax);

  if (theMaxMC > theMax) theMax = theMaxMC;

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

//.........这里部分代码省略.........
  TH1* ewkVV = (TH1*)f->Get("EWK_VVx");
  ewktau->Add(ewkDY);
  ewktau->Add(ewkVV);
  
  //TH1* qcd = (TH1*)f->Get("QCDInv");
  TH1* qcd = (TH1*)f->Get("QCD");
  ci = TColor::GetColor("#ffcc33");
  qcd->SetFillColor(ci);
  qcd->SetLineWidth(0);
  TH1* fakett = (TH1*)f->Get("fake_tt");
  ci = TColor::GetColor("#669900");
  fakett->SetFillColor(ci);
  fakett->SetLineWidth(0);
  TH1* fakeW = (TH1*)f->Get("fake_W");
  ci = TColor::GetColor("#cc3300");
  fakeW->SetFillColor(ci);
  fakeW->SetLineWidth(0);
  TH1* faket = (TH1*)f->Get("fake_t");

  TH1F *hFrame = new TH1F("hFrame","",20,0,400);
  hFrame->SetMinimum(ymin);
  if (log)
    hFrame->SetMaximum(ymax*1.5);
  else
    hFrame->SetMaximum(ymax);
  hFrame->SetDirectory(0);
  hFrame->SetStats(0);
  hFrame->SetLineStyle(0);
  hFrame->SetMarkerStyle(20);
  hFrame->SetXTitle("Transverse mass (#tau jet, E_{T}^{miss}), (GeV/c^{2})");
  if (paperStatus)
    hFrame->SetXTitle("Transverse mass (#tau_{h}, E_{T}^{miss}), (GeV/c^{2})");
  hFrame->SetYTitle("Events / 20 GeV/c^{2}");
  hFrame->GetXaxis()->SetTitleSize(0);
  hFrame->GetXaxis()->SetLabelSize(0);
  hFrame->GetYaxis()->SetTitleFont(43);
  hFrame->GetYaxis()->SetTitleSize(27);
  hFrame->GetYaxis()->SetTitleOffset(1.3);
  

  // signal
  hh->Scale(br*br);
  hw->Scale(2*br*(1.0-br));
  TH1* signal = (TH1*)hh->Clone();
  signal->Add(hw);

  ci = TColor::GetColor("#ff3399");
  signal->SetLineColor(ci);
  signal->SetLineStyle(2);
  signal->SetLineWidth(2);

  // Fakes
  TH1* fakes = (TH1*)(fakett->Clone());
  fakes->Add(fakeW);
  fakes->Add(faket);

  // stacked backgrounds
  THStack *exp = new THStack();
  exp->SetName("exp");
  exp->SetTitle("exp");
  exp->Add(fakes);
  exp->Add(ewktau);
  exp->Add(qcd);
  exp->Add(signal);
  
  TH1* hExpBkg = (TH1*)fakes->Clone();

//.........这里部分代码省略.........
  allmc->SetFillStyle  (   3345);
  allmc->SetLineColor  (kGray+2);
  allmc->SetMarkerColor(kGray+2);
  allmc->SetMarkerSize (      0);

  for (UInt_t ibin=1; ibin<=allmc->GetNbinsX(); ibin++) {

    Double_t binValue = 0;
    Double_t binError = 0;

    for (UInt_t ip=0; ip<nProcesses; ip++) {

      if (ip == iData) continue;

      Double_t binContent = hist[ip]->GetBinContent(ibin);
      
      binValue += binContent;
      binError += (hist[ip]->GetBinError(ibin) * hist[ip]->GetBinError(ibin));

      //We need to calculate systematic uncertainty for ggH case
//      if (_dataDriven)
//	binError += (systError[ip]*binContent * systError[ip]*binContent);
    }
    
    binError = sqrt(binError);

    allmc->SetBinContent(ibin, binValue);
    allmc->SetBinError  (ibin, binError);
  }


  // Axis labels
  //----------------------------------------------------------------------------
  TAxis* xaxis = hist[iData]->GetXaxis();
  TAxis* yaxis = hist[iData]->GetYaxis();

  TString ytitle = Form("entries / %s.%df", "%", precision);

  xaxis->SetTitle(xtitle);
  yaxis->SetTitle(Form(ytitle.Data(), hist[iData]->GetBinWidth(0)));
  yaxis->SetTitleOffset(1.6);

  if (!units.Contains("NULL")) {
    
    xaxis->SetTitle(Form("%s [%s]", xaxis->GetTitle(), units.Data()));
    yaxis->SetTitle(Form("%s %s",   yaxis->GetTitle(), units.Data()));
  }


  // Draw
  //----------------------------------------------------------------------------
  xaxis->SetRangeUser(xmin, xmax);

  hist[iData]->Draw("ep");
  hstack     ->Draw("hist,same");
  allmc      ->Draw("e2,same");
  hist[iData]->Draw("ep,same");


  // Adjust scale
  //----------------------------------------------------------------------------
  Float_t theMax   = GetMaximumIncludingErrors(hist[iData], xmin, xmax);
  Float_t theMaxMC = GetMaximumIncludingErrors(allmc,       xmin, xmax);

  if (theMaxMC > theMax) theMax = theMaxMC;

//------------------------------------------------------------------------
void DrawCFD()
{
  Int_t npeaks = 20;
  Int_t sigma=3.;
  Bool_t down=false;
  Int_t index[20];
   
  TCanvas *c1 = new TCanvas("c1", "CFD C side",0,48,1280,951);
  c1->Divide(4,3);
  gStyle->SetOptFit(1111);
  //c1->Divide(2,2);
  Char_t buf1[10];
  for (Int_t i=0; i<12; i++)
    {
      c1->cd(i+1);
      sprintf(buf1,"T0_C_%i_CFD",i+1);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      //    cfd->Draw();
      
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.05);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0){
	Float_t *xpeak = s->GetPositionX();
  	TMath::Sort(nfound, xpeak, index,down);
	Float_t xp = xpeak[index[0]];
	cout<<" index[0] "<<index[0]<<endl;
		//	Float_t xp = xpeak[1];
	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
	Float_t yp = cfd->GetBinContent(xbin);
	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
	Float_t hmax = xp+3*sigma;
	Float_t hmin = xp-3*sigma;
	cout<<hmin<< " "<<hmax<<endl;
	cfd->GetXaxis()->SetRange(hmin-10,hmax+10);
	cout<<" cfd range "<<hmin-10<<" "<<hmax+10<<endl;
	cfd->GetXaxis()->SetLabelSize(0.03);
	TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
		cfd->Fit("g1","R");
		//	cfd->Draw();
	
      }

    }
  TCanvas *c2 = new TCanvas("c2", "CFD A side",0,48,1280,951);
  c2->Divide(4,3);
  gStyle->SetOptFit(1111);
  //c1->Divide(2,2);
  Char_t buf1[10];
  for (Int_t i=0; i<12; i++)
    {
      c2->cd(i+1);
      sprintf(buf1,"T0_A_%i_CFD",i+1);
      TH1F *cfd = (TH1F*) gFile->Get(buf1);
      //    cfd->Draw();
      
      TSpectrum *s = new TSpectrum(2*npeaks,1);
      Int_t nfound = s->Search(cfd,sigma," ",0.05);
      cout<<"Found "<<nfound<<" peaks sigma "<<sigma<<endl;;
      if(nfound!=0){
	Float_t *xpeak = s->GetPositionX();
  	TMath::Sort(nfound, xpeak, index,down);
	Float_t xp = xpeak[index[0]];
	cout<<" index[0] "<<index[0]<<endl;
		//	Float_t xp = xpeak[1];
	Int_t xbin = cfd->GetXaxis()->FindBin(xp);
	Float_t yp = cfd->GetBinContent(xbin);
	cout<<"xbin = "<<xbin<<"\txpeak = "<<xpeak[1]<<"\typeak = "<<yp<<endl;
	Float_t hmax = xp+3*sigma;
	Float_t hmin = xp-3*sigma;
	cout<<hmin<< " "<<hmax<<endl;
	cfd->GetXaxis()->SetRange(hmin-10,hmax+10);
	cout<<" cfd range "<<hmin-10<<" "<<hmax+10<<endl;
	cfd->GetXaxis()->SetLabelSize(0.03);
	TF1 *g1 = new TF1("g1", "gaus", hmin, hmax);
		cfd->Fit("g1","R");
		//	cfd->Draw();
	
      }

    }

}

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

}

//.........这里部分代码省略.........
       << delim << setw(left) << Form("%.3f %s %.4f",effdata1,pm,effdata1err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effdata2,pm,effdata2err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effdata3,pm,effdata3err) << setw(4) 
       << delim << setw(left) << Form("%.3f %s %.4f",effdata4,pm,effdata4err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effdata5,pm,effdata5err) << setw(4) 
       << delimend << endl;

  cout << delimstart << setw(10) << "mc" << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effmc1,pm,effmc1err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effmc2,pm,effmc2err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effmc3,pm,effmc3err) << setw(4) 
       << delim << setw(left) << Form("%.3f %s %.4f",effmc4,pm,effmc4err) << setw(4)
       << delim << setw(left) << Form("%.3f %s %.4f",effmc5,pm,effmc5err) << setw(4) 
       << delimend << endl;

  cout << delimstart << setw(10) << "data/mc" << setw(4)
       << delim << setw(left) << Form("%.2f %s %.2f",ratio1,pm,ratio1err) << setw(4)
       << delim << setw(left) << Form("%.2f %s %.2f",ratio2,pm,ratio2err) << setw(4)
       << delim << setw(left) << Form("%.2f %s %.2f",ratio3,pm,ratio3err) << setw(4) 
       << delim << setw(left) << Form("%.2f %s %.2f",ratio4,pm,ratio4err) << setw(4)
       << delim << setw(left) << Form("%.2f %s %.2f",ratio5,pm,ratio5err) << setw(4) 
       << delimend << endl;

  //--------------------------------------
  // draw stuff
  //--------------------------------------

  hdata->Scale(1.0/hdata->Integral());
  hmc->Scale(1.0/hmc->Integral());

  if( log ) hmc->GetYaxis()->SetRangeUser(0.0001,5);  
  else      hmc->GetYaxis()->SetRangeUser(0.0,1);  

  hmc->GetXaxis()->SetTitle(xtitle);
  hmc->SetLineColor(2);
  hmc->SetMarkerColor(2);
  hmc->DrawNormalized("hist");
  hmc->DrawNormalized("sameE1");
  hdata->SetLineColor(4);
  hdata->SetMarkerColor(4);
  hdata->Draw("sameE1");

  grdata->SetLineColor(6);
  grmc->SetLineColor(7);
  //grdata->Draw("sameP");
  //grmc->Draw("sameP");

  TLegend *leg = new TLegend(0.6,0.7,0.8,0.9);
  leg->AddEntry(hdata , "data" , "lp");
  leg->AddEntry(hmc   , "MC"   , "lp");
  leg->SetBorderSize(0);
  leg->SetFillColor(0);			       
  leg->Draw();

  TLatex *t = new TLatex();
  t->SetNDC();

  if( TString(plottitle).Contains("el") ) t->DrawLatex(0.6,0.6,"electrons");
  if( TString(plottitle).Contains("mu") ) t->DrawLatex(0.6,0.6,"muons");

  if( TString(plottitle).Contains("0j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 0");
  if( TString(plottitle).Contains("1j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 1");
  if( TString(plottitle).Contains("2j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 2");
  if( TString(plottitle).Contains("3j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 3");
  if( TString(plottitle).Contains("4j") ) t->DrawLatex(0.6,0.5,"n_{jets} #geq 4");