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

void offline(const char* FileName="test")
{
    TH1F::SetDefaultSumw2();
    // Set Style parameters for this macro
    gStyle->SetOptTitle(1); // Show Title (off by default for cleanliness)
    gErrorIgnoreLevel = kError; // Set Verbosity Level (kPrint shows all)

    // Set Output options
    Int_t number;
    Bool_t fPaintAll = checkPaintAllTrigs();
    checkBatchMode();
    Bool_t makePDF = checkMakePDF();
    Bool_t makeROOT= checkMakeRoot();

    // Open ROOT File
    char name[1000];
    sprintf(name,"/Users/zach/Research/rootFiles/run12NPEhPhi/%s.root",FileName);
    TFile *f = new TFile(name,"READ");
    if (f->IsOpen()==kFALSE)
    {   std::cout << "!!! File Not Found !!!" << std::endl;
        exit(1);
    }
    // f->ls(); // - DEBUG by printing all objects in ROOT file

    char fname[100];
    TFile* file;
    if(makeROOT) {
        sprintf(fname,"/Users/zach/Research/rootFiles/run12NPEhPhi/%s_processed.root",FileName);
        file = new TFile(fname,"RECREATE");
        if (file->IsOpen()==kFALSE)
        {
            std::cout << "!!! Outfile Not Opened !!!" << std::endl;
            makeROOT = kFALSE;
        }
    }

    const Int_t numPtBins = anaConst::nPtBins;
    Float_t lowpt[numPtBins],highpt[numPtBins];
    for(Int_t c=0; c< numPtBins; c++) {
        lowpt[c] = anaConst::lpt[c];
        highpt[c] = anaConst::hpt[c];
    }
    Float_t hptCut=anaConst::hptCut;
    const Int_t numTrigs = 4;
    Double_t epsilon[numPtBins] = {0.593164, 0.626663, 0.655916, 0.674654, 0.685596, 0.700600, 0.716682, 0.724638, 0.713977, 0.730550, 0.735204, 0.744336, 0.761323, 0.758423};
    Float_t hptMax=25; // Set max above range to allow overflow
    Float_t lowPhi=anaConst::lowPhi, highPhi=anaConst::highPhi;
    Double_t pu[2][numPtBins][numTrigs]; // To store fit parameters for later use
    Double_t hhNorm, HHScale, hadPur;

    TH1D * LSIM[numPtBins][numTrigs];
    TH1D * USIM[numPtBins][numTrigs];
    TH1D * USIMNP[numPtBins][numTrigs];
    TH1D * LSIMNP[numPtBins][numTrigs];
    TH1D * INCLNP[numPtBins][numTrigs];
    TH1D * INCL[numPtBins][numTrigs];
    TH1D * LSIM2[numPtBins][numTrigs];
    TH1D * USIM2[numPtBins][numTrigs];
    TH1D * INCL2[numPtBins][numTrigs];
    TH1D * LSMM[numPtBins][numTrigs];
    TH1D * USMM[numPtBins][numTrigs];
    TH1D * HHDP[numPtBins][numTrigs];
    TH1D * NSPI[numPtBins][numTrigs];
    TH2F * MON[numTrigs];
    TH1F * AVGTRKS[numTrigs];
    TH2F * MON2[numTrigs];
    TH1F * AVGTRKS2[numTrigs];
    TH2F* mh2PhiQPt[numTrigs];
    TH2F* mh2nSigmaEPt[numTrigs];
    TH2F* mh2nSigmaEPt_eID[numTrigs];
    TH3F* mh3DelPhiIncl[numTrigs];
    TH3F* mh3DelPhiPhotLS[numTrigs];
    TH3F* mh3DelPhiPhotUS[numTrigs];
    TH3F* mh3DelPhiPhotUSNP[numTrigs];
    TH3F* mh3DelPhiPhotLSNP[numTrigs];
    TH3F* mh3DelPhiInclNP[numTrigs];
    TH3F* mh3DelPhiInclWt[numTrigs];
    TH3F* mh3DelPhiPhotLSWt[numTrigs];
    TH3F* mh3DelPhiPhotUSWt[numTrigs];
    TH2F* mh2InvMassPtLS[numTrigs];
    TH2F* mh2InvMassPtUS[numTrigs];
    TH2F* mh2PtEZdcx[numTrigs];
    TH2F* mh2PtEZdcxUS[numTrigs];
    TH2F* mh2PtEZdcxLS[numTrigs];
    TH2F* mh2PtEZdcxHad[numTrigs];
    TH3F* mh3nTracksZdcx[numTrigs];
    TH3F* mh3nTracksZdcxUS[numTrigs];
    TH3F* mh3nTracksZdcxLS[numTrigs];
    TH3F* mh3nTracksZdcxHad[numTrigs];
    TH3F* mh3MixedDelPhi;
    TH3F* mh3MixedDelEta;
    TH3F* mh3MixedEtaPhi;
    TH3F* mh3DelPhiHadHad[numTrigs];
    TH2F* mh2nSigmaPion[numTrigs];
    TH1F* mh1PtHadTracks[numTrigs];
    TH1F* mh1PtETracks[numTrigs];
    TH1D* projHPhi[numPtBins][numTrigs];
    TH1D* projnSigmaE[numPtBins][numTrigs];
    TH1D* projnSigmaE_eID[numPtBins][numTrigs];
    TH1D* projDelPhiIncl[numPtBins][numTrigs];
//.........这里部分代码省略.........

void AllInOne_combineLepton(const std::string& dataset,int btag=0, double timesX=1, bool useNewShape=false){

  std::string useNewShapeText = (useNewShape) ? "_newshape" : "";

  char wsFileName[700];
  sprintf(wsFileName,"datacards/400/hzz2l2q_ee%db.input.root",btag);
  //sprintf(wsFileName,"PROVA/datacards_%s%s/hzz2l2q_ee%db.input.root",dataset.c_str(), useNewShapeText.c_str(),btag);

  gSystem->Load("libRooFit");
  gSystem->Load("libFFTW");

  string histoName[3];
  histoName[0]="mZZ_kinfit_hiMass_0btag";
  histoName[1]="mZZ_kinfit_hiMass_1btag";
  histoName[2]="mZZ_kinfit_hiMass_2btag";
  string btagName[3]={"0b","1b","2b"};

  double LumiScale=0.;
  if( dataset=="Run2011A_FULL" ) LumiScale = 2100.;
  else if( dataset=="LP11" ) LumiScale = 1600.;
  else {
    std::cout << "Unknown dataset '" << dataset << "'. Exiting." << std::endl;
  }

  RooDataSet *data_bkg;  
  RooDataSet *data_temp;
  TFile *file;

  string cutString[3];
  cutString[0]="nBTags==0 && (mZjj>75 && mZjj<105) && mZZ>183";
  cutString[1]="nBTags==1 && (mZjj>75 && mZjj<105) && mZZ>183";
  cutString[2]="nBTags==2 && (mZjj>75 && mZjj<105) && mZZ>183";

  int binWidth=20;
  int highBin=750;
  int lowBin=150;

//double muonEff[3]={.576,.548,.489};

//double expSig[3];
//expSig[0]=5.65*LumiScale ;
//expSig[1]=4.89*LumiScale ;
//expSig[2]=2.37*LumiScale ;

  char alphaFileName[200];
  sprintf( alphaFileName, "alphaFile_%s_%dbtag_ALL.root", dataset.c_str(), btag);
  TFile* alphaFile = TFile::Open(alphaFileName);
  TTree* treeSidebandsDATA_alphaCorr = (TTree*)alphaFile->Get("sidebandsDATA_alpha");
  TH1D* h1_mZZ_sidebands_alpha = new TH1D("mZZ_sidebands_alpha", "", 65, 150., 800.);
  char sidebandsCut_alpha[500];
  sprintf(sidebandsCut_alpha, "eventWeight_alpha*(isSidebands && nBTags==%d)", btag);
  treeSidebandsDATA_alphaCorr->Project("mZZ_sidebands_alpha", "mZZ", sidebandsCut_alpha);
  float expBkg = h1_mZZ_sidebands_alpha->Integral();
  std::cout <<  "++++ expBkg: " << expBkg << std::endl;



  stringstream convert;

  // --------------------- measurable (ZZ invariant mass) ----------------
  string temp;
  temp="m_{ZZ}";
  RooRealVar mZZ("mZZ",temp.c_str(),lowBin,highBin);
  RooRealVar nBTags("nBTags","nBTags",-1.,3.);
  RooRealVar leptType("leptType","leptType",0,1);
  RooRealVar mZjj("mZjj","mZjj",0,200.);

  // ----------------- get parameters from data cards! -----------------
  
  TFile *wsFile = new TFile(wsFileName);
  RooWorkspace *ws = (RooWorkspace*) wsFile->Get("w");

  // ==================== defining bkg PDF ==========================
  // ------------------------ fermi ------------------------------
  RooRealVar cutOff("cutOff","position of fermi",ws->var("cutOff_BKG")->getVal());
  cutOff.setConstant(kTRUE);
  RooRealVar beta("beta","width of fermi",ws->var("beta_BKG")->getVal());
  beta.setConstant(kTRUE);
	     		       
  RooFermi fermi("fermi","fermi function",mZZ,cutOff,beta);
  // -------------------- double gauss ---------------------------
  temp="CMS_hzz2l2q_bkg"+btagName[btag]+"p1";
  RooRealVar m("m","m",ws->var(temp.c_str())->getVal());
  m.setConstant(kTRUE);
  temp="CMS_hzz2l2q_bkg"+btagName[btag]+"p2";
  RooRealVar wdth("wdth","wdth",ws->var(temp.c_str())->getVal());
  wdth.setConstant(kTRUE);
  temp="CMS_hzz2l2q_bkg"+btagName[btag]+"p3";
  RooRealVar n("n","n",ws->var(temp.c_str())->getVal());
  n.setConstant(kTRUE);
  temp="CMS_hzz2l2q_bkg"+btagName[btag]+"p4";
  RooRealVar alpha("alpha","alpha",ws->var(temp.c_str())->getVal());
  alpha.setConstant(kTRUE);
  temp="CMS_hzz2l2q_bkg"+btagName[btag]+"p5";
  RooRealVar theta("theta","theta",ws->var(temp.c_str())->getVal());
  theta.setConstant(kTRUE);

  RooCB CB("CB","Crystal ball",mZZ,m,wdth,alpha,n,theta);
  
  RooProdPdf background("background","background",RooArgSet(fermi,CB));
//.........这里部分代码省略.........

void ROCoverlays() {
  TFile *SoftLepton = TFile::Open("FullVariables/histos/AllHistograms.root"); 
  TFile *NoSLInfo = TFile::Open("../../NewVariables/DUSG/LDefault2/histos/AllHistograms.root");

  // loop over pT ranges
  for (Int_t i=7; i<8; i++){
    // Discriminators with SL information
    TH1F *SLC = (TH1F*)SoftLepton->Get(Form("histBDTG_C_Pt%i",i));
    TH1F *SLL = (TH1F*)SoftLepton->Get(Form("histBDTG_light_Pt%i",i));
    
    // Discriminators without SL information
    TH1F *DefC = (TH1F*)NoSLInfo->Get(Form("histBDTG_C_Pt%i",i));
    TH1F *DefL = (TH1F*)NoSLInfo->Get(Form("histBDTG_light_Pt%i",i));
    
    //ROC CvDUSG
    TGraph *SLCvL = (TGraph*)SoftLepton->Get(Form("ROC_C_light_%i",i));
    TGraph *noSLCvL = (TGraph*)NoSLInfo->Get(Form("ROC_C_light_%i",i));
  
    // Normalize Discriminators
    SLC->Scale(1./SLC->Integral());
    SLL->Scale(1./SLL->Integral());
    DefC->Scale(1./DefC->Integral());
    DefL->Scale(1./DefL->Integral());
    
    //Set colors
    noSLCvL->SetLineColor(kBlack);
    SLCvL->SetLineColor(kRed);

    noSLCvL->SetLineWidth(2); 
    SLCvL->SetLineWidth(2);     
    DefC->SetLineWidth(2);
    DefL->SetLineWidth(2);
    SLL->SetLineWidth(2); 
    SLC->SetLineWidth(2);
    
    TCanvas *cv = 0;
    TLegend *l = 0;
    TLatex *tex = 0;

    // Default Discriminators
    DefC->SetLineColor(kBlack);
    DefL->SetLineColor(kRed)

    cv = new TCanvas("cv","cv",800,800);
    cv->SetGridx();
    cv->SetGridy();
    DefL->SetTitle("No Soft Lepton Information");
    DefL->GetYaxis()->SetRangeUser(0,0.03);
    DefL->GetYaxis()->SetTitleOffset(1.6);
    DefL->GetYaxis()->SetTitle("Normalized Events");  
    DefL->GetXaxis()->SetTitle("BDT Output");
    DefL->Draw();
    DefC->Draw("same");
    gStyle->SetOptStat(0);

    l = new TLegend(0.55,0.8,0.90,0.9);
    l->SetBorderSize(0);
    l->SetFillStyle(0);
    l->AddEntry(DefC,"C Discriminator","l");
    l->AddEntry(DefL,"DUSG Discriminator","l");
    l->Draw();
    cv->Update();
    cv->SaveAs(Form("Default_Discriminators_%i.png",i));
    DefL->SetLineColor(kBlack);   

    // Soft Lepton Discriminators
    SLC->SetLineColor(kBlack);
    SLL->SetLineColor(kRed);

    cv = new TCanvas("cv","Reco",800,800);
    cv->SetGridx();
    cv->SetGridy();
    SLC->SetLineColor(kBlack);
    SLL->SetTitle("Soft Lepton Information");
    SLL->GetYaxis()->SetRangeUser(0,0.03);
    SLL->GetYaxis()->SetTitleOffset(1.6);
    SLL->GetYaxis()->SetTitle("Normalized Events");
    SLL->GetXaxis()->SetTitle("BDT Output");
    SLL->Draw();
    SLC->Draw("same");
    gStyle->SetOptStat(0);

    l = new TLegend(0.55,0.8,0.90,0.9);
    l->SetBorderSize(0);
    l->SetFillStyle(0);
    l->AddEntry(SLC,"C Discriminator","l");
    l->AddEntry(SLL,"DUSG Discriminator","l");
    l->Draw();
    cv->Update();
    cv->SaveAs(Form("SL_Discriminators_%i.png",i));
    SLC->SetLineColor(kBlue);


    // C Discriminators
    cv = new TCanvas("cv","C",800,800);
    cv->SetGridx();
    cv->SetGridy();
    SLC->Draw();
    SLC->SetTitle("Charm Discriminators");
    SLC->GetYaxis()->SetRangeUser(0,0.03);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    selection->SetAxisRange(MinX, MaxX, "X");
    selection->GetYaxis()->SetRangeUser(YRangeMin, YRangeMax);
    selection->SetMarkerStyle(20);
    selection->SetMarkerSize(0.9);
    selection->SetMarkerColor(kBlack);
    selection->SetXTitle(xTitle);
    selection->SetYTitle(yTitle);

    prediction->SetAxisRange(MinX, MaxX, "X");
    prediction->GetYaxis()->SetRangeUser(YRangeMin, YRangeMax);
    prediction->SetFillColor(c_LightGray);
    prediction->SetTitle("");
    prediction->SetXTitle(xTitle);
    prediction->SetYTitle(yTitle);

    background->SetAxisRange(MinX, MaxX, "X");
    background->GetYaxis()->SetRangeUser(YRangeMin, YRangeMax);
    background->SetTitle("");
    background->SetLineColor(kRed);
    background->SetLineWidth(2);
    background->SetXTitle(xTitle);
    background->SetYTitle(yTitle);

    TCanvas *c = new TCanvas("ca", "Comparison and ratio of two histos", 700, 700);

    TPad *pad1 = new TPad("pad1a", "pad1a", 0, 0.35, 1, 1);
    pad1->SetLogy();
    pad1->SetBottomMargin(0);
    pad1->Draw();
    pad1->cd();

    prediction->DrawCopy("hist");
    selection->Draw("same");
    prediction->SetFillColor(kAzure-3);
    prediction->SetFillStyle(3354);
    prediction->DrawCopy("e2same");
    background->Scale(36.1/32.6);
    background->Draw("same");

    prediction->SetFillStyle(1001);
    //prediction->SetFillColor(c_LightBrown);
    prediction->SetFillColor(c_LightGray);

    //TLegend* leg1 = new TLegend(0.48, 0.63, 0.95, 0.83);
    TLegend* leg1 = new TLegend(0.44, 0.63, 0.91, 0.83);
    leg1->SetFillStyle(0);
    leg1->SetLineStyle(1);
    leg1->SetTextFont(42);
    //leg1->SetTextSize(0.04);
    leg1->SetTextSize(0.045);
    leg1->AddEntry(prediction, titlePrediction, "lf");
    leg1->AddEntry(selection, titleSelection, "lep");
    leg1->AddEntry(background, titleBackground, "l");
    leg1->Draw("same");

    TPaveText* pt = new TPaveText(0.11, 0.98, 0.95, 0.86, "NDC");
    pt->SetBorderSize(0);
    pt->SetFillStyle(0);
    pt->SetTextAlign(12);
    pt->SetTextSize(0.045);
    pt->AddText(Title);
    pt->AddText(LumiTitle);
    pt->Draw();

    c->cd();
    TPad *pad2 = new TPad("pad2a", "pad2a", 0, 0, 1, 0.35);
    pad2->SetTopMargin(0);
    pad2->Draw();
    pad2->cd();
    TH1F* r = new TH1F(*prediction);
    r->SetTitle("");
    r->SetLabelSize(0.08, "XYZ");
    r->SetLabelOffset(0.01, "XYZ");
    // r->SetTitleSize(0.09, "XYZ");
    r->SetTitleSize(0.125, "XYZ");
    r->SetTitleOffset(0.95, "X");
    r->SetTitleOffset(0.53, "Y");
    // r->SetTitleOffset(0.65, "Y");
    r->SetTickLength(0.05);
    r->SetYTitle(RatioTitle);
    r->SetStats(0);
    r->SetMarkerStyle(20);
    r->SetMarkerSize(0.9);
    r->SetMarkerColor(kBlack);
    r->Reset();
    r->Add(prediction, 1);
    r->Add(background, 1);
    r->Add(selection, -1);
    r->Divide(selection);
    r->SetMaximum(2.2);
    r->SetMinimum(-2.2);
    r->Draw("ep");
    TLine l;
    l.DrawLine(MinX, 0., MaxX+BinWidth, 0.);
    c->cd();
    
    c->SaveAs("compare.pdf");
	
	return 0;
}

void pythia8(Int_t nev  = 100, Int_t ndeb = 1)
{
   const char *p8dataenv = gSystem->Getenv("PYTHIA8DATA");
   if (!p8dataenv) {
      const char *p8env = gSystem->Getenv("PYTHIA8");
      if (!p8env) {
         Error("pythia8.C",
               "Environment variable PYTHIA8 must contain path to pythia directory!");
         return;
      }
      TString p8d = p8env;
      p8d += "/xmldoc";
      gSystem->Setenv("PYTHIA8DATA", p8d);
   }

   const char* path = gSystem->ExpandPathName("$PYTHIA8DATA");
   if (gSystem->AccessPathName(path)) {
         Error("pythia8.C",
               "Environment variable PYTHIA8DATA must contain path to $PYTHIA8/xmldoc directory !");
      return;
   }

// Load libraries
#ifndef G__WIN32 // Pythia8 is a static library on Windows
   if (gSystem->Getenv("PYTHIA8")) {
      gSystem->Load("$PYTHIA8/lib/libpythia8");
   } else {
      gSystem->Load("libpythia8");
   }
#endif
   gSystem->Load("libEG");
   gSystem->Load("libEGPythia8");
// Histograms
   TH1F* etaH = new TH1F("etaH", "Pseudorapidity", 120, -12., 12.);
   TH1F* ptH  = new TH1F("ptH",  "pt",              50,   0., 10.);


// Array of particles
   TClonesArray* particles = new TClonesArray("TParticle", 1000);
// Create pythia8 object
   TPythia8* pythia8 = new TPythia8();

// Configure
   pythia8->ReadString("HardQCD:all = on");


// Initialize

   pythia8->Initialize(2212 /* p */, 2212 /* p */, 14000. /* TeV */);

// Event loop
   for (Int_t iev = 0; iev < nev; iev++) {
      pythia8->GenerateEvent();
      if (iev < ndeb) pythia8->EventListing();
      pythia8->ImportParticles(particles,"All");
      Int_t np = particles->GetEntriesFast();
// Particle loop
      for (Int_t ip = 0; ip < np; ip++) {
         TParticle* part = (TParticle*) particles->At(ip);
         Int_t ist = part->GetStatusCode();
         // Positive codes are final particles.
         if (ist <= 0) continue;
         Int_t pdg = part->GetPdgCode();
         Float_t charge = TDatabasePDG::Instance()->GetParticle(pdg)->Charge();
         if (charge == 0.) continue;
         Float_t eta = part->Eta();
         Float_t pt  = part->Pt();

         etaH->Fill(eta);
         if (pt > 0.) ptH->Fill(pt, 1./(2. * pt));
      }
   }

   pythia8->PrintStatistics();

   TCanvas* c1 = new TCanvas("c1","Pythia8 test example",800,800);
   c1->Divide(1, 2);
   c1->cd(1);
   etaH->Scale(5./Float_t(nev));
   etaH->Draw();
   etaH->SetXTitle("#eta");
   etaH->SetYTitle("dN/d#eta");

   c1->cd(2);
   gPad->SetLogy();
   ptH->Scale(5./Float_t(nev));
   ptH->Draw();
   ptH->SetXTitle("p_{t} [GeV/c]");
   ptH->SetYTitle("dN/dp_{t}^{2} [GeV/c]^{-2}");
 }

//.........这里部分代码省略.........
  TFile * file = new TFile(era+"/"+TStrName+TString(".root"),"update");
  file->cd("");

  TH1D * inputEventsH = new TH1D("inputEventsH","",1,-0.5,0.5);
  TH1D * histWeightsH = new TH1D("histWeightsH","",1,-0.5,0.5);


  TH1D * metAll  = new TH1D("metAll","",400,0,4000);

  TH1D * hDiJetmet = new TH1D("hDiJetmet","",400,0,4000);
  TH1D * hDiJetmass = new TH1D("hDiJetmass","",400,0,4000);
  TH1D * hDiJet1mass = new TH1D("hDiJet1mass","",400,0,4000);
  TH1D * hDiJet2mass = new TH1D("hDiJet2mass","",400,0,4000);
  TH1D * hHT_ = new TH1D("hHT_","",400,0,4000);
  TH1D * PUweightsOfficialH = new TH1D("PUweightsOfficialH","PU weights w/ official reweighting",1000, 0, -1);

 // PILE UP REWEIGHTING - OPTIONS

  if (applyPUreweighting_vertices and applyPUreweighting_official) 
	{std::cout<<"ERROR: Choose only ONE PU reweighting method (vertices or official, not both!) " <<std::endl; exit(-1);}

  // reweighting with vertices

  // reading vertex weights
  TFile * fileDataNVert = new TFile(TString(cmsswBase)+"/src/"+dataBaseDir+"/"+vertDataFileName);
  TFile * fileMcNVert   = new TFile(TString(cmsswBase)+"/src/"+dataBaseDir+"/"+vertMcFileName);

  TH1F * vertexDataH = (TH1F*)fileDataNVert->Get(TString(vertHistName));
  TH1F * vertexMcH   = (TH1F*)fileMcNVert->Get(TString(vertHistName));

  float normVertexData = vertexDataH->GetSumOfWeights();
  float normVertexMc   = vertexMcH->GetSumOfWeights();

  vertexDataH->Scale(1/normVertexData);
  vertexMcH->Scale(1/normVertexMc);


 // reweighting official recipe 
  	// initialize pile up object
  PileUp * PUofficial = new PileUp();
  

  
  if (applyPUreweighting_official) {
    TFile * filePUdistribution_data = new TFile(TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/PileUpDistrib/Data_Pileup_2015D_Nov17.root","read"); 
    TFile * filePUdistribution_MC = new TFile (TString(cmsswBase)+"/src/DesyTauAnalyses/NTupleMaker/data/PileUpDistrib/MC_Spring15_PU25_Startup.root", "read"); 
    TH1D * PU_data = (TH1D *)filePUdistribution_data->Get("pileup");
    TH1D * PU_mc = (TH1D *)filePUdistribution_MC->Get("pileup");
    PUofficial->set_h_data(PU_data); 
    PUofficial->set_h_MC(PU_mc);
  }

  

  int nFiles = 0;
  int nEvents = 0;
  int selEventsAllMuons = 0;
  int selEventsIdMuons = 0;
  int selEventsIsoMuons = 0;



  int nTotalFiles = 0;
  std::string dummy;
  // count number of files --->
  while (fileList0 >> dummy) nTotalFiles++;

//.........这里部分代码省略.........
  //   restriction in centrality
  //*********************************************************************************
  ULong64_t nbytes = 0;
  cout << nEntriesTrack/100000000 << " * 10^8 tracks have to be processed" << endl;
  for (ULong_t i=0; i<nEntriesTrack;i++) {
    nbytes += Track->GetEvent(i);
    /// To do: here you need to add something

		
		
		
		
		
		
		
		
		
		
  }

	
  //*********************************************************************************
  //	Definition of canvas
  //*********************************************************************************
  TCanvas *c1 = new TCanvas("c1","Pt");
  TCanvas *c2 = new TCanvas("c2","nTracksTPC");
  TCanvas *c3 = new TCanvas("c3","nTracksTPCvsCent");
  //*********************************************************************************
  // * Axis labeling for multiplicity histogram and plotting it to an eps file with
  //   the name dependent on the collision system and centrality
  //*********************************************************************************
  hNTPC->Sumw2();
  hNTPC->GetYaxis()->SetTitleOffset(1.3);
  hNTPC->Scale(1./nEntriesPerCent);
  hNTPC->SetXTitle("number TPC tracks");
  hNTPC->SetYTitle("Counts/N_{evt}");	
  c2->cd();
  c2->SetLogy(1);
  HistoSetMarkerAndColor( hNTPC, 20, 1., kBlue, kBlue );
  hNTPC->Draw("e1,p");
  c2->SaveAs(Form("nTracksTPC_%i-%i.eps",startCentrality,endCentrality));	
  
  
  hNTPCvsCent->SetXTitle("number TPC tracks");
  hNTPCvsCent->SetYTitle("Centrality");
  hNTPCvsCent->GetYaxis()->SetTitleOffset(1.3);
  c3->cd(); 
  c3->SetLogy(0);
  c3->SetLogz(1);
  c3->SetRightMargin(0.1);
  hNTPCvsCent->Draw("colz");
  c3->SaveAs(Form("nTracksTPCvsCent_%i-%i.eps",startCentrality,endCentrality));	

	
	
  //*********************************************************************************
  // * Normalisation of the pt-histogram to the number of entries in this 
  //   centrality bin and to the bin with in pt by deviding through the created 
  //   histogram for deltaPt (fDeltaPt)
  // * Labeling of the histograms and plotting it to an eps file with the name
  //   dependent on the collision system and the centralities
  //*********************************************************************************
  cout << "Entries total: " << nEntriesEvent << "\t Entries per Centrality: "  << nEntriesPerCent << endl;
  
  /// To do:The following commands need to be uncommented and you need to add the plotting of the 
  ///        pt-histogram, don't forget to label your axis.

//------------------------------------------------------------------------------
// DrawHistogram
//------------------------------------------------------------------------------
void DrawHistogram(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, 550, 720);
  TCanvas* canvas = new TCanvas(hname, hname, 800, 800);

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

  pad1->SetTopMargin   (0.08);
  pad1->SetBottomMargin(0.02);
  pad1->Draw();
      
  pad2->SetTopMargin   (0.08);
  pad2->SetBottomMargin(0.35);
  pad2->Draw();


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

  pad1->SetLogy(_setLogy);

  THStack* hstack = new THStack(hname, hname);

  TH1F* hist[nProcesses];

  //Save histograms to root file
  TFile* outfile;
  //TString fname = Form("files/%s_%djet.root", hname.Data(),_njet);
  TString fname = "files/0jet_"+hname+".root";
  if(_njet==1) fname = "files/1jet_"+hname+".root";
  outfile = new TFile(fname, "create");
  
  TH1F* data;
  TH1F* top;
  TH1F* tW;
  TH1F* WW;
  TH1F* WZ;
  TH1F* ZZ;
  TH1F* Wg;
  TH1F* WgSMu;
  TH1F* WgSEl;
  TH1F* Wjets;
  TH1F* Zjets;
  TH1F* DYtau;
  TH1F* Zgamma;
  TH1F* ggH;



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

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

    if(ip == iData)   data   = (TH1F*)hist[iData]->Clone("Data");     //data   -> Sumw2();
    if(ip == itt)     top    = (TH1F*)hist[itt]->Clone("top");        //top    -> Sumw2();
    if(ip == itW)     tW     = (TH1F*)hist[itW]->Clone("tW");         //tW     -> Sumw2();
    if(ip == iWW)     WW     = (TH1F*)hist[iWW]->Clone("WW");         //WW     -> Sumw2();
    if(ip == iWZ)     WZ     = (TH1F*)hist[iWZ]->Clone("WZ");   //VV     -> Sumw2();
    if(ip == iZZ)     ZZ     = (TH1F*)hist[iZZ]->Clone("ZZ");         //ZZ     -> Sumw2();
    if(ip == iWg)     Wg     = (TH1F*)hist[iWg]->Clone("Wg");         //Wg     -> Sumw2();
    if(ip == iWgSMu){
      WgSMu  = (TH1F*)hist[iWgSMu]->Clone("WgSMu");         //WgSMu     -> Sumw2();
      hist[iWgSMu]->Scale(1.5);
    }
    if(ip == iWgSEl){
      WgSEl  = (TH1F*)hist[iWgSEl]->Clone("WgSEl");         //WgSel     -> Sumw2();
      hist[iWgSEl]->Scale(1.5);         //WgSel     -> Sumw2();
    }
    if(ip == iWj)     Wjets  = (TH1F*)hist[iWj]->Clone("W+jets");     //Wjets  -> Sumw2();
    if(ip == iDY)     Zjets  = (TH1F*)hist[iDY]->Clone("Z+jets");     //Zjets  -> Sumw2();
    if(ip == iDYtau)  DYtau  = (TH1F*)hist[iDYtau]->Clone("DYtau");   //DYtau  -> Sumw2();
    if(ip == iZgamma) Zgamma = (TH1F*)hist[iZgamma]->Clone("Zgamma"); //Zgamma -> Sumw2();
    if(ip == iH125)   ggH    = (TH1F*)hist[iH125]->Clone("ggH");      //ggH    -> Sumw2();
    
    if (moveOverflow) MoveOverflowBins  (hist[ip], xmin, xmax);
    else              ZeroOutOfRangeBins(hist[ip], xmin, xmax);

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

    if (ip == iWg) {
      //hist[ip]->Scale(0.01);
    }
    
    if (ip == iData) {
//.........这里部分代码省略.........

int ntuAnalyzer(std::string fileName)
{
  setGlobalStyle();
  
  //###############################
  //## run274200 ##
  unsigned int HT250Calo  = 9; //index of DST_HT250_CaloScouting_v   old:1 ref:9
  float HT250Calo_rate = 1928;

  unsigned int HT410PF = 7; //index of DST_HT410_PFScouting_v   old:3 ref:7
  float HT410PF_rate = 294;

  unsigned int MJJ200Calo = 12; //index of DST_DiCaloWideJetMass200_CaloScouting_v
  
  unsigned int HTT200 = 0; //index if L1_HTT200
  unsigned int HTT240 = 1; //index if L1_HTT240
  unsigned int HTT270 = 2; //index if L1_HTT270
  unsigned int HTT280 = 3; //index if L1_HTT280
  unsigned int DoubleJetC100 = 7; //index if L1_DoubleJetC100
  unsigned int DoubleJetC112 = 8; //index if L1_DoubleJetC112
  unsigned int DoubleIsoTau28er = 11; //index if L1_DoubleJetC112

  float instLumi = 0.4; //E34
  float targetLumi = 1; //E34
  float lumiScaleFactor = targetLumi/instLumi;

  float PDRate = 59300; //unprescaled rate of HLTPhysics accordingly to: https://cmswbm2.web.cern.ch/cmswbm2/cmsdb/servlet/DatasetSummary?RUN=274200 and prescale of 9000
  

  unsigned int L1scenario = HTT240;
  //###############################

  TChain* tt = new TChain("MyAnalysis/HLTree");
  tt->Add(fileName.c_str());

  //set branches
  TBranch* b_lumi;
  
  TBranch* b_caloMjj;
  TBranch* b_PFMjj;

  TBranch* b_hltAccept;
  TBranch* b_l1Accept;
  TBranch* b_l1Names;

  TBranch* b_caloJet1Pt;
  TBranch* b_caloJet2Pt;
  TBranch* b_caloJet1Eta;
  TBranch* b_caloJet2Eta;
  TBranch* b_caloDeltaEta;

  TBranch* b_PFJet1Pt;
  TBranch* b_PFJet2Pt;
  TBranch* b_PFJet1Eta;
  TBranch* b_PFJet2Eta;
  TBranch* b_PFDeltaEta;

  int lumi = 0;
  float caloMjj = 0;
  float PFMjj = 0;

  float caloJet1Pt_ = 0;
  float caloJet2Pt_ = 0;
  float caloJet1Eta_ = -999;
  float caloJet2Eta_ = -999;
  float caloDeltaEta_ = -999;

  float PFJet1Pt_ = 0;
  float PFJet2Pt_ = 0;
  float PFJet1Eta_ = -999;
  float PFJet2Eta_ = -999;
  float PFDeltaEta_ = -999;

  
  std::vector<int>* hltAccept = 0;
  std::vector<int>* l1Accept = 0;
  std::vector<string>* l1Names = 0;

  tt->SetBranchAddress("lumi", &lumi, &b_lumi);

  tt->SetBranchAddress("caloMjj", &caloMjj, &b_caloMjj);
  tt->SetBranchAddress("PFMjj", &PFMjj, &b_PFMjj);

  tt->SetBranchAddress("caloJet1Pt", &caloJet1Pt_, &b_caloJet1Pt);
  tt->SetBranchAddress("caloJet2Pt", &caloJet2Pt_, &b_caloJet2Pt);
  tt->SetBranchAddress("caloJet1Eta", &caloJet1Eta_, &b_caloJet1Eta);
  tt->SetBranchAddress("caloJet2Eta", &caloJet2Eta_, &b_caloJet2Eta);
  tt->SetBranchAddress("caloDeltaEta", &caloDeltaEta_, &b_caloDeltaEta);

  tt->SetBranchAddress("PFJet1Pt", &PFJet1Pt_, &b_PFJet1Pt);
  tt->SetBranchAddress("PFJet2Pt", &PFJet2Pt_, &b_PFJet2Pt);
  tt->SetBranchAddress("PFJet1Eta", &PFJet1Eta_, &b_PFJet1Eta);
  tt->SetBranchAddress("PFJet2Eta", &PFJet2Eta_, &b_PFJet2Eta);
  tt->SetBranchAddress("PFDeltaEta", &PFDeltaEta_, &b_PFDeltaEta);

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

  int nentries = tt->GetEntries();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
      if (category == "wGwoGSF") CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf<0.5";
      if (category == "wGwGSFwoPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput<-0.1";
      if (category == "wGwGSFwPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput>-0.1";
    }

    if(discriminator == "-AntiEMed"){
      if (category == "NoEleMatch") CategorySelection = "Tau_GsfEleMatch<0.5"; 
      if (category == "woG") CategorySelection = "Tau_NumGammaCands<0.5"; 
      if (category == "wGwoGSF") CategorySelection = "Tau_NumGammaCands>0.5 && (Tau_HasGsf<0.5 || (Tau_HasGsf>0.5 && Elec_PFMvaOutput>-0.1))";
      if (category == "wGwGSFwoPFMVA")CategorySelection = "Tau_NumGammaCands>0.5 && Tau_HasGsf>0.5 && Elec_PFMvaOutput<-0.1";
    }

  TCut ElecSelection = CategorySelection && PUSelection && ElecPtSelection && ElecAbsEtaSelection && ElecMatchSelection ;
  TCut TauSelection = CategorySelection && PUSelection && TauPtSelection && TauAbsEtaSelection && TauMatchSelection ;
  TCut Selection;
  if (variable.find("Elec")!=std::string::npos)Selection = ElecSelection;
  if (variable.find("Tau")!=std::string::npos)Selection = TauSelection;
  

  TH1F* hVariable   = new TH1F( "hVariable" ,"" , nBins ,xMin, xMax);
  hVariable->SetXTitle(Form("%s",variable.data()));

  if (matching->find("EleMatch")!=std::string::npos){
//     hVariable->SetFillColor(kRed);
//     hVariable->SetFillStyle(3345);
    hVariable->SetLineColor(kRed);
    hVariable->SetLineWidth(2);
  }
  if (matching->find("HadMatch")!=std::string::npos){
//     hVariable->SetFillColor(kBlue);
//     hVariable->SetFillStyle(3354);
    hVariable->SetLineColor(kBlue);
    hVariable->SetLineWidth(2);
  }  
  inputTree->Draw(Form("%s>>hVariable",variable.data()));

  cout<<"Variable plotted : "<<variable<<endl;
  cout<<"Matching applied : "<<matching->data()<<endl;
  cout<<"  Total number of Candidates : "<<hVariable->GetEntries()<<endl;
  inputTree->Draw(Form("%s>>hVariable",variable.data()),Selection);
  cout<<"  Number of Cantidates after selection: "<<hVariable->GetEntries()<<endl;
  hVariable->Scale(1./hVariable->Integral());
  leg->AddEntry(hVariable,Form("%s",matching->data()));

  histograms.push_back(hVariable);
  c1->Clear();
  }
//   double yMin = +1.e+6;
//   double yMax = -1.e+6;
  TH1* refHistogram = histograms.front();
  refHistogram->SetStats(false);
  refHistogram->SetTitle("");
//   refHistogram->SetMinimum(yMin);
//   refHistogram->SetMaximum(yMax);


  if (xAxisTitle == "HoHplusE" ) {
    refHistogram->SetMaximum(1.0);
    refHistogram->SetMinimum(0.01);
    c1->SetLogy();
  }

  if(xAxisTitle == "E_{#gamma}/(P_{in}-P_{out})" ){
    refHistogram->SetMaximum(0.03);
    refHistogram->SetMinimum(0.0);
  }

  if(xAxisTitle == "HadrMva(#tau)" ){
    refHistogram->SetMaximum(0.25);
    refHistogram->SetMinimum(0.0);
  }

  TAxis* xAxis = refHistogram->GetXaxis();
  xAxis->SetTitle(xAxisTitle.Data());
  xAxis->SetTitleOffset(1.15);
  //if(variable.find("AbsEta")!=std::string::npos)xAxis->SetLimits(AbsEtaMin, AbsEtaMax);
  TAxis* yAxis = refHistogram->GetYaxis();
  yAxis->SetTitle(yAxisTitle.Data());
  yAxis->SetTitleOffset(1.30);

  int numHistograms = histograms.size();
  float YMax = 0;
  for ( int iHistogram = 0; iHistogram < numHistograms; ++iHistogram ) {
    TH1* histogram = histograms[iHistogram];
    if(histogram->GetMaximum()>YMax) YMax = histogram->GetMaximum();
  }
  for ( int iHistogram = 0; iHistogram < numHistograms; ++iHistogram ) {
    TH1* histogram = histograms[iHistogram];
    yAxis->SetRangeUser(0.,YMax+0.10*YMax);
    std::string drawOption = "hist";
    if ( iHistogram > 0 ) drawOption.append("same");
    histogram->Draw(drawOption.data());
    leg->Draw();

  }//loop matchings
  string outputName = Form("plots/plotVariablesAntiEMVA/%s/plotVariablesAntiEMVA_v4_%s_%s_%s",category.Data(),category.Data(),variable.data(),Region.Data());
  c1->Print(std::string(outputName).append(".png").data());
  c1->Print(std::string(outputName).append(".pdf").data());

}

//.........这里部分代码省略.........
  vardraw = var+">>"+"Other";
  tree3->Draw(vardraw.c_str(),mthcut.c_str());
  InitHist(other, xtitle.c_str(), ytitle.c_str(),  TColor::GetColor(222,90,106), 1001);
  
  delete canv0;

  //----------------------------------------------------------------------------
  //Print out the yields
  /*  Double_t error=0.0;
  ofstream outfile;
  outfile.open("yields_test.txt");
  outfile << "Yields for the signal region." << std::endl;
  outfile << "VBF   "  << vbf->IntegralAndError(0,vbf->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "TTbar   "  << ttbar->IntegralAndError(0,ttbar->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(0,ztt->GetNbinsX(),error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "other   "  << other->IntegralAndError(0,other->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << "S/sqrt(B)    "  << vbf->Integral()/(other->Integral()+ztt->Integral()+ttbar->Integral()) << endl;
  //--------------------------------------------------------------------------
  //continue outputing
  //outfile << "Ewk total    "  << ewk->IntegralAndError(0,ewk->GetNbinsX(),error) << "+/-" << error << endl;
  outfile << endl << endl << endl;
  outfile << "In the signal region (100,150GeV)  " <<endl;
  outfile << "VBF   "  << vbf->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "TTbar    "  << ttbar->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "Ztt    "  << ztt->IntegralAndError(5,11,error) << "+/-" << error << endl;
  //outfile << "ewk    "  << ewk->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile << "other    "  << other->IntegralAndError(5,11,error) << "+/-" << error << endl;
  outfile.close();*/
  //-----------------------------------------------------------------------
  //Draw the histograms
  TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);
  canv->cd();
  //ewk->Add(other); ttbar->Add(ewk); 
  //ztt->Add(ttbar); vbf->Add(ztt);
  other->Add(ztt);
  //Error band stat
  TH1F* errorBand = (TH1F*)vbf ->Clone("errorBand");
  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;
  //     }
  //}
  Float_t n=other->GetEntries();
  other->Scale(1.0/n);
  n=vbf->GetEntries();
  vbf->Scale(1.0/n);
  //for (Int_t i=1; i<nbins+1; i++) {
  //cout << "i: " << i << " " <<  other->GetBinCenter(i) << " " << other->Integral(0,i) << " " << other->Integral(i,nbins+1) << endl;
  //}
  cout << " other " << other->Integral(10,15) << endl;
  cout << " htt " << vbf->Integral(10,15) << endl;
  //n=ztt->GetEntries();
  //ztt->Scale(1.0/n);
  other->SetMaximum(1.2*std::max(maximum(vbf, 0), maximum(other, 0)));
  //blind(vbf,75,150);
  //ttbar->Draw("histsame");
  //ewk->Draw("histsame");
  other->Draw("hist");
  //ztt->Draw("histsame");
  vbf->Draw("histsame");
  //errorBand->Draw("e2same");
  canv->RedrawAxis();
  //---------------------------------------------------------------------------
  //Adding a legend
  TLegend* leg = new TLegend(0.53, 0.65, 0.95, 0.90);
  SetLegendStyle(leg);
  leg->AddEntry(vbf , "HH"             , "F");
  //leg->AddEntry(ztt  , "bjj-vbf"  , "F" );
  //leg->AddEntry(ttbar, "t#bar{t}"              , "F" );
  //leg->AddEntry(ewk  , "Electroweak"           , "F" );
  leg->AddEntry(other, "Other"                 , "F" );
  //leg->AddEntry(errorBand,"bkg. uncertainty","F");
  leg->Draw();
  //---------------------------------------------------------------------------

  //CMS preliminary 
  const char* dataset = "CMS Preliminary, H#rightarrow#tau#tau, 3.0 ab^{-1} at 14 TeV";
  const char* category = "";
  CMSPrelim(dataset, "#mu#tau_{h}", 0.17, 0.835);
  //CMSPrelim(dataset, "", 0.16, 0.835);
  TPaveText* chan     = new TPaveText(0.52, 0.35, 0.91, 0.55, "tlbrNDC");
  chan->SetBorderSize(   0 );
  chan->SetFillStyle(    0 );
  chan->SetTextAlign(   12 );
  chan->SetTextSize ( 0.05 );
  chan->SetTextColor(    1 );
  chan->SetTextFont (   62 );
  chan->AddText(category);
  chan->Draw();
  //-------------------------------------------------------------------------
  //Save histograms
  canv->Print((var+"_test.png").c_str());

}

void draw_prefit_Sample(std::string inputF, std::string channel, int MaxY, std::string xTitle, std::string nameHisto) {

    //    gStyle->SetOptStat(0);


    SetStyle();
    //    InitSubPad
    TCanvas *canv = MakeCanvas("canv", "histograms", 600, 600);
    float SIGNAL_SCALE = 10;
    bool scaled = true;
    bool log = false;

    TFile* input = new TFile(inputF.c_str());
    //cout<<"1";

    const char* dataset;


    //    std::string channel = "muTau_btag/";
    //    std::string channel = "eleTau_inclusive/";
    THStack hs("hs", "");
    TH1F* data = (TH1F*) input->Get((channel + "data_obs").c_str());
    TH1F* zero = (TH1F*) data->Clone("zero");


    TH1F* QCD = (TH1F*) input->Get((channel + "QCD").c_str());
    InitHist(QCD, "", "", TColor::GetColor(250, 202, 255), 1001);
    hs.Add(QCD);



    TH1F* W = (TH1F*) input->Get((channel + "W").c_str());
    InitHist(W, "", "", 46, 1001);

    TH1F* ZJ = (TH1F*) input->Get((channel + "ZJ").c_str());
//    InitHist(ZJ, "", "", TColor::GetColor(100, 182, 232), 1001);
    W->Add(ZJ);

    TH1F* ZL = (TH1F*) input->Get((channel + "ZL").c_str());
//    InitHist(ZL, "", "", TColor::GetColor(100, 182, 232), 1001);
    W->Add(ZL);


    TH1F* VV = (TH1F*) input->Get((channel + "VV").c_str());
//    InitHist(VV, "", "", TColor::GetColor(100, 182, 232), 1001);
    W->Add(VV);

    hs.Add(W);
    ////    TH1F* ZLL = (TH1F*) input->Get((channel +"ZLL");
    //    InitHist(ttbar, "", "", TColor::GetColor(155, 152, 204), 1001);

    TH1F* TT = (TH1F*) input->Get((channel + "TT").c_str());
    InitHist(TT, "", "", 9, 1001);
    hs.Add(TT);

    TH1F* ZTT = (TH1F*) input->Get((channel + "ZTT").c_str());
    InitHist(ZTT, "", "", TColor::GetColor(248, 206, 104), 1001);

    TH1F* ZTTLow = (TH1F*) input->Get((channel + "ZTT_lowMass").c_str());
    InitHist(ZTTLow, "", "", TColor::GetColor(248, 206, 104), 1001);

    ZTT->Add(ZTTLow);

    hs.Add(ZTT);
    

    InitData(data);

    TH1F* signal = (TH1F*) input->Get((channel + "bba150").c_str());
    signal->Scale(10);
    InitSignal(signal);
    //    signal->SetFillColor(kGreen);
    //    signal->SetLineColor(kGreen);
    hs.Add(signal);


    canv->cd();

    //    const char * MMM = xTitle.c_str();
    //    hs.GetXaxis()->SetLabelSize(9);
    zero->Scale(0);
    zero->GetXaxis()->SetRangeUser(0,60);
    zero->GetXaxis()->SetTitle(xTitle.c_str());
    zero->SetMaximum(MaxY);
    zero->Draw();

    hs.Draw("hsame");
    data->SetBinContent(1,0);
    data->SetBinContent(2,0);
    data->SetBinContent(3,0);
    data->SetBinContent(4,0);
    data->SetBinContent(5,0);
    data->SetBinContent(6,0);
    data->SetBinContent(7,0);
    data->SetBinContent(8,0);
    data->SetBinContent(9,0);
    data->SetBinContent(10,0);
    data->SetBinContent(11,0);
    data->SetBinContent(12,0);
    data->SetBinContent(13,0);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		hST_s->SetBinContent(maxbin,      hST_s->GetBinContent(maxbin-1)+hST_s->GetBinContent(maxbin) );
		hST_t->SetBinContent(maxbin,      hST_t->GetBinContent(maxbin-1)+hST_t->GetBinContent(maxbin) );
		hST_tW->SetBinContent(maxbin,     hST_tW->GetBinContent(maxbin-1)+hST_tW->GetBinContent(maxbin) );
	}
	//syst.
	if (UseOverflow) {
		int maxbin=hTTjetsS->GetNbinsX();
		hTTjetsS->SetBinContent(maxbin,    hTTjetsS->GetBinContent(maxbin+1)+hTTjetsS->GetBinContent(maxbin) );
		hWjetsS->SetBinContent(maxbin,     hWjetsS->GetBinContent(maxbin+1)+hWjetsS->GetBinContent(maxbin) );
		hWjetsSFast->SetBinContent(maxbin, hWjetsSFast->GetBinContent(maxbin+1)+hWjetsSFast->GetBinContent(maxbin) );
		hZjetsS->SetBinContent(maxbin,     hZjetsS->GetBinContent(maxbin+1)+hZjetsS->GetBinContent(maxbin) );
		hZjetsSFast->SetBinContent(maxbin, hZjetsSFast->GetBinContent(maxbin+1)+hZjetsSFast->GetBinContent(maxbin) );
		hQCDS->SetBinContent(maxbin,       hQCDS->GetBinContent(maxbin+1)+hQCDS->GetBinContent(maxbin) );
		hST_sS->SetBinContent(maxbin,      hST_sS->GetBinContent(maxbin+1)+hST_sS->GetBinContent(maxbin) );
		hST_tS->SetBinContent(maxbin,      hST_tS->GetBinContent(maxbin+1)+hST_tS->GetBinContent(maxbin) );
		hST_tWS->SetBinContent(maxbin,     hST_tWS->GetBinContent(maxbin+1)+hST_tWS->GetBinContent(maxbin) );
	}
	if (UseUnderflow) {
		//underflow bin
		int maxbin=1;
		hTTjetsS->SetBinContent(maxbin,    hTTjetsS->GetBinContent(maxbin-1)+hTTjetsS->GetBinContent(maxbin) );
		hWjetsS->SetBinContent(maxbin,     hWjetsS->GetBinContent(maxbin-1)+hWjetsS->GetBinContent(maxbin) );
		hWjetsSFast->SetBinContent(maxbin, hWjetsSFast->GetBinContent(maxbin-1)+hWjetsSFast->GetBinContent(maxbin) );
		hZjetsS->SetBinContent(maxbin,     hZjetsS->GetBinContent(maxbin-1)+hZjetsS->GetBinContent(maxbin) );
		hZjetsSFast->SetBinContent(maxbin,     hZjetsSFast->GetBinContent(maxbin-1)+hZjetsSFast->GetBinContent(maxbin) );
		hQCDS->SetBinContent(maxbin,       hQCDS->GetBinContent(maxbin-1)+hQCDS->GetBinContent(maxbin) );
		hST_sS->SetBinContent(maxbin,      hST_sS->GetBinContent(maxbin-1)+hST_sS->GetBinContent(maxbin) );
		hST_tS->SetBinContent(maxbin,      hST_tS->GetBinContent(maxbin-1)+hST_tS->GetBinContent(maxbin) );
		hST_tWS->SetBinContent(maxbin,     hST_tWS->GetBinContent(maxbin-1)+hST_tWS->GetBinContent(maxbin) );
	}
	
	// scale histograms to 20/pb

	hTTjets->Scale(0.0081); // madgraph
	//hTTjets->Scale(0.0777);//Tauola
	hWjets->Scale(0.0883);
	//hWjetsFast->Scale(0.0091); //fastsim
	hWjetsFast->Scale(hWjets->Integral() / hWjetsFast->Integral()); // scale to FullSim
		
	hZjets->Scale(0.0731);
	hZjetsFast->Scale(hZjets->Integral()/hZjetsFast->Integral()); //scale to FullSim
	hQCD->Scale(0.4003);
	hQCD_WFast = (TH1F*) hWjetsFast->Clone("hQCD_WFast"); //take shape from Wjets
	hQCD_WFast->Scale(hQCD->Integral()/hQCD_WFast->Integral()); //scale to FullSim
	hST_t->Scale(0.003);
	hST_s->Scale(0.0027);
	hST_tW->Scale(0.0034);

	hTTjetsS->Scale(0.0081); // 
	//hTTjetsS->Scale(0.0008); // for fastsim
	hWjetsS->Scale(0.0883);
	//hWjetsS->Scale(0.0091);// from fastsim
	//hWjetsSFast->Scale(hWjetsS->Integral() / hWjetsSFast->Integral()); // scale to FullSim
	//hWjetsSFast->Scale(0.6642); // scaleUP
	//hWjetsSFast->Scale(0.8041); // scaleDown
	//hWjetsSFast->Scale(0.0605); // threshold 5gev
	hWjetsSFast->Scale(0.042); // threshold 20gev
	
	hZjetsS->Scale(0.0731);
	//hZjetsS->Scale(0.0085);// from fastsim
	hZjetsSFast->Scale(hZjetsS->Integral() / hZjetsSFast->Integral()); // scale to FullSim
	hQCDS->Scale(0.4003);
	//hQCDS_WFast = (TH1F*) hWjetsS->Clone("hQCDS_WFast");
	//hQCDS_WFast->Scale(hQCDS->Integral()/hQCDS_WFast->Integral());
	hST_tS->Scale(0.003);
	hST_sS->Scale(0.0027);

//================================================
void DeltaZVsPos(const Int_t save = 0)
{
  THnSparseF *hn = (THnSparseF*)f->Get(Form("mhTrkDzDy_%s",trigName[kTrigType]));
  TList *list = new TList;

  // dz vs BL
  TH2F *hTrkDzVsBL = (TH2F*)hn->Projection(1,3);
  c = draw2D(hTrkDzVsBL,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName[30];
  TH1F *hTrkDzInBL[30];
  Int_t counter = 0;
  for(Int_t i=0; i<30; i++)
    {
      hTrkDzInBL[i] = (TH1F*)hTrkDzVsBL->ProjectionY(Form("hDeltaZ_BL%d",i+1),i+1,i+1);
      if(hTrkDzInBL[i]->GetEntries()>0)
	{
	  legName[counter] = Form("Module %d",i+1);
	  hTrkDzInBL[i]->SetLineColor(color[counter]);
	  list->Add(hTrkDzInBL[i]);
	  counter ++;
	}
    }
  c = drawHistos(list,"TrkDzInBL",Form("%s: #Deltaz of matched track-hit pairs in backleg;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInBL[1]->GetMaximum(),kFALSE,kTRUE,legName,kTRUE,"",0.15,0.25,0.2,0.88,kFALSE,0.04,0.04,kFALSE,1,kTRUE,kFALSE);
  TLine *line = GetLine(0,0,0,1.1*hTrkDzInBL[1]->GetMaximum(),1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_BL_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  // dz vs Mod
  TH2F *hTrkDzVsMod = (TH2F*)hn->Projection(1,4);
  c = draw2D(hTrkDzVsMod,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_vs_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  TH1F *hMthMod = (TH1F*)hTrkDzVsMod->ProjectionX("hMthMod");
  hMthMod->Sumw2();
  hMthMod->Scale(1./hMthMod->Integral());
  TH2F *hMtdHitMap = (TH2F*)f->Get(Form("mhMtdHitMap_%s",trigName[kTrigType]));
  TH1F *htmp = (TH1F*)hMtdHitMap->ProjectionY("hHitMod_finebin");
  htmp->Rebin(12);
  TH1F *hMtdHitMod = new TH1F(Form("hMtdHitMod_%s",trigName[kTrigType]),"# of MTD hits per module;module",5,1,6);
  for(int i=0; i<hMtdHitMod->GetNbinsX(); i++)
    {
      hMtdHitMod->SetBinContent(i+1,htmp->GetBinContent(i+1));
      hMtdHitMod->SetBinError(i+1,htmp->GetBinError(i+1));
    }
  hMtdHitMod->Scale(1./hMtdHitMod->Integral());
  list->Clear();
  list->Add(hMthMod);
  list->Add(hMtdHitMod);
  TString legName3[2] = {"Matched good hits","All good hits"};
  c = drawHistos(list,"MtdHitMod",Form("%s: MTD hits per module;module;probability",trigName[kTrigType]),kFALSE,0,5,kTRUE,0,0.5,kFALSE,kTRUE,legName3,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sCompMtdHitMod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }

  list->Clear();
  TString legName2[5];
  TH1F *hTrkDzInMod[5];
  for(Int_t i=0; i<5; i++)
    {
      hTrkDzInMod[i] = (TH1F*)hTrkDzVsMod->ProjectionY(Form("hDeltaZ_Mod%d",i+1),i+1,i+1);
      legName2[i] = Form("Module %d",i+1);
      list->Add(hTrkDzInMod[i]);
    }
  c = drawHistos(list,"TrkDzInMod",Form("%s: #Deltaz of matched track-hit pairs in module;#Deltaz (cm)",trigName[kTrigType]),kTRUE,-100,100,kTRUE,0,1.2*hTrkDzInMod[3]->GetMaximum(),kFALSE,kTRUE,legName2,kTRUE,"",0.15,0.25,0.6,0.88,kTRUE);
  TLine *line = GetLine(0,0,0,hTrkDzInMod[3]->GetMaximum()*1.05,1);
  line->Draw();
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.pdf",run_type,run_cfg_name.Data(),trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/%sDeltaZ_in_Mod_%s.png",run_type,run_cfg_name.Data(),trigName[kTrigType]));
    }
}

//.........这里部分代码省略.........
	if(j>1) num2 += data->GetBinContent(j); // ... +1,2,3,4... jets
	if(j>2) num3 += data->GetBinContent(j); // ... +2,3,4... jets
	if(j>3) num4 += data->GetBinContent(j); // ..    if(str=="jet_pT"){
	if(j>4) num5 += data->GetBinContent(j); // ... +4... jets
      }
      cout << "\n";
      cout << data->GetNbinsX() <<" Number of bins of the zYieldVsjets histo\n";
      printf("Number of Z+n jet %i --- >1 %i --- >2 %i --- >3 %i --- >4 %i \n",num1,num2,num3,num4,num5);
      cout << "\n";
    }    

    //======================
    // Z + jets signal
    mcf->cd("validationJEC");
    if (isMu) mcf->cd("validationJECmu/");

    if (isAngularAnalysis) {
      mcf->cd("validationJEC/");
      if (isMu) mcf->cd("validationJECmu/");
    }

    TH1F* mc;
    gDirectory->GetObject(plot.c_str(),mc);
    TH1F * hsum;
    if(mc){
      hsum =  (TH1F*) mc->Clone();
      hsum->SetTitle("hsum");
      hsum->SetName("hsum");
      hsum->Reset();

      Double_t mcint = mc->Integral();
      mc->SetFillColor(kRed);
      mc->Sumw2();
      if(lumiweights==0) mc->Scale(dataint/mcint);
		
      // Blocco da propagare negli altri MC
      if(zNumEvents>0.){
	if(lumiweights==1) {
	  if (WholeStat){
	    if (lumiPixel) mc->Scale( dataLumi2011pix / (zNumEvents / zjetsXsect));
	    else mc->Scale( dataLumi2011 / (zNumEvents / zjetsXsect));
	  }
	  else{
	    if (RunA){
	      if (lumiPixel) mc->Scale( dataLumi2011Apix / (zNumEvents / zjetsXsect));
	      else mc->Scale( dataLumi2011A / (zNumEvents / zjetsXsect));
	    }
	    if (!RunA){
	      if (lumiPixel) mc->Scale( dataLumi2011Bpix / (zNumEvents / zjetsXsect));
	      else mc->Scale( dataLumi2011B / (zNumEvents / zjetsXsect));
	    }
	  }
	}
      }
      else {
	if(lumiweights==1) mc->Scale(zjetsScale);
      }

      // fin qui
      if(lumiweights==1) mc->Scale(1./zwemean);  // perche' i Weights non fanno 1...
      mc->Rebin(rebin);
      if(lumiweights==0) mc->Draw("HISTO SAMES");
      hsum->Rebin(rebin);
      hsum->Add(mc);
      legend->AddEntry(mc,"Z+jets","f");
    }