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

//.........这里部分代码省略.........
  jetpt+= "&&pt_3rd_leadJet>";
  jetpt+= min3rdJetPt;
  jetpt+= "&&";
  cutsSig += jetpt;
  cutsSL += jetpt;
  cutsLDP += jetpt;


  TH2F* h_susy_sig[10] ;
  TH2F* h_susy_sl[10] ;
  TH2F* h_susy_ldp[10] ;
  for ( int bi=0; bi<nBinsBjets; bi++ ) {
     char hname[1000] ;
     sprintf( hname, "h_susy_sig_%db", bi+1 ) ;
     h_susy_sig[bi] = new TH2F( hname, hname, nBinsMET, Mbins, nBinsHT, Hbins ) ;
     h_susy_sig[bi] -> Sumw2() ;
     sprintf( hname, "h_susy_sl_%db", bi+1 ) ;
     h_susy_sl[bi] = new TH2F( hname, hname, nBinsMET, Mbins, nBinsHT, Hbins ) ;
     h_susy_sl[bi] -> Sumw2() ;
     sprintf( hname, "h_susy_ldp_%db", bi+1 ) ;
     h_susy_ldp[bi] = new TH2F( hname, hname, nBinsMET, Mbins, nBinsHT, Hbins ) ;
     h_susy_ldp[bi] -> Sumw2() ;
  }


  stringstream njcut ; njcut << nJetsCut;
  TString cutsNjets = "&&nJets>=";
  cutsNjets += njcut.str();

  float xsec = -1.;
  float xsec8TeV = -1. ;
  for ( int mGl = minGlMass ; mGl < maxGlMass ; mGl = mGl + 100 ) {

    int theBin = gluinoxsec->FindBin( mGl ) ;					      
    if ( theBin <=0 || theBin > gluinoxsec->GetNbinsX() ) {			      
       printf("\n\n *** can't find bin for mgl=%d.  Returned %d\n\n", mGl, theBin ) ; 
       return ; 								      
    }										      
    xsec = gluinoxsec->GetBinContent( theBin ) ;				      

    int theBin8TeV = gluinoxsec8TeV->FindBin( mGl ) ;					      
    if ( theBin8TeV <=0 || theBin8TeV > gluinoxsec8TeV->GetNbinsX() ) {			      
       printf("\n\n *** can't find bin for mgl=%d.  Returned %d\n\n", mGl, theBin ) ; 
       return ; 								      
    }										      
    xsec8TeV = gluinoxsec8TeV->GetBinContent( theBin8TeV ) ;				      

    printf("\n\n  SUSY Xsecs:  7 TeV = %f,  8 TeV = %f\n\n", xsec, xsec8TeV ) ;

    ////// for ( int mLsp = 50 ; mLsp < ( mGl - 25 ) ; mLsp = mLsp + 25 ) {
    for ( int mLsp = 300 ; mLsp < 710 ; mLsp = mLsp + 400 ) {

      inFile << mGl << " " << mLsp << " " << dummyEvts << " " ;
      printf(" mGl=%4d, mLsp=%4d\n", mGl, mLsp ) ; cout << flush ;


      printf("\n\n") ;
      for (int k = 0 ; k < nBinsBjets ; k++) {

         TString cutSMS = "mgluino>";
         cutSMS += mGl-1;
         cutSMS += "&&mgluino<";
         cutSMS += mGl+1;
         cutSMS += "&&mlsp>";
         cutSMS += mLsp-1;
         cutSMS += "&&mlsp<";

//.........这里部分代码省略.........
  EWK2 ->Add(EWK1 );
//EWK3 ->Add(EWK2 );
//EWK  ->Add(EWK3 );
  EWK  ->Add(EWK2 );
//  ttbar->Add(EWK  );
  Ztt  ->Add(EWK);
  //if(log){
//#ifdef MSSM
 //   ggH->Add(bbH);
//#else
 //   qqH->Add(VH );
  //  ggH->Add(qqH);
//#endif
 // }
  //else{
//#ifdef MSSM    
 //   bbH->Add(Ztt);
  //  ggH->Add(bbH);
//#else
 //   VH ->Add(Ztt);
  //  qqH->Add(VH );
   // ggH->Add(qqH);
//#endif
 // }

  /*
    Mass plot before and after fit
  */
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);

  canv->cd();
  if(log){ canv->SetLogy(1); }
#if defined MSSM
  if(!log){ data->GetXaxis()->SetRange(200, data->FindBin(UPPER_EDGE)); } else{ data->GetXaxis()->SetRange(200, data->FindBin(UPPER_EDGE)); };
#else
  data->GetXaxis()->SetRange(200, data->FindBin(UPPER_EDGE));
#endif
  data->SetNdivisions(505);
  data->SetMinimum(min);
  data->SetMaximum(max>0 ? max : std::max(std::max(maximum(data, log), maximum(Ztt, log)), maximum(ggHTohhTo2Tau2B, log)));
  data->Draw("e");

  TH1F* errorBand = (TH1F*)Ztt ->Clone();
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(13);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
    if(errorBand->GetBinContent(idx)>0){
      std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
      break;
    }
  }
  if(log){
    Ztt  ->Draw("histsame");
//    ttbar->Draw("histsame");
    EWK  ->Draw("histsame");
    Fakes->Draw("histsame");
    ttbar->Draw("histsame");
//#ifdef MSSM
//    VH_SM125->Draw("histsame");
//#endif
    $DRAW_ERROR
    ggHTohhTo2Tau2B  ->Draw("histsame");
  }
  else{

void AlignTopAndBottomTOFs(string ArrayPart="NV",Float_t Start=-5,Float_t End=25){
  
  for (int bar=1;bar<13;bar++){
    stringstream nameTop;
    stringstream nameBottom;

    nameTop<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"_TopCutTOFPr";
    TH1F * CurrentTop = (TH1F*)gDirectory->Get(nameTop.str().c_str());


    nameBottom<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"_BottomCutTOFPr";
    TH1F * CurrentBottom = (TH1F*)gDirectory->Get(nameBottom.str().c_str());
  
    // TCanvas *c=new TCanvas("c");
    // c->cd(1);
    // CurrentTop->Draw();
    // CurrentBottom->Draw("same");

    int binsTop = CurrentTop->GetNbinsX();

    int binsBottom= CurrentBottom->GetNbinsX();
  
    if ( binsTop !=binsBottom){
      cout<<"Nubmer of bins for the top channel not the same as for the bottom channel"<<endl;
      cout<<"This was for "<<CurrentTop->GetName()<<" and "<<CurrentBottom->GetName()<<endl;
      return;
    }
  
    int bins=binsTop;
  
    double topLow = CurrentTop->GetBinLowEdge(1);//ROOT histograms start at bin 1. bin 0 is underflow
    double topHigh= CurrentBottom->GetBinLowEdge(binsTop) +CurrentBottom->GetBinWidth(binsTop);//ROOT histograms end at bin N. N+1 is overflow bin.

    double NanoSecsPerBin = ((topHigh-topLow)*4)/binsTop;
  
    //    cout<<topLow<<" "<<topHigh<<endl;
    //    cout<<"NanoSecsPerBin "<<NanoSecsPerBin<<endl;
  
 
    int ForthOfBins = TMath::Floor(0.25* bins);
    int EighthOfBins = TMath::Floor( (1.0/8)* bins);

    int NumberOfShifts = 2*EighthOfBins; // 1/8 of the number of bins *2 for left/right shifts
  
    int TopZeroBin = CurrentTop->GetMaximumBin();
  
    int StartBin = CurrentTop->FindBin(Start);
    int EndBin = CurrentTop->FindBin(End);

    //Find shift That moves Bottom -> Top
    vector <double> TheChi2s;
    vector <int> TheShifts;
    for (int i=0;i<NumberOfShifts;i++){
      double binShift = (i - NumberOfShifts/2);

      double chi2=0;
      // cout<<"Looking in bin range "<<TopZeroBin-EighthOfBins<<" "<<TopZeroBin+EighthOfBins<<endl;
      // cout<<"That is from "<<CurrentTop->GetBinCenter(TopZeroBin-EighthOfBins)<<" "<<CurrentTop->GetBinCenter(TopZeroBin+EighthOfBins)<<endl;

      for (int bin=StartBin ;bin<EndBin;bin++){
	double b = CurrentBottom->GetBinContent(bin + binShift);
	double t = CurrentTop->GetBinContent(bin);
	if (b !=0 && t!=0){
	  double temp =((b-t)*(b-t))/( TMath::Sqrt(t) );
	  chi2+=temp;
	}
      }
    
      TheChi2s.push_back(chi2);
      TheShifts.push_back(binShift);

    }
  
    ///Find miminum by linear search 
    double min = 9999999999.0;
    int MinSpot=-1;
    for (int i=0;i<NumberOfShifts;i++){
      if ( TheChi2s[i] < min){
	min=TheChi2s[i];
	MinSpot=i;
      }
    }

    cout<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"B  "<<"1  0  "<<fixed<<setw(6)<<setprecision(4)<<TheShifts[MinSpot]/4.0*NanoSecsPerBin<<endl;
    cout<<ArrayPart<<setfill('0')<<setw(2)<<bar<<"T  "<<"1  0  0"<<endl;

  }


  cout<<"There were "<<NanoSecsPerBin<<" nano secs per bin"<<endl;
  // TGraph * graph = new TGraph();
  
  // cout<<"Size of Chi2s "<<TheChi2s.size()<<endl;
  // cout<<"Num "<<NumberOfShifts<<endl;

  // for (int i=0;i<NumberOfShifts;i++){
  //   graph->SetPoint(i,TheShifts[i],TheChi2s[i]);
  // }

  // TCanvas *c2 = new TCanvas("c2");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    
    
    
    
    Printf("SSD - QA");
    

    MPVL5=0;
    MPVErrL5=0;
    MPVL6=0;
    MPVErrL6=0;
    ChargeRatioL5=0;
    ChargeRatioErrL5=0;
    ChargeRatioL6=0;
    ChargeRatioErrL6=0;
    EmptyModulesSDD=0;
    
    
    TH2F* QAchargeRatio=(TH2F*)SSDList->FindObject("QAChargeRatio");
    
    if(QAchargeRatio->GetEntries()==0){
        printf("Run %d QAchargeRatio EMPTY -> Return\n",iRun);
    }
    
    TH2F* QAcharge=(TH2F*)SSDList->FindObject("QACharge");
    
    if(QAcharge->GetEntries()==0){
        printf("Run %d QAcharge EMPTY -> Return\n",iRun);
        
    }
    
    if((QAcharge)&&(QAchargeRatio)&&(QAcharge->GetEntries()>10)&&(QAchargeRatio->GetEntries()>10)){

        Int_t biny = QAcharge->GetXaxis()->FindBin(747);
        Int_t maxy = QAcharge->GetXaxis()->GetXmax();

        Int_t  contEmpty=0;
        Int_t  contFull=0;

        TH1D *hChargeL5=QAcharge->ProjectionY("hChargeL5",0,biny);
        TH1D *hChargeL6=QAcharge->ProjectionY("hChargeL6",biny,maxy);



        TH1D *hChargeRatioL5=QAchargeRatio->ProjectionY("hChargeRatioL5",0,biny);
        TH1D *hChargeRatioL6=QAchargeRatio->ProjectionY("hChargeRatioL6",biny,maxy);


        if(QAcharge->GetEntries()< 45000)
            contEmpty=1;

        else{
            for(Int_t i =0;i<1698;i++){

                TString tmpQ("Q");
                tmpQ+=i;

                TH1D* fHist1DQ= QAcharge->ProjectionY(tmpQ,i+1,i+1);
                Double_t mean=fHist1DQ->GetMean();

                if(TMath::Abs(mean)<1.0 ||fHist1DQ->GetEntries()<10)
                    contEmpty++;

                else
                    contFull++;

//.........这里部分代码省略.........
    // Select events with a generated photon in mid-rapidity
    bool genPhotonFlag=false;
    if ( isMC )  
      genPhotonFlag = true;
    else {
      for ( int g=0 ; g< c->genp.nPar ; g++) {
	if ( c->genp.id[g] != 22 )
	  continue;
	if ( fabs( c->genp.momId[g] ) > 22 )
	  continue;
	if ( fabs( c->genp.status[g] ) != 1 )
	  continue;
	if ( fabs( c->genp.eta[g] ) > 1.44 ) 
	  continue;
	if ( c->genp.et[g] < 35 )
	  continue;
	genPhotonFlag = true;
      }
    }
    if ( !genPhotonFlag) 
      continue;
    
    c->GetEntry(jentry);
    evt.clear();
    evt.run   = c->evt.run;
    evt.evt = c->evt.evt;
    evt.hf4Pos = c->evt.hiHFplusEta4;
    evt.hf4Neg = c->evt.hiHFminusEta4;
    evt.hf4Sum = evt.hf4Pos + evt.hf4Neg;
    evt.cBin = -99;
    evt.pBin   = -99 ;
    if ((colli==kHIDATA)||(colli==kHIMC))   {
      evt.cBin = c->evt.hiBin;
      evt.pBin   = hEvtPlnBin->FindBin( c->evt.hiEvtPlanes[theEvtPlNumber] ) ;
    }
    else if ((colli==kPADATA)||(colli==kPAMC))   {
      evt.cBin =  getHfBin(evt.hf4Sum);
      if (  ((evt.cBin) < 0) || (evt.cBin) > 18 )
        cout << " Check the pA centrality..  cbin = " << evt.cBin << endl;
    }

    evt.trig = 0;
    evt.offlSel = (c->skim.pcollisionEventSelection > 0);
    evt.noiseFilt = (c->skim.pHBHENoiseFilter > 0);
    evt.anaEvtSel = c->selectEvent() && evt.trig;
    evt.vz = c->evt.vz;


    if ( ( (colli==kHIDATA)||(colli==kHIMC)||(colli==kPADATA)||(colli==kPAMC) ) && ( c->selectEvent() == 0 ))
      continue;
    if ( ( (colli==kPADATA)||(colli==kPPDATA) ) && ( c->skim.pVertexFilterCutGplus ==0 ) ) // No Pile up events
      continue;

    eSel++;      // OK.  This event is a collisional and no-noise event.

    int cBin = evt.cBin;
    int vzBin = hvz->FindBin(evt.vz)  ;
    hvz->Fill(evt.vz) ;
    if ( (vzBin<1) || ( vzBin > nVtxBin) )
      continue;

 
    // Reweight for vertex and centrality of MC 
    evt.vtxCentWeight = 1;
    double wVtx=1;
    double wCent=1;

//.........这里部分代码省略.........
  EWK2 ->Add(EWK1 );
//EWK3 ->Add(EWK2 );
//EWK  ->Add(EWK3 );
  EWK  ->Add(EWK2 );
  ttbar->Add(EWK  );
  Ztt  ->Add(ttbar);
  if(log){
#ifdef MSSM
    ggH  ->Add(bbH);
#else
    qqH  ->Add(VH );
    ggH  ->Add(qqH);
#endif
  }
  else{
#ifdef MSSM    
    bbH  ->Add(Ztt);
    ggH  ->Add(bbH);
#else
    VH   ->Add(Ztt);
    qqH  ->Add(VH );
    ggH  ->Add(qqH);
#endif
  }

  /*
    Mass plot before and after fit
  */
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);

  canv->cd();
  if(log){ canv->SetLogy(1); }
#if defined MSSM
  if(!log){ data->GetXaxis()->SetRange(0, data->FindBin(500)); } else{ data->GetXaxis()->SetRange(0, data->FindBin(1000)); };
#else
  data->GetXaxis()->SetRange(0, data->FindBin(350));
#endif
  data->SetNdivisions(505);
  data->SetMinimum(min);
  data->SetMaximum(max>0 ? max : std::max(maximum(data, log), maximum(Ztt, log)));
  data->Draw("e");

  TH1F* errorBand = (TH1F*)Ztt ->Clone();
  errorBand  ->SetMarkerSize(0);
  errorBand  ->SetFillColor(1);
  errorBand  ->SetFillStyle(3013);
  errorBand  ->SetLineWidth(1);
  for(int idx=0; idx<errorBand->GetNbinsX(); ++idx){
    if(errorBand->GetBinContent(idx)>0){
      std::cout << "Uncertainties on summed background samples: " << errorBand->GetBinError(idx)/errorBand->GetBinContent(idx) << std::endl;
      break;
    }
  }
  if(log){
    Ztt  ->Draw("histsame");
    ttbar->Draw("histsame");
    EWK  ->Draw("histsame");
    Fakes->Draw("histsame");
    $DRAW_ERROR
    ggH  ->Draw("histsame");
  }
  else{
    ggH  ->Draw("histsame");
    Ztt  ->Draw("histsame");
    ttbar->Draw("histsame");
    EWK  ->Draw("histsame");

void 
postfit(const char* inputfile, const char* analysis = "SM", const char* dataset = "2011+2012", const char* extra="", const char* extra2="", float min=0.1, float max=-1., bool log=true)
{
  // defining the common canvas, axes pad styles
  SetStyle(); gStyle->SetLineStyleString(11,"20 10");
  // switch for MSSM/SM
  bool MSSM = std::string(analysis) == std::string("MSSM");
  // determine label
  if (std::string(dataset) == std::string("2011"     )){ dataset = "CMS Preliminary,  H#rightarrow#tau#tau, 4.9 fb^{-1} at 7 TeV"; }
  if (std::string(dataset) == std::string("2012"     )){ 
	if (std::string(extra) == std::string("#mu#mu") ) dataset = "MS Preliminary,  H#rightarrow#tau#tau, 18.7 fb^{-1} at 8 TeV"; 
	else dataset = "MS Preliminary,  H#rightarrow#tau#tau, 19.4 fb^{-1} at 8 TeV";
  }
  if (std::string(dataset) == std::string("2011+2012")){ 
	if (std::string(extra) == std::string("#mu#mu") ) dataset = "CMS Preliminary,  H#rightarrow#tau#tau,  4.9 fb^{-1} at 7 TeV, 18.6 fb^{-1} at 8 TeV"; 
 	else dataset = "CMS Preliminary,  H#rightarrow#tau#tau,  4.9 fb^{-1} at 7 TeV, 19.4 fb^{-1} at 8 TeV";
	if (MSSM) dataset = "CMS Preliminary,  H#rightarrow#tau#tau,  4.9 fb^{-1} at 7 TeV, 12.1 fb^{-1} at 8 TeV";
  }
  // determine category tag
  const char* category_extra = "";
  if(std::string(extra2) == std::string("0jet_low"  )){ category_extra = "0 jet, low p_{T}";  }
  if(std::string(extra2) == std::string("0jet_high" )){ category_extra = "0 jet, high p_{T}"; }
  if(std::string(extra2) == std::string("0jet"      )){ category_extra = "0 jet";             }
  if(std::string(extra2) == std::string("boost_low" )){ category_extra = "1 jet, low p_{T}";  }
  if(std::string(extra2) == std::string("boost_high")){ category_extra = "1 jet, high p_{T}"; }
  if(std::string(extra2) == std::string("boost"     )){ category_extra = "1 jet";             }
  if(std::string(extra2) == std::string("vbf"       )){ category_extra = "2 jet (VBF)";       }
  if(std::string(extra2) == std::string("nobtag"    )){ category_extra = "No B-Tag";          }
  if(std::string(extra2) == std::string("btag"      )){ category_extra = "B-Tag";             }

  TFile* input = new TFile(inputfile);
  TH1F* Fakes  = refill((TH1F*)input->Get("Fakes"   ), "Fakes/QCD"); 
  TH1F* EWK    = refill((TH1F*)input->Get("EWK"     ), "EWK"      ); 
  TH1F* ttbar  = refill((TH1F*)input->Get("ttbar"   ), "ttbar"    ); 
  TH1F* Ztt    = refill((TH1F*)input->Get("Ztt"     ), "Ztt"      ); 
  TH1F* Zmm    = refill((TH1F*)input->Get("Zmm"     ), "Zmm"      ); 
  TH1F* ggH    = refill((TH1F*)input->Get("ggH"     ), "ggH"      ); 
  TH1F* data   = (TH1F*)input->Get("data_obs"); 
  // determine channel for etau Z->ee (EWK) will be shown separated from the rest (EWK1)
  TH1F* EWK1   = 0;
  if(std::string(extra) == std::string("e#tau_{h}")){
    EWK1 = refill((TH1F*)input->Get("EWK1"),  "EWK1");
  }
  TH1F* errorBand = (TH1F*)input->Get("errorBand");

  /* 
    mass plot before and after fit
  */
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);
  if(log) canv->SetLogy(1);
  // reduce the axis range if necessary for linea plots and SM
  if(MSSM && !log){ data->GetXaxis()->SetRange(0, data->FindBin(350)); } else{ data->GetXaxis()->SetRange(0, data->FindBin(490)); };
  if(!MSSM){ data->GetXaxis()->SetRange(0, data->FindBin(350)); }

  data->SetNdivisions(505);
  data->SetMinimum(min);
  if(Zmm){
    data->SetMaximum(max>0 ? max : std::max(maximum(data, log), maximum(EWK, log)));
  }
  else{
    data->SetMaximum(max>0 ? max : std::max(maximum(data, log), maximum(Ztt, log)));
  }
  data->Draw("e");

  if(log){
    if(Zmm){
      EWK  ->Draw("same");
      ttbar->Draw("same");
      Fakes->Draw("same");
      Zmm  ->Draw("same");
      Ztt  ->Draw("same");
    }
    else{
      Ztt  ->Draw("same");
      ttbar->Draw("same");
      EWK  ->Draw("same");
      if(EWK1){
	EWK1->Draw("same");
      }
      if(Fakes){ Fakes->Draw("same"); }
    }
    if(ggH) ggH  ->Draw("histsame");
  }
  else{
    if(ggH) ggH  ->Draw("histsame");
    if(Zmm){
      EWK->Draw("same");
      Fakes->Draw("same");
      ttbar->Draw("same");
      Zmm->Draw("same");
      Ztt->Draw("same");
    }
    else{
      Ztt  ->Draw("same");
      ttbar->Draw("same");
      EWK  ->Draw("same");
      if(EWK1){
	EWK1->Draw("same");
      }
      if(Fakes){ Fakes->Draw("same"); }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    }
    TH2D *LED[3][3][20];
    TH2D *PED[3][3];
    for(int iSpig = 0; iSpig < 3; iSpig++)for(int bb = 1; bb < 4; bb++){
        for(int ii=0; ii<Npoints; ii++){
           sprintf(spename,"spetest/spigot_%d/bb_%d/LED_HVset_%d_sp_%d_BB_%d",iSpig,bb,ii,iSpig,bb);
           LED[iSpig][bb-1][ii]=(TH2D *)tf->Get(spename);
        }
        sprintf(spename,"spetest/spigot_%d/bb_%d/PED_sp_%d_BB_%d",iSpig,bb,iSpig,bb);
        PED[iSpig][bb-1]=(TH2D *)tf->Get(spename);
    }


    for(int ii=0; ii<Npoints; ii++) {
      drawflag=false;
      int HV=0;
      for (int iSpig = 0; iSpig < 3; iSpig++) {
	for(int i = 0; i < 24; i++) {
            int bb=BB_MAP[i];
            int pmt=PMT_MAP[i];
            sprintf(spehistname,"led %d %d %d",ii,iSpig,i);
            TH1D *hspe_temp = (TH1D *)LED[iSpig][bb-1][ii]->ProjectionX(spehistname,pmt,pmt,"")->Clone();
            sprintf(spehistname,"ped %d %d %d",ii,iSpig,i);
            TH1D *hped = (TH1D *)PED[iSpig][bb-1]->ProjectionX(spehistname,pmt,pmt,"")->Clone();
            sscanf(&hspe_temp->GetTitle()[7],"%d",&HV);
	    hspe->Reset();
            sprintf (spehistname, "SumLED%d_sp_%d_ch_%d", HV, iSpig, i);
	    hspe->SetTitle(spehistname);

	    //combine bins of original SPE histogram
	    for(int ib=1; ib<=hspe_temp->GetNbinsX(); ib++) {
	      double bin_center = hspe_temp->GetBinCenter(ib);
	      if(bin_center>hspe->GetXaxis()->GetXmax()) continue;
	      int newbin = hspe->FindBin(bin_center);
	      double new_content = hspe->GetBinContent(newbin) + hspe_temp->GetBinContent(ib);
	      double new_error = sqrt(pow(hspe->GetBinError(newbin),2)+pow(hspe_temp->GetBinError(ib),2));
	      hspe->SetBinContent(newbin,new_content);
	      hspe->SetBinError(newbin,new_error);
	    }
	    TH1F* hspe_unscaled = (TH1F*)hspe->Clone("hspe_unscaled");

	    //renormalize bins of new SPE histogram
	    for(int ib=1; ib<=hspe->GetNbinsX(); ib++) {
	      double new_content = hspe->GetBinContent(ib)/hspe->GetXaxis()->GetBinWidth(ib)*hspe_temp->GetXaxis()->GetBinWidth(1);
	      double new_error = hspe->GetBinError(ib)/hspe->GetXaxis()->GetBinWidth(ib)*hspe_temp->GetXaxis()->GetBinWidth(1);
	      hspe->SetBinContent(ib,new_content);
	      hspe->SetBinError(ib,new_error);
	    }
	    
	    if(hspe_temp->Integral()==0) continue;
	    else drawflag=true;	  
    
            Nev = hspe_temp->Integral()*hspe_temp->GetXaxis()->GetBinWidth(1); 
	    
	    TF1 *fped = new TF1("fped","gaus",0, 80);
	    hped->Fit(fped,"NQR");
	    double pploc = fped->GetParameter(1), ppwidth = fped->GetParameter(2);
	    //cout<<"Ped only: ped mean "<<fped->GetParameter(1)<<", ped width "<<fped->GetParameter(2)<<" normalization "<<fped->GetParameter(0)<<endl;
	    hspe->Fit(fped, "NQ", "", pploc - 3*ppwidth, pploc + ppwidth);
	    //cout<<"SPE distribution: ped mean "<<fped->GetParameter(1)<<", ped width "<<fped->GetParameter(2)<<" normalization "<<fped->GetParameter(0)<<endl;
	    
	    
	    //estimate SPE peak location
	    int max_SPE_bin, maxbin, Nbins;
	    double max_SPE_height=0, minheight, max_SPE_location;
	    bool minflag = false;

//.........这里部分代码省略.........
    ptVeto  = (veto>=0) ? (*jets)[veto].Pt() : -99;
    etaVeto = (veto>=0) ? (*jets)[veto].Eta(): -99;
    phiVeto = (veto>=0) ? (*jets)[veto].Phi(): -99;
 
    isVetoInJets = 0;
    for(int l = 0 ; l < indexes.size() ; l++){
      if(lead>=0 && trail>=0 && (l!= lead && l!= trail) &&
	 (*jets)[indexes[l]].Pt()>PtVETO && ((*jets)[indexes[l]].Eta() - eta1)*((*jets)[indexes[l]].Eta() - eta2)<=0 )
	isVetoInJets = 1;
    }

    diTauSVFitMass = (*diTauSVfitP4)[0].M();
    diTauVisMass   = (*diTauVisP4)[0].M();
    diTauVisPt     = (*diTauVisP4)[0].Pt();
    diTauVisEta    = (*diTauVisP4)[0].Eta();
    diTauVisPhi    = (*diTauVisP4)[0].Phi();
    
    ptL1     = (*diTauLegsP4)[0].Pt();
    ptL2     = (*diTauLegsP4)[1].Pt();
    etaL1    = (*diTauLegsP4)[0].Eta();
    etaL2    = (*diTauLegsP4)[1].Eta();
    phiL1    = (*diTauLegsP4)[0].Phi();
    phiL2    = (*diTauLegsP4)[1].Phi();
    dPhiL1L2 =  abs((*diTauLegsP4)[0].Phi()-(*diTauLegsP4)[1].Phi()) > TMath::Pi() ? 
      -abs( (*diTauLegsP4)[0].Phi()-(*diTauLegsP4)[1].Phi() ) + 2*TMath::Pi()  :
      abs( (*diTauLegsP4)[0].Phi()-(*diTauLegsP4)[1].Phi() ) ;

    diTauCharge_    = diTauCharge;
    
    

    ////////////////////////////////////////////////////////////////////

    float scalarSumPt     = ( *diTauLegsP4)[0].Pt() + (*METP4)[0].Pt();
    float vectorSumPt     = ((*diTauLegsP4)[0] +      (*METP4)[0]).Pt() ;

    MtLeg1_     = TMath::Sqrt( scalarSumPt*scalarSumPt - vectorSumPt*vectorSumPt ) ;
    MEt     = (*METP4)[0].Et();
    combRelIsoLeg1DBeta    = (chIsoLeg1+ std::max( nhIsoLeg1+phIsoLeg1-0.5*(nhIsoPULeg1),double(0.0)))/(*diTauLegsP4)[0].Pt();

    tightestHPSDBWP_ = tightestHPSDBWP;
    visibleTauMass_  = visibleTauMass;
    decayMode_       = decayMode;
    numPV_           = numPV;
    sampleWeight     = scaleFactor; 
    puWeight         = (std::string(sample.Data())).find("Data")!=string::npos ? 1.0 : pileupWeight(nPUVertices);   

    HqTWeight = histo!=0 ? histo->GetBinContent( histo->FindBin( (*genVP4)[0].Pt() ) ) : 1.0;

    if((std::string(sample.Data())).find("Data")!=string::npos){
      
      HLTx =  float((*triggerBits)[0]);  //HLT_IsoMu15_LooseIsoPFTau15_v8
      bool isTriggMatched = (*tauXTriggers)[0] && (*tauXTriggers)[1] ; //hltSingleMuIsoL3IsoFiltered15 && hltOverlapFilterIsoMu15IsoPFTau15
      HLTmatch = isTriggMatched ? 1.0 : 0.0;
      HLTweightTau = 1.0;
      HLTweightMu  = 1.0;
      SFTau        = 1.0;
      SFMu         = 1.0;
      
    } 
    else{

      HLTx  =  float((*triggerBits)[0]); //HLT_IsoMu15_LooseIsoPFTau15_v9
      bool isTriggMatched = (*tauXTriggers)[0] && (*tauXTriggers)[1] ; //hltSingleMuIsoL3IsoFiltered15 && hltOverlapFilterIsoMu15IsoPFTau15
      HLTmatch = isTriggMatched ? 1.0 : 0.0;

      HLTweightTau  = TMath::Abs((*diTauLegsP4)[1].Eta())<1.5 ?  
	ratioTauHLTBL->Eval( (*diTauLegsP4)[1].Pt() ) : ratioTauHLTEC->Eval( (*diTauLegsP4)[1].Pt() );
      HLTweightMu  = TMath::Abs((*diTauLegsP4)[0].Eta())<1.5 ?  
	ratioMuHLTBL->Eval( (*diTauLegsP4)[0].Pt() ) : ratioMuHLTEC->Eval( (*diTauLegsP4)[0].Pt() );
      SFTau  = 1.0;
      SFMu   =  TMath::Abs((*diTauLegsP4)[0].Eta())<1.5 ?
	ratioMuIDIsoBL->Eval( (*diTauLegsP4)[0].Pt() ): 
	ratioMuIDIsoEC->Eval( (*diTauLegsP4)[0].Pt() );
      
    }
   
    muFlag_          = muFlag;
    genDecay_        = genDecay ;
    event_           = event;
    run_             = run;
    lumi_            = lumi;
    
    outTreePtOrd->Fill();
  }


  file->Close();
  
  if(SAVE) outFile->Write();
  outFile->Close();

  delete jets; delete diTauLegsP4; delete diTauVisP4; delete diTauSVfitP4; delete genDiTauLegsP4;
  delete tauXTriggers; delete triggerBits;
  delete METP4; delete genVP4; delete genMETP4;
  delete HqT;
  
  return;

}

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

    for( unsigned int ibin = 1 ; ibin <= nx ; ibin++ ){
      for( unsigned int jbin = 1 ; jbin <= ny ; jbin++ ){

	float mg = heff[i]->GetXaxis()->GetBinCenter(ibin);
	float ml = heff[i]->GetYaxis()->GetBinCenter(jbin);

	float eff    = heff[i]->GetBinContent(ibin,jbin);
	// float effup  = heffup[i]->GetBinContent(ibin,jbin);
	// float effdn  = heffdn[i]->GetBinContent(ibin,jbin);

	if( eff   < 1e-20 ) continue;

	// float dup    = effup/eff-1;
	// float ddn    = 1-effdn/eff;
	// float djes   = 0.5 * (dup+ddn);
	// hjes[i]->SetBinContent(ibin,jbin,djes);

	// float toterr = sqrt( 0.04*0.04 + 0.05*0.05 + 0.05*0.05 + djes*djes );

	//float this_ul = getObservedLimit( cuts.at(i) , toterr );
	float this_ul = getObservedUpperLimit( cuts.at(i) );
	float xsecul  = this_ul / ( lumi * eff );

	//float this_ul_exp = getExpectedLimit( cuts.at(i) , toterr );
	float this_ul_exp = getExpectedUpperLimit( cuts.at(i) );
	float xsecul_exp  = this_ul_exp / ( lumi * eff );

	if( eff > 0 ){
	  hxsec[i]->SetBinContent(ibin,jbin, xsecul );
	  hxsec_exp[i]->SetBinContent(ibin,jbin, xsecul_exp );
	}

	int   bin = refxsec->FindBin(mg);
	float xsec = refxsec->GetBinContent(bin);

	hexcl[i]->SetBinContent(ibin,jbin,0);
	if( xsec > xsecul )   hexcl[i]->SetBinContent(ibin,jbin,1);
	//cout << "ibin jbin mg xsec " << ibin << " " << jbin << " " << mg << " " << xsec << endl;
      }
    }
  }

  delete ctemp;

  cout << endl << endl;

  //--------------------------------------------------
  // make pretty pictures
  //--------------------------------------------------
  
  TLatex *t = new TLatex();
  t->SetNDC();
  t->SetTextSize(0.04);

  TCanvas* can[nsig];

  for( unsigned int i = 0 ; i < nsig ; ++i ){
  
    //can[i] = new TCanvas(Form("can_%i",i),Form("can_%i",i),1200,600);
    //can[i]->Divide(2,1);
    //can[i] = new TCanvas(Form("can_%i",i),Form("can_%i",i),1800,600);
    //can[i]->Divide(3,1);
    can[i] = new TCanvas(Form("can_%i",i),Form("can_%i",i),1000,1000);
    can[i]->Divide(2,2);

//.........这里部分代码省略.........
  }
  if (flags & 1<<5) userScale *= ((TParameter<float> *)input->Get("trgEff"))->GetVal();
  if (flags & 1<<4) userScale *= ((TParameter<float> *)input->Get("trgDataMcScaleFactor"))->GetVal();
  if (flags & 1<<1) userScale *= ((TParameter<float> *)input->Get("muScaleFactor"))->GetVal();
  float eleScaleFactorEB = ((TParameter<float> *)input->Get("eleScaleFactorEB"))->GetVal();
  float eleScaleFactorEE = ((TParameter<float> *)input->Get("eleScaleFactorEE"))->GetVal();
  float lumiScaleFactorEB = ((TParameter<float> *)input->Get("lumiScaleFactorEB"))->GetVal();
  float lumiScaleFactorEE = ((TParameter<float> *)input->Get("lumiScaleFactorEE"))->GetVal();

  //cout << "Scalefactor = " << userScale << endl;

  // get the tree
  TTree *tree;
  tree = (TTree *)input->Get(treeName);

  // get branches
  float var;
  bool passTrg;
  bool passHeep;
  float puWeight = 1.;
  int eCharge;
  int muCharge;
  int evtRegion;
  float cutVar = 0.;
  float fakeRate = 0.;
  tree->SetBranchStatus("*",0); //disable all branches
  tree->SetBranchStatus(brName,1);
  tree->SetBranchAddress(brName, &var);
  if (flags & 1<<8) {
    tree->SetBranchStatus("passTrg",1);
    tree->SetBranchAddress("passTrg", &passTrg);
  }
  if (signs != 0) {
    tree->SetBranchStatus("eCharge",1);
    tree->SetBranchStatus("muCharge",1);
    tree->SetBranchAddress("eCharge", &eCharge);
    tree->SetBranchAddress("muCharge", &muCharge);
  }
  if (region < 2) {
    tree->SetBranchStatus("evtRegion",1);
    tree->SetBranchAddress("evtRegion", &evtRegion);
  }
  if (flags & 1) {
    tree->SetBranchStatus("puWeight",1);
    tree->SetBranchAddress("puWeight", &puWeight);
  }
  if (cutVariable[0] != '\0') {
    tree->SetBranchStatus(cutVariable,1);
    tree->SetBranchAddress(cutVariable, &cutVar);
  }
  if (flags & 1<<9) {
    tree->SetBranchStatus("passHeep",1);
    tree->SetBranchStatus("fakeRate",1);
    tree->SetBranchAddress("passHeep", &passHeep);
    tree->SetBranchAddress("fakeRate", &fakeRate);
  }

  Long64_t nEntries = (*tree).GetEntries();
  for (unsigned int i = 0; i < nEntries; ++i) {
    tree->GetEntry(i);

    // trigger fired?
    if ((flags & 1<<8) && passTrg == false) continue;

    // select electron region
    if (evtRegion == 0 && region == 1) continue;
    if (evtRegion == 1 && region == 0) continue;

    float scaleFactor = userScale;
    // set lumi and electron scalefactor according to detector region
    if (evtRegion == 0 && flags & 1<<2) scaleFactor *= eleScaleFactorEB;
    if (evtRegion == 1 && flags & 1<<2) scaleFactor *= eleScaleFactorEE;
    if (evtRegion == 0 && flags & 1<<3) scaleFactor *= lumiScaleFactorEB;
    if (evtRegion == 1 && flags & 1<<3) scaleFactor *= lumiScaleFactorEE;

    // PU reweight
    if (flags & 1) scaleFactor *= puWeight;

    // get only the desired charge combination. Scheme emu -3 to +3: -+, +-, OS, ALL, SS, ++, --
    if (signs < 0 && (eCharge * muCharge) > 0) continue; // OS
    if (signs > 0 && (eCharge * muCharge) < 0) continue; // SS
    if (abs(signs) == 3 && eCharge > 0) continue; // e-mu+ or e-mu-
    if (abs(signs) == 2 && eCharge < 0) continue; // e+mu- or e+mu+

    // user defined cut
    if (cutVariable[0] != '\0')
      if (cutVar < cutLow || cutVar >= cutHigh) continue;
    
    if (flags & 1<<9) {
      if (!passHeep) scaleFactor *= fakeRate / (1 - fakeRate);
      else continue;
    }

    if (normToBinWidth) scaleFactor /= histo->GetBinWidth(histo->FindBin(var));
    histo->Fill(var, scaleFactor);
  }

  //cout << "integral: " << histo->Integral() << "       overflow: " << histo->GetBinContent(histo->GetNbinsX() + 1) << endl;
  return histo;
}

void fitPtOverMCJLST(int mass = 125, int LHCsqrts = 7, int whichtype = 1, 
		     bool correctErrors = false, /* string changeParName = "", */
		     bool showErrorPDFs = false, string systString = "Default")

// whichtype
// 0 - gg Signal
// 1 - VBF Signal
// 2 - ZZ
// 3 - ZX
// 4 - ggZZ
// 5 - WH
// 6 - ZH
// 7 - ttH

{

  string changeParName = "";
  if (systString == "Default") changeParName = "up";

  string nameSample[8] = {"gg","vbf","zz","zx","ggzz","wh","zh","tth"};
  float maxType[8] = {2.4,3.2,1.6,1.6,1.6,3.2,3.2,3.2};   
  float rebinType[8] = {1,2,1,1,4,10,10,40};
  
  for (int t = 0; t < 8; t++) {
    if (mass > 150) maxType[t] = int(117.90*maxType[t]/sqrt(mass-10.91))/10.;
  }

  char fileToOpen[200];
  sprintf(fileToOpen,"selRootFiles/PToverM_%s%d_SEL_%dTeV.root",nameSample[whichtype].c_str(),mass,LHCsqrts);
  // if (whichtype == 3) sprintf(fileToOpen,"PTOVERM_%s_SEL_allTeV.root",nameSample[whichtype].c_str());

  RooRealVar* ptoverm = new RooRealVar("ptoverm","p_{T}/M^{4l}",0.,maxType[whichtype],"GeV/c");
 
  TFile input(fileToOpen);
  // if (systString == "Mass" || systString == "Mela") {
  //  sprintf(fileToOpen,"ptovermH_%sUp",systString.c_str());
  // } else {
  sprintf(fileToOpen,"ptovermH_%s",systString.c_str());
  //}
  TH1F* ptovermH = (TH1F*)input.Get(fileToOpen);
  
  if (rebinType[whichtype] > 1) ptovermH->Rebin(rebinType[whichtype]);
  if (maxType[whichtype] < ptovermH->GetBinLowEdge(ptovermH->GetNbinsX() + 1) - ptovermH->GetBinWidth(1)) {
    int theBin = ptovermH->FindBin(maxType[whichtype]);
    ptovermH->GetXaxis()->SetRange(1,theBin-1);
  }

  gROOT->ProcessLine(".L mytdrstyle.C");
  gROOT->ProcessLine("setTDRStyle()");
  
  // cout << endl << "Signal " << endl;   
  ptoverm->setBins(ptovermH->GetNbinsX());

  RooDataHist* rdh = new RooDataHist("rdh","Some dataset",RooArgList(*ptoverm),Import(*ptovermH,kFALSE));
 
  // fit definitions
  // RooWorkspace *ws = new RooWorkspace("ws");

  RooRealVar m("m","emme", 1.,0.01, 40.);
  RooRealVar n("n","enne", 0.93, 0.05, 15.);
  RooRealVar n2("n2","enne2", 0.75, 0.5, 15.);
  RooRealVar bb("bb","bibi",0.02, 0.000005, 20.0);
  RooRealVar T("T","tti",0.2,0.00000005,1.);
  RooRealVar bb2("bb2","bibi2",0.02, 0.0005, 10.0);
  RooRealVar fexp("fexp","f_exp",0.02, 0.0, 1.0);
  if (whichtype == 0) {
    if (LHCsqrts == 8) {
      m.setVal(3.319);   // m.setConstant(kTRUE);
      n.setVal(0.7606);    if (systString != "Default" || mass != 125) n.setConstant(kTRUE); 
      n2.setVal(0.8061);   n2.setConstant(kTRUE);
      bb.setVal(3.728);   // bb.setConstant(kTRUE);
      T.setVal(0.00333);   // T.setConstant(kTRUE);
      bb2.setVal(1.7172);    // bb2.setConstant(kTRUE);
      fexp.setVal(0.002144);   if (systString != "Default" || mass != 125) fexp.setConstant(kTRUE);
    } else {
      m.setVal(0.061);   // m.setConstant(kTRUE);
      n.setVal(1.6141);   if (systString == "Resummation" || systString == "TopMass")  n.setConstant(kTRUE);
      n2.setVal(1.3294);   n2.setConstant(kTRUE);
      bb.setVal(4.2761);   // bb.setConstant(kTRUE);
      T.setVal(0.0361);   // T.setConstant(kTRUE);
      bb2.setVal(1.6643);   bb2.setConstant(kTRUE);
      fexp.setVal(0.0004);   // fexp.setConstant(kTRUE);
    }
  } else if (whichtype == 1) {
    m.setVal(1.006);   // m.setConstant(kTRUE);
    n.setVal(10.939);   n.setConstant(kTRUE);
    n2.setVal(1.1448);   n2.setConstant(kTRUE);
    bb.setVal(0.02048);   bb.setConstant(kTRUE);
    T.setVal(0.16115);   if (systString.find("Mela") != string::npos) T.setConstant(kTRUE); // T.setConstant(kTRUE);
    bb2.setVal(1.0024);   bb2.setConstant(kTRUE);
    fexp.setVal(0.005);   fexp.setConstant(kTRUE);
    if (mass > 300) {
      fexp.setVal(0.0);   fexp.setConstant(kFALSE);
    }
    if (mass > 500) {
      bb2.setVal(5.0);  //  bb2.setConstant(kFALSE);
    }
    if (mass > 500) {
      bb.setVal(15.0);  //  bb.setConstant(kFALSE);
    }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
	hE->Fill(myBSDLatePE/CALSum);
      
      // Fill Efrac for pions
      if( ptype == 2 && (penergy >= 250 || penergy <= 350) )
	hPi->Fill(myBSDLatePE/CALSum);
      // Increment counts, since everything has to be equally weighted
      if( ptype == 2 && penergy == 250 ) NPi250++;
      if( ptype == 2 && penergy == 300 ) NPi300++;
      if( ptype == 2 && penergy == 350 ) NPi350++;
      
      if( ptype == 1 && penergy == 75 ) NE75++;
      if( ptype == 1 && penergy == 100 ) NE100++;
      if( ptype == 1 && penergy == 125 ) NE125++;
      if( ptype == 1 && penergy == 150 ) NE150++;
      if( ptype == 1 && penergy == 175 ) NE175++;
    }
  
  hE->SetLineColor(kRed);
  hE->SetTitle("Efrac for 250-350 GeV Pions (Blue) and 75-175 GeV Electrons (Red) and CALSum Between 6000 and 9000");
  hE->GetXaxis()->SetTitle("Efrac (\"detector\" units)");
  hE->GetXaxis()->CenterTitle();
  hE->Scale(1./hE->GetEntries());
  hE->Draw();
  hPi->Scale(1./hPi->GetEntries());
  hPi->Draw("sames");

  Double_t reg_power[1000];
  Double_t e_accept[1000];
  Double_t pi_accept[1000];
  Double_t endbin[1000];
  Int_t Npts = 0;
  Int_t Nskip = 0;

  Int_t start = hE->FindBin(0);
  Int_t nbins = 40;
  // ToDo: Add error bars
  cout << "Start = " << start << endl;

  Double_t x_err[1000];
  Double_t y_err[1000];

  Double_t num = 0;
  Double_t den = 0;
  Double_t num_err = 0;
  Double_t den_err = 0;
  Double_t Ne = 0;
  Double_t Ne_total = 0;
  Double_t Npi = 0;
  Double_t Npi_total = 0;
  Double_t Ne_entries = hE->GetEntries();
  Double_t Npi_entries = hPi->GetEntries();
  for( int i = 0; i < nbins; i++ )
    {
      cout << "Efrac = " << hE->GetBinLowEdge(start+i) << endl;
      if( hPi->Integral(start,start+i) > 10./Ne_entries && 
	  hE->Integral(start,start+i) > 10./Npi_entries )
	{
	  num_err = 0;
	  den_err = 0;
	  Ne = 0;
	  Npi = 0;
	  num = 0;
	  den = 0;
	  Npi_total = 0;
	  Ne_total = 0;
	  

//.........这里部分代码省略.........
	xmin=-3.8;
	xmax= 3.8;
	binwidth=0.01;
      }
    }
    // fill data in binned TH1F
    hist= new TH1F ( variable, variable, (int)((xmax-xmin)/binwidth), xmin, xmax);
    TH1F* ori=(TH1F*)hist->Clone("original points");
    std::cout << "fine binned theory prediction has " << hist->GetNbinsX() << " bins" << std::endl;
    std::cout << "loaded values from .dat file: " << std::endl;
    // list all data values loaded and the corresponding bin
    for(int bin=1; bin<=Nbins; ++bin){
      double x=-999;
      double y=-999;
      // NB: choose if loaded data is interpreted as points with nothing 
      //     between (like kidonakis) or as integrated over the bin range (like ahrens)
      if(points){
	// check if you are still inside the array
	//std::cout << "data point " << bin-1 << "/" << sizeof(Xvalues)/sizeof(double) << std::endl;
	if(rawHist->GetPoint(bin-1, x, y)!=-1){
	  //x=Xvalues[bin-1]; // get value from data points
	  x+=0.5*binwidth;  // add half the binwidth to get the center of the bin
	}
      }
      else{
	x=hist->GetBinCenter(bin);  // get bin center
	y=rawHist->GetY()[bin-1];
      }
      if(x!=-999){
	std::cout << "data point: " << bin;
	std::cout << "(<x>=" << x << ")-> bin";
	// get bin in target (fine binned) plot
	int bin2=bin;
	if(points) bin2=hist->FindBin(x);
	//double y=Yvalues[bin2-1];
	std::cout << bin2 << " (";
	std::cout << hist->GetBinLowEdge(bin2) << ".." << hist->GetBinLowEdge(bin2+1);
	std::cout << "): " << y << std::endl;
	// fill target (fine binned) plot
	if(!points) hist->SetBinContent(bin2, y);
	// mark bins coming from the original prediction
	ori->SetBinContent(bin2, 1.);
	// -------------------------------
	// fit range without data entries
	// -------------------------------
	// NB: needed if loaded data is interpreted as points with nothing 
	//     between (like kidonakis)
	if(points){
	  // perform a linear fit wrt previous point
	  // get the two points (this bin and the previous one)
	  int binPrev= (bin==1&&variable=="topPt") ? 0 : hist->FindBin(Xvalues[bin-2]+0.5*binwidth);
	  double x2=-1;
	  double x1= 0;
	  double y2=-1;
	  double y1= 0.;
	  rawHist->GetPoint(bin-1, x2, y2);
	  x2+=0.5*binwidth;
	  if(bin==1&&variable=="topPt"){
	    y1=0;
	    x1=0;
	  }
	  else{
	    rawHist->GetPoint(bin-2, x1, y1);
	    x1+=0.5*binwidth;
	  }
	  // calculate linear funtion

void GenerateInputFile_ge4b( double mgl=-1., double mlsp=-1., double target_susy_all0lep=-1. ) {

  TChain* dyTree = new TChain("treeZ") ;
  int nAdded = dyTree->Add("files15fb_8TeV_old1/DY.root") ;
  if ( nAdded <= 0 ) {
     printf("\n\n\n *** No treeZ in files15fb_8TeV_old1/DY.root\n\n\n") ;
     return ;
  }

  double t1bbbbWeight(0.) ;
  TChain chainT1bbbb("tree") ;
  char susycutstring[1000] ;
  sprintf( susycutstring, "&&mgluino==%.0f&&mlsp==%.0f", mgl, mlsp ) ;
  TString susycut( susycutstring ) ;
  if ( mgl>0. && mlsp>0. ) {
     nAdded = chainT1bbbb.Add("files5fb_MT/T1bbbb.root") ;
     if ( nAdded <= 0 ) {
        printf("\n\n\n *** No tree in files5fb_MT/T1bbbb.root\n\n\n") ;
        return ;
     }
     TFile f("referenceXSecs.root") ;
     TH1F* xsechist = (TH1F*) f.Get("gluino8TeV_NLONLL") ;
     if ( xsechist==0x0 ) { printf("\n\n *** can't find reference Xsec histogram in referenceXSecs.root.\n\n") ; return ; }
     int theBin = xsechist->FindBin( mgl ) ;
     if ( theBin <=0 || theBin > xsechist->GetNbinsX() ) {
        printf("\n\n *** can't find bin for mgl=%g.  Returned %d\n\n", mgl, theBin ) ;
        return ;
     }
     double xsec = xsechist->GetBinContent( theBin ) ;
     printf("\n\n T1bbbb xsec for mgl=%g is %g\n\n", mgl, xsec ) ;
     t1bbbbWeight = 1.5*xsec ;  //in pb. scan has 10k events, so nScan*1.5*xsec = events in 15fb-1
     //////  t1bbbbWeight = 0.1*xsec ;  //in pb. T1tttt scan has 50k events, so nScan*0.1*xsec = events in 5fb-1
     printf("\n\n Susy ttree cut: %s\n\n", susycutstring ) ;
  }


  TChain chainQCD("tree") ;
   //--- these have high weight
//chainQCD.Add("files15fb_8TeV/QCD-50to80.root");
//chainQCD.Add("files15fb_8TeV/QCD-80to120.root");
  chainQCD.Add("files15fb_8TeV/QCD-120to170.root");
  chainQCD.Add("files15fb_8TeV/QCD-170to300.root");
   //--- below here, these have weight less than one.
  chainQCD.Add("files15fb_8TeV/QCD-300to470.root");
  chainQCD.Add("files15fb_8TeV/QCD-470to600.root");
  chainQCD.Add("files15fb_8TeV/QCD-600to800.root");
  chainQCD.Add("files15fb_8TeV/QCD-800to1000.root");
  chainQCD.Add("files15fb_8TeV/QCD-1000to1400.root");
  chainQCD.Add("files15fb_8TeV/QCD-1400to1800.root");
  chainQCD.Add("files15fb_8TeV/QCD-1800.root");

  TChain chainTT("tree") ;
  chainTT.Add("files15fb_8TeV/TT.root") ;

  TChain chainZnn("tree") ;
  chainZnn.Add("files15fb_8TeV/Zinv-100to200.root") ;
  chainZnn.Add("files15fb_8TeV/Zinv-200to400.root") ;
  chainZnn.Add("files15fb_8TeV/Zinv-400.root") ;

  TChain chainWJets("tree") ;
  chainWJets.Add("files15fb_8TeV/WJets-250to300.root") ;
  chainWJets.Add("files15fb_8TeV/WJets-300to400.root") ;
  chainWJets.Add("files15fb_8TeV/WJets-400.root") ;
  chainWJets.Add("files15fb_8TeV/T-s.root") ;
  chainWJets.Add("files15fb_8TeV/T-t.root") ;
  chainWJets.Add("files15fb_8TeV/T-tW.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-s.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-t.root") ;
  chainWJets.Add("files15fb_8TeV/Tbar-tW.root") ;



      char qcdinputfile[9][1000] = {
        "files15fb_8TeV/QCD-120to170.root"
       ,"files15fb_8TeV/QCD-170to300.root"
       ,"files15fb_8TeV/QCD-300to470.root"
       ,"files15fb_8TeV/QCD-470to600.root"
       ,"files15fb_8TeV/QCD-600to800.root"
       ,"files15fb_8TeV/QCD-800to1000.root"
       ,"files15fb_8TeV/QCD-1000to1400.root"
       ,"files15fb_8TeV/QCD-1400to1800.root"
       ,"files15fb_8TeV/QCD-1800.root"
      } ;

      char qcdsamplename[9][100] = {
        "qcd_0120_to_0170"
       ,"qcd_0170_to_0300"
       ,"qcd_0300_to_0470"
       ,"qcd_0470_to_0600"
       ,"qcd_0600_to_0800"
       ,"qcd_0800_to_1000"
       ,"qcd_1000_to_1400"
       ,"qcd_1400_to_1800"
       ,"qcd_1800_to_9999"
      } ;



  gROOT->Reset();

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