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

void drawHistos(TCanvas * C, TString filename, TString category, TTree* Tmine, TTree* Tother,TString var, int nbins, float xmin, float xmax, TString selection, TString myGroup, TString myRootFile, TString group, TString groupRootFile,TString mySel="1",TString groupSel="1"){
  TH1F* Hmine = new TH1F(TString("Hmine")+var,"",nbins,xmin,xmax); 
  Hmine->GetYaxis()->SetTitle(category);
  Hmine->GetXaxis()->SetTitle(var);
  Hmine->SetLineColor(1);
  Hmine->SetStats(0);
  TH1F* Hother = new TH1F(TString("Hother")+var,"",nbins,xmin,xmax); 
  Hother->GetYaxis()->SetTitle(category);
  Hother->GetXaxis()->SetTitle(var);
  Hother->SetLineColor(2);
  Hother->SetStats(0);

  TText TXmine;
  TXmine.SetTextColor(1);
  TXmine.SetTextSize(.04);
  TText TXother;
  TXother.SetTextColor(2);
  TXother.SetTextSize(.04);

  Tmine->Draw(var+">>"+Hmine->GetName(),selection+"*("+mySel+")");
  Tother->Draw(var+">>"+Hother->GetName(),selection+"*("+groupSel+")");

  ////Draw one histogram on top of the other
  C->Clear();
  //Hmine->Scale(1./Hmine->Integral());
  //Hother->Scale(1./Hother->Integral()); 
  //Hother->Scale(968134./688134.); //GGH e-tau
  if(Hmine->GetMaximum()>Hother->GetMaximum())
    Hmine->GetYaxis()->SetRangeUser(0,Hmine->GetMaximum()*1.1);
  else Hmine->GetYaxis()->SetRangeUser(0,Hother->GetMaximum()*1.1);
  Hmine->Draw("hist");
  Hother->Draw("histsame");

//   ///Draw the difference of the historgrams
//   TH1F*HDiff=(TH1F*)Hmine->Clone("HDiff");
//   HDiff->Add(Hother,-1);
//   int max= abs(HDiff->GetMaximum())>abs( HDiff->GetMinimum()) ?   abs(HDiff->GetMaximum()): abs( HDiff->GetMinimum());
//   HDiff->GetYaxis()->SetRangeUser(-2*(max>0?max:1),2*(max>0?max:1));
//   HDiff->Draw("hist");
//   TLine line;
//   line.DrawLine(HDiff->GetXaxis()->GetXmin(),0,HDiff->GetXaxis()->GetXmax(),0);

  //Print the integrals of the histograms a the top
  //TXmine.DrawTextNDC(.2,.965,myGroup+"_"+myRootFile+": "+(long)(Hmine->Integral(0,Hmine->GetNbinsX()+1)));
  //TXother.DrawTextNDC(.2,.93,group+"_"+groupRootFile+": "+(long)(Hother->Integral(0,Hother->GetNbinsX()+1)));
  TXmine.DrawTextNDC(.2,.965,myGroup+" : "+(long)(Hmine->Integral(0,Hmine->GetNbinsX()+1)));
  TXother.DrawTextNDC(.2,.93,group+": "+(long)(Hother->Integral(0,Hother->GetNbinsX()+1)));
  C->Print(filename);

  delete Hmine;
  delete Hother;
}

void pv_dist(const string& fFile, const string& fTitle, const string& fNameExt) {
  
   TH1F *h;
    
   TFile file(fFile.c_str());
   TDirectoryFile *dir = (TDirectoryFile*)file.Get("offsetAnalysis");
   TDirectoryFile *subDir = (TDirectoryFile*)dir->Get("PrimaryVertices");

   h = (TH1F*)subDir->Get("h_NofPVs"); 
   
   string name = h->GetName();
   string fileName = name + "__" + fNameExt + ".png";
   
   h->SetTitle(fTitle.c_str());
   
   TCanvas *c = new TCanvas("c","",1120,800);
   c->cd();
     
   h->SetLineWidth(2);
   h->Draw();
   
   c->SetLogy();
   c->SaveAs(fileName.c_str());
   
   delete c;
}

void KVElasticCountRates::PrintResults(Double_t beam_intensity)
{
   // Print mean energy deposit & counting rate for given beam intensity in particles per second

   TIter it(&fHistos);
   TH1F* h;
   fRates.clear();

   std::vector<count_rate> count_rates;

   while ((h = (TH1F*)it())) {
      TString name = h->GetName();
      if (!name.EndsWith("_dW") && !name.EndsWith("_map")) {
         TH2F* map = (TH2F*)fHistos.FindObject(name + "_map");
         double rate = h->Integral() * fAtomicDensity * beam_intensity * fVolume / fNtirages;
         double emean = h->GetMean();
         KVDetector* det = gMultiDetArray->GetDetector(name);
         double fluence = rate / det->GetEntranceWindowSurfaceArea();
         double dissipation = emean * rate / det->GetEntranceWindowSurfaceArea();
         count_rates.push_back(
            count_rate(name, rate, emean, map->GetMean(), map->GetMean(2), fluence, dissipation)
         );
         fRates[name.Data()] = KVElasticCountRate(rate, emean, fluence, dissipation);
      }
   }
   std::sort(count_rates.begin(), count_rates.end(), compare_count_rates);

   for (std::vector<count_rate>::iterator it = count_rates.begin(); it != count_rates.end(); ++it) {
      it->print();
   }
}

void make_histos_syst_rawyield(TString file_syst, TString file_default, TString out_tag){

  TFile *f1 = new TFile(file_syst.Data(),"read");
  TFile *f2 = new TFile(file_default.Data(),"read");
  TDirectoryFile *dir1 = (TDirectoryFile*)(f1->Get("effunf"));
  TDirectoryFile *dir2 = (TDirectoryFile*)(f2->Get("effunf"));
  TList *list = dir1->GetListOfKeys();

  TFile *f = new TFile(Form("plots/ratiosyst_%s.root",out_tag.Data()),"recreate");

  for (int i=0; i<list->GetSize(); i++){
    TString name = dir1->GetListOfKeys()->At(i)->GetName();
    if (!(name.Contains("hreco_"))) continue;
    TObject *obj1 = dir1->Get(name.Data());
    assert(obj1);
    TObject *obj2 = dir2->Get(name.Data());
    assert(obj2);
    TString newname = name;
    newname.Append("_ratiosyst");
    if (name.EndsWith("_0")) newname.ReplaceAll("_0_","_EBEB_");
    if (name.EndsWith("_1")) newname.ReplaceAll("_1_","_EBEE_");
    if (name.EndsWith("_2")) newname.ReplaceAll("_2_","_EEEE_");
    TH1F *h = (TH1F*)(((TH1F*)obj1)->Clone(newname.Data()));
    h->SetTitle(h->GetName());
    h->Divide((TH1F*)obj2);
    for (int j=0; j<h->GetNbinsX(); j++) h->SetBinError(j+1,0);
    for (int j=0; j<h->GetNbinsX(); j++) h->SetBinContent(j+1,1+fabs(1-h->GetBinContent(j+1)));
    f->cd();
    h->Write();
  }

}

void example(double E0 = 50, int nevents = 100000)
{
   TStopwatch timer;

   // compound nucleus = carbon
   KVNucleus CN(6, 12);
   CN.SetExcitEnergy(E0);

   // decay products
   KVEvent decay;
   KVNucleus* n = decay.AddParticle();
   n->SetZandA(1, 2);
   n = decay.AddParticle();
   n->SetZandA(2, 4);
   n = decay.AddParticle();
   n->SetZandA(3, 6);

   MicroStat::mdweight gps;
   Double_t etot = E0 + decay.GetChannelQValue();
   Double_t total_mass = decay.GetSum("GetMass");

   if (etot <= 0) {
      printf("Break-up channel is not allowed\n");
      return;
   }
   gps.SetWeight(&decay, etot);
   gps.initGenerateEvent(&decay);

   std::cout << "Edisp = " << etot << " MeV" << std::endl;
   KVHashList histos;
   TH1F* h;

   while ((n = decay.GetNextParticle())) {
      Double_t kappa = total_mass / (total_mass - n->GetMass());
      std::cout << n->GetSymbol() << " : max KE = " << 1. / kappa << " * " << etot << " MeV" << std::endl;
      std::cout << n->GetSymbol() << " : m/M = " << n->GetMass() / total_mass << " k = " << kappa << std::endl;
      histos.Add(h = new TH1F(n->GetSymbol(), Form("Kinetic energy of %s", n->GetSymbol()), 200, 0, etot));
      h->Sumw2();
   }
   KVEvent event;

   while (nevents--) {
      gps.GenerateEvent(&decay, &event);
      while ((n = event.GetNextParticle()))((TH1F*)histos.FindObject(n->GetSymbol()))->Fill(n->GetE());
      gps.resetGenerateEvent();
   }

   TIter it(&histos);

   while ((h = (TH1F*)it())) {

      KVNucleus part(h->GetName());
      new TCanvas;
      FitEDist(h, etot, decay.GetMult(), total_mass, part.GetMass());

   }

   timer.Print();
}

//================================================
void DeltaZVsPt(const Int_t save = 0)
{
  THnSparseF *hn = (THnSparseF*)f->Get(Form("mhTrkDzDy_%s",trigName[kTrigType]));
  TH2F *hTrkDzVsPt = (TH2F*)hn->Projection(1,0);
  c = draw2D(hTrkDzVsPt,Form("%s: #Deltaz of matched track-hit pairs",trigName[kTrigType]));
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/DeltaZ_vs_pt_%s.pdf",run_type,trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/DeltaZ_vs_pt_%s.png",run_type,trigName[kTrigType]));
    }

  Double_t pt_cut = 1;
  hTrkDzVsPt->GetXaxis()->SetRangeUser(pt_cut+0.1,100);
  TH1F *hMthDz = (TH1F*)hTrkDzVsPt->ProjectionY(Form("hTrkDzVsPt_%s_proj",trigName[kTrigType]));
  hMthDz->SetTitle(Form("%s: #Deltaz of matched track-hit pairs (p_{T}>%1.1f GeV/c);#Deltaz (cm)",trigName[kTrigType],pt_cut));
  TH1F *hClone = (TH1F*)hMthDz->Clone(Form("%s_clone",hMthDz->GetName()));
  c = draw1D(hClone,"",kFALSE,kFALSE);
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/DeltaZ_%s.pdf",run_type,trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/DeltaZ_%s.png",run_type,trigName[kTrigType]));
    }

  Double_t range = 50;
  TF1 *func = new TF1("func","gaus(0)+gaus(3)",-1*range,range);
  func->SetParameters(10000,0,10,1000,0,40);
  c = FitDeltaZ(hMthDz,func,range,20.);
  if(save) 
    {
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/FitDz_Pt%1.0f_%s.pdf",run_type,pt_cut,trigName[kTrigType]));
      c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/FitDz_Pt%1.0f_%s.png",run_type,pt_cut,trigName[kTrigType]));
    }

  // pt dependence
  Double_t pt_cuts[5] = {1,2,3,5,20};
  for(Int_t i=0; i<4; i++)
    {
      hTrkDzVsPt->GetXaxis()->SetRangeUser(pt_cuts[i]+0.1,pt_cuts[i+1]-0.1);
      TH1F *htmp = (TH1F*)hTrkDzVsPt->ProjectionY(Form("hTrkDz_pt%1.0f-%1.0f_%s",pt_cuts[i],pt_cuts[i+1],trigName[kTrigType]));
      htmp->SetTitle(Form("%s: #Deltaz of matched track-hit pairs (%1.0f < p_{T} < %1.0f GeV/c);#Deltaz (cm)",trigName[kTrigType],pt_cuts[i],pt_cuts[i+1]));

      TF1 *func = new TF1(Form("func_pt%1.0f-%1.0f",pt_cuts[i],pt_cuts[i+1]),"gaus(0)+gaus(3)",-1*range,range);
      if(i==0) func->SetParameters(100,0,100,1000,0,10);
      if(i==1) func->SetParameters(1000,0,15,1000,0,60);
      if(i==2) func->SetParameters(1000,0,15,1000,0,60);
      if(i==3) func->SetParameters(1000,0,60,1000,0,15);
      c = FitDeltaZ(htmp,func,range,20.);
      if(save) 
	{
	  c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/FitDz_Pt%1.0f_%1.0f_%s.pdf",run_type,pt_cuts[i],pt_cuts[i+1],trigName[kTrigType]));
	  c->SaveAs(Form("~/Work/STAR/analysis/Plots/%s/ana_Match/FitDz_Pt%1.0f_%1.0f_%s.png",run_type,pt_cuts[i],pt_cuts[i+1],trigName[kTrigType]));
	}
    }
 
}

void superimposeHistos()
{
  TFile* bFile = TFile::Open("Electron_In_Jets_900GeV_bJets.root");
  TFile* cFile = TFile::Open("Electron_In_Jets_900GeV_cJets.root");
  TFile* udsgFile = TFile::Open("Electron_In_Jets_900GeV_udsgJets.root");
  
  TIter next(bFile->GetListOfKeys());
  TFile* newFile = new TFile("testFile.root", "RECREATE");
  while(TKey* key = (TKey*)next())
  {
    TH1F* bHist = (TH1F*)bFile->Get(key->GetName());
    bHist->SetFillColor(2);
    
    TH1F* cHist = (TH1F*)cFile->Get(key->GetName());
    cHist->SetFillColor(3);
    
    TH1F* udsgHist = (TH1F*)udsgFile->Get(key->GetName());
    udsgHist->SetFillColor(4);
    
    THStack* stack = new THStack(bHist->GetName(), bHist->GetTitle());
    stack->Add(udsgHist, "hist ][");
    stack->Add(cHist, "hist ][");
    stack->Add(bHist, "hist ][");
    
    TLegend* legend = new TLegend(0.5, 0.68, 0.88, 0.88);
    legend->AddEntry(bHist, "b-Jets");
    legend->AddEntry(cHist, "c-Jets");
    legend->AddEntry(udsgHist, "udsg-Jets");
    
    TCanvas* canvas = new TCanvas(bHist->GetName());
    stack->Draw();
    stack->GetXaxis()->SetTitle(bHist->GetXaxis()->GetTitle());
    legend->Draw();
    canvas->Write(canvas->GetName());
  }
  
  newFile->Close();
  bFile->Close();
  cFile->Close();
  udsgFile->Close();
}

//loads single track pT efficiency from root file
TH1F *GetEfficiencyPt(TString effLoc){
  TFile *fIn = 0;
  TH1F *hEffPt = 0; 

  if(!fIn)fIn = TFile::Open(effLoc.Data());
  if(!fIn)Printf("%s%d no input data",(char*)__FILE__,__LINE__);
  if(!hEffPt)hEffPt = (TH1F*)fIn->Get("hSingleTrackEffPt");
  if(!hEffPt)Printf("%s%d no single track efficiency spectrum available",(char*)__FILE__,__LINE__);
  gROOT->cd();

  TH1F *hEffPtClone = (TH1F*)hEffPt->Clone(hEffPt->GetName()); 
  fIn->Close();
  return hEffPtClone;


}

void PileUpMaker() {
  // Access data file
  TFile *file = TFile::Open("ZmumuGammaNtuple_Full2012_MuCorr.root");
  TTree *Datatree = (TTree*)file->Get("ZmumuGammaEvent");

  // Set up histogram for data pile up
  TH1F *pileupraw = new TH1F("pileup","",40,0,80);
 
  // Fill pile up distribution
  Int_t nentries = Datatree->GetEntriesFast();
  Int_t nbytes = 0; 
  for (Int_t i=0; i<nentries;i++) {
    nbytes += Datatree->GetEvent(i);
    UInt_t NVertices;
    Datatree->SetBranchAddress("NVertices",&NVertices);
    pileupraw->Fill(float(NVertices)/0.7);
  }

  // Normalize the distribution
  TH1F* PU = (TH1F*)pileupraw->Clone("pileup");
  PU->Scale(1./PU->Integral());

  // Plot both the normalized and unnormalized pile up
  TCanvas *cv = 0;
  cv= new TCanvas("cv","cv",800, 600);
  pileupraw->Draw("");
  cv->SaveAs("DataPU.pdf");

  cv= new TCanvas("cv","cv",800,600);
  PU->Draw("");
  cv->SaveAs("DataPU_normalized.pdf");

  // Save the normalized distribution to a root file for weighting MC
  TFile *file = TFile::Open("2012_PhosphorData_PileUp.root","UPDATE");
  file->cd();
  file->WriteTObject(PU,PU->GetName(),"WriteDelete");
  file->Close();
  delete file;
}

// Called just after the main event loop
// Can be used to write things out, dump a summary etc
// Return non-zero to indicate a problem
int IslandAmplitude::AfterLastEntry(TGlobalData* gData,const TSetupData *setup){

  // Print extra info if we're debugging this module:
  if(Debug()){
     cout<<"-----IslandAmplitude::AfterLastEntry(): I'm debugging!"<<endl;
  }

  double run_norm = fAmpNorm->Integral(0,-1);
  for(mapSH_t::iterator it = fAmpHist.begin(); it != fAmpHist.end(); ++it)
    {
      TH1F* h = it->second;
      TObject* obj = h->Clone((std::string(h->GetName()) + "_RunNorm").c_str());
      TH1* hn = static_cast<TH1*>(obj);
      hn->SetTitle((std::string(h->GetTitle()) + " (run normalized)").c_str());
      hn->Scale(1.0/run_norm);
    }

  for(mapSH_t::iterator it = fAmpHistNorm.begin(); it != fAmpHistNorm.end(); ++it)
    it->second->Scale(1.0/fNProcessed);

  return 0;
}

void plotProducer(const string& fFile, const string& fPlot, const string& fNameExt, const string& fExt, const string& fSwitch) {
  
   TH1F *h;
    
   TFile file(fFile.c_str());
   TDirectoryFile *dir = (TDirectoryFile*)file.Get("offsetAnalysis");
   TDirectoryFile *subDir = (TDirectoryFile*)dir->Get("PrimaryVertices");

   h = (TH1F*)subDir->Get(fPlot.c_str());

   string name = h->GetName();
   string fileName = name + "__" + fNameExt + "." + fExt;
   
   TCanvas *c = new TCanvas("c","",1120,800);
   c->cd();
     
   //h->SetLineWidth(2);
   h->Draw();
   
   if ( fSwitch.find("log") != string::npos ) c->SetLogy();
   c->SaveAs(fileName.c_str());
   
   delete c;
}

int main(int argc, char* argv[])
{
        string htitle = "Events rate, cloud ";
        htitle = htitle+string(argv[1]);
	string drootpath;
    
        if (strcmp(argv[2],"ttreecache")==0)
	{
                htitle = htitle + ", TTreeCache enabled, protocol: WebDAV";
		drootpath = "../../../Histograms/TDavix/";
        }
	else if (strcmp(argv[2],"xrootd")==0)
        {
		htitle = htitle + ", TTreeCache enabled xrootd";
		drootpath = "../../../Histograms/xrootd/";
        }
	else
        {
	        htitle = htitle + ", TTreeCache disabled, protocol: WebDAV ";
		drootpath = "../../../Histograms/Davix/";
	}

	string cloudpath = argv[1] ;
	cloudpath = cloudpath + "/";
	DIR* dcloudrep = opendir(string2char(drootpath+cloudpath));
	struct dirent* dcloudfichier = NULL;
	string allgrid;
	allgrid = string(argv[3]);

	vector<ratioMatrixValue> th2array;
	
	while ((dcloudfichier = readdir(dcloudrep))!=0)
	{
		if (dcloudfichier->d_name[0]=='.')
			continue;
		string dcloudpath = cloudpath + dcloudfichier->d_name + "/";
		//DIR* drootfile = opendir(string2char(drootpath+dcloudpath));
		//struct dirent* drootfichier = NULL;

		string dfilepath = drootpath + dcloudpath + allgrid + "/" + dcloudfichier->d_name + ".root";
		ratioMatrixValue column;

		TFile *dfr = TFile::Open(string2char(dfilepath));
		if (dfr==NULL)
			continue;

		column.analy = dcloudfichier->d_name;

		for (int i=0;i<dfr->GetListOfKeys()->GetEntries();i++)
		{
			string dstrRead = string2char(dfr->GetListOfKeys()->At(i)->GetName());
			TH1F* dhist = (TH1F*)dfr->Get(string2char(dstrRead));
			if (dhist==NULL)
				continue;
			string otherTitle = string(dhist->GetName()) + "No Zero";
			TH1F* histSansZeros = new TH1F(string2char(otherTitle),dhist->GetTitle(),3000.,0,3000);
			for ( int l=2;l<3000;l++)
				for (int k=0;k<dhist->GetBinContent(l);k++)
					histSansZeros->Fill(l-1);
			column.datadisk.push_back(string2char(dstrRead));
			column.ratio.push_back(histSansZeros->GetMean());
			column.error.push_back(histSansZeros->GetRMSError());
		}
		th2array.push_back(column);
	}



        TCanvas * c = new TCanvas("c", "c");
        c->SetLeftMargin(0.20);
        c->SetBottomMargin(0.15);
        c->SetWindowSize(2100,1500);

	/*for (unsigned int i=0;i<th2array.size();i++)
	{
		cout << th2array[i].analy << "\t";
		for (unsigned int j=0;j<th2array[i].datadisk.size();j++)
			cout << th2array[i].datadisk[j] << "\t" ;
		cout << endl;
	}*/

	cout << "--------------------------" << endl;

	th2array = ordenate(th2array,string(argv[1]));

	cout << "--------------------------" << endl;

	/*for (unsigned int i=0;i<th2array.size();i++)
	{
		cout << th2array[i].analy << "\t";
		for (unsigned int j=0;j<th2array[i].datadisk.size();j++)
			cout << th2array[i].datadisk[j] << "\t" ;
		cout << endl;
	}*/

	unsigned int datadisk_number=0;
	unsigned int analy_number=0;
	for (unsigned int i=0;i<th2array.size();i++)
	{
		if (datadisk_number < th2array[i].datadisk.size())
//.........这里部分代码省略.........

int BDTInterpolation(std::string inFileName,bool Diagnose=false, bool doNorm=true, bool doSidebands=false){

  std::cout << getTime() << std::endl;

  system("rm -r plots");
  system("mkdir plots");

  //gROOT->SetBatch();
  gStyle->SetOptStat(0);

  // input flags
  bool all=0;
 
  if (Diagnose) {
    std::cout << "Diagnostics turned on \n Output plots available in \"plots\" directory \n Diagnostic log available in \"plots/BDTInterpolationDiagnostics.txt\"" << std::endl;
    all = true;
  }
  else std::cout << "Diagnostics turned off" << std::endl;
  if (doNorm) std::cout << "Normalization turned on" << std::endl;
  else std::cout << "Normalization turned off" << std::endl;
  if (doSidebands) std::cout << "Background model from sidebands turned on" << std::endl;
  else std::cout << "Background model from sidebands turned off" << std::endl;

  TFile *inFile = new TFile(inFileName.c_str());
  //TFile *inFile = new TFile("/vols/cms02/nw709/hgg/src_cvs/oct13/CMSSW_4_2_8/src/HiggsAnalysis/HiggsTo2photons/h2gglobe/Macros/CMS-HGG_1658pb_mva.root");
  //TFile *inFile = new TFile("RefWorkspaces/CMS-HGG_1658pb_mva.root");
  TFile *outFile = new TFile(Form("%s_interpolated.root",inFileName.c_str()),"RECREATE");
  TFile *systTestF = new TFile("systTest.root","RECREATE");
  ofstream diagFile("plots/BDTInterpolationDiagnostics.txt");

  const int nBDTs=2;
  const int nMasses=8;
  const int nProds=4;
  std::string BDTtype[nBDTs] = {"ada","grad"};
  std::string BDTmasses[nMasses] = {"110.0","115.0","120.0","125.0","130.0","135.0","140.0","150.0"};
  std::string productionTypes[nProds] = {"ggh","vbf","wzh","tth"};
// ----- else can just get rebinned histograms straight out of workspace

  std::cout << "Extracting histograms from workspace........." << std::endl;

  diagFile << "Diagnostics for Signal Interpolation run at " << getTime() << std::endl;
  diagFile << "---------------------------------------------------------------------" << std::endl;
  diagFile << "Following orginal histograms rewritten into new workspace: " << std::endl;
  
  // ------ stuff for interpolation systematic -------
  TH1F *systHists120[2][3];
  TH1F *systHists135[2][3];

  std::string syst120mass[3] = {"120.0","115.0","125.0"};
  std::string syst135mass[3] = {"135.0","130.0","140.0"};
  // ---------------------------------------------------
  
  // make plots of background model from sidebands
  if (doSidebands){
    for (int bdt=0; bdt<nBDTs; bdt++){
      for (double mass=110.; mass<=150.0; mass+=0.5){
        TList *bkgModelList = new TList();
        TH1F *bkgModel[7];
        bkgModel[0] = (TH1F*)inFile->Get(Form("th1f_bkg_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[1] = (TH1F*)inFile->Get(Form("th1f_bkg_3low_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[2] = (TH1F*)inFile->Get(Form("th1f_bkg_2low_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[3] = (TH1F*)inFile->Get(Form("th1f_bkg_1low_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[4] = (TH1F*)inFile->Get(Form("th1f_bkg_1high_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[5] = (TH1F*)inFile->Get(Form("th1f_bkg_2high_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        bkgModel[6] = (TH1F*)inFile->Get(Form("th1f_bkg_3high_%s_%3.1f_cat0",BDTtype[bdt].c_str(),mass));
        for (int i=0; i<7; i++) {
          bkgModelList->Add(bkgModel[i]);
        }
        std::string name(Form("%3.1f",mass));
        plotBkgModel(bkgModelList,name);
      }
    }
  }
  
  // write original histograms in out file
  TList *HistList = inFile->GetListOfKeys();
  for (int j=0; j<HistList->GetSize(); j++){
    TH1F *temp = (TH1F*)inFile->Get(HistList->At(j)->GetName());
    TString name = temp->GetName();

    // store stuff for interpolation systematic
    for (int bdt=0; bdt<nBDTs; bdt++){
      for (int syst=0; syst<3; syst++){
        if (name.Contains(("th1f_sig_"+BDTtype[bdt]+"_"+productionTypes[0]+"_"+syst120mass[0]+"_"+syst120mass[syst]).c_str()) && !name.Contains("sigma")){
          systHists120[bdt][syst] = (TH1F*)inFile->Get(name.Data());
          systHists120[bdt][syst]->SetLineColor(syst+2);
        }
        if (name.Contains(("th1f_sig_"+BDTtype[bdt]+"_"+productionTypes[0]+"_"+syst135mass[0]+"_"+syst135mass[syst]).c_str()) && !name.Contains("sigma")){
          systHists135[bdt][syst] = (TH1F*)inFile->Get(name.Data());
          systHists135[bdt][syst]->SetLineColor(syst+2);
        }
      }
    }
    // ---------------------------------------------
    std::string tName = temp->GetName();
    for (int i=0; i<nMasses; i++){
      int ind = tName.find(BDTmasses[i]+"_"+BDTmasses[i]);
      if (ind>0) tName.replace(ind,11,BDTmasses[i]);
    }
    temp->SetName(tName.c_str());
//.........这里部分代码省略.........

void plotSystFracs(TList* HistList, TH1F* compT, std::string name){

  gROOT->SetBatch();
  system("mkdir -p plots/ada/systs");
  system("mkdir -p plots/grad/systs");

  std::string bdt;
  TString str = compT->GetName();
  if (str.Contains("ada")) bdt="ada";
  if (str.Contains("grad")) bdt="grad";

  int nHists=HistList->GetEntries();
  TH1F *comp = linearBin(compT);
  
  gROOT->SetStyle("Plain");
  gROOT->ForceStyle();
  gStyle->SetOptStat(0);
  TCanvas *canv = new TCanvas("","",700,700);
  
  TLegend *leg = new TLegend(0.6,0.6,0.88,0.88);
  leg->SetLineColor(0);
  leg->SetFillColor(0);
  
  TF1 *line = new TF1("line","0.0",0.,comp->GetBinLowEdge(comp->GetNbinsX()+1));
  line->SetLineColor(kBlack);
  TF1 *line1 = new TF1("line","10.0",0.,comp->GetBinLowEdge(comp->GetNbinsX()+1));
  line1->SetLineColor(kGray+2);
  line1->SetLineStyle(2);
  TF1 *line2 = new TF1("line","-10.0",0.,comp->GetBinLowEdge(comp->GetNbinsX()+1));
  line2->SetLineColor(kGray+2);
  line2->SetLineStyle(2);
  
  int colors[10] = {kBlue,kMagenta,kGreen,kCyan,kRed,kBlue+3,kOrange+1,kSpring-1,kMagenta+3,kGreen+3};
  int color=0;
  for (int i=0; i<nHists; i++){
    TH1F *systHist = linearBin((TH1F*)HistList->At(i));
    systHist->Add(comp,-1);
    systHist->Divide(comp);
    systHist->Scale(100.);
    std::string systStr = systHist->GetName();
    int ind = systStr.rfind("cat0");
    std::string systName = systStr.substr(ind+5,systStr.size());
    systHist->SetLineColor(colors[color]);
    systHist->SetTitle(Form("%s",name.c_str()));
    systHist->GetYaxis()->SetTitle("Difference over nominal %");
    systHist->GetYaxis()->SetTitleOffset(1.4);
    systHist->GetXaxis()->SetTitle("BDT output");
    systHist->GetYaxis()->SetRangeUser(-100.,250);
    if (int(systName.find("Up"))>0){
      systHist->SetLineStyle(1);
      systName = systName.substr(0,systName.rfind("Up"));
      leg->AddEntry(systHist,systName.c_str(),"l");
      color++;
    }
    else if (int(systName.find("Down"))>0){
      systHist->SetLineStyle(2);
      systName = systName.substr(0,systName.rfind("Down"));
    }
    if (i==0) systHist->DrawCopy("hist");
    else systHist->DrawCopy("same hist");
  }
  leg->Draw("same");
  line->Draw("same");
  line1->Draw("same");
  line2->Draw("same");
  canv->Print(("plots/"+bdt+"/systs/"+name+".png").c_str(),"png");
  
  delete canv;
  
  systCalls++;
  
}

void changeAxis(bool armed, TString iExprMatch,float iMin, const char* filename, const char* pattern="", unsigned int debug=0)
{
  std::vector<std::string> histnames; histnames.clear();
  if( debug>0 ){
    std::cout << "file  = " << filename             << std::endl;
    std::cout << "old   = " << iExprMatch.Data()    << std::endl;    
    std::cout << "armed = " << armed                << std::endl;
  }
  TFile* file = new TFile(filename, "update");
  TIter nextDirectory(file->GetListOfKeys());
  std::vector<std::string> buffer;
  TKey* idir;
  while((idir = (TKey*)nextDirectory())){
    buffer.clear();
    if( idir->IsFolder() ){
      file->cd(); // make sure to start in directory head 
      if( debug>0 ){ std::cout << "Found directory: " << idir->GetName() << std::endl; }
      if( file->GetDirectory(idir->GetName()) ){
	file->cd(idir->GetName()); // change to sub-directory
	buffer = signalList(idir->GetName(), pattern,iExprMatch.Data(), debug);
      }
      // append to the vector of histograms to be rescaled
      for(std::vector<std::string>::const_iterator elem=buffer.begin(); elem!=buffer.end(); ++elem){
	histnames.push_back(*elem);
      }
      if(debug>1){
	std::cout << "added " << buffer.size() << " elements to histnames [" << histnames.size() << "] for directory " << idir->GetName() << std::endl;
      }
    }
  }
  // pick up files which are not kept in an extra folder
  file->cd(); buffer.clear();
  buffer = signalList("", pattern,iExprMatch.Data(), debug);
  // append to the vector of histograms to be rescaled
  for(std::vector<std::string>::const_iterator elem=buffer.begin(); elem!=buffer.end(); ++elem){
    histnames.push_back(*elem);
  }
  if(debug>1){
    std::cout << "added " << buffer.size() << " elements to histnames [" << histnames.size() << "] for file head" << std::endl;
  }

  for(std::vector<std::string>::const_iterator hist=histnames.begin(); hist!=histnames.end(); ++hist){
    file->cd();
    TH1F* h = (TH1F*)file->Get(hist->c_str());
    std::string histName;
    if(hist->find("/")!=std::string::npos){
      histName = hist->substr(hist->find("/")+1);
    }
    else{
      histName = *hist;
    }
    TH1F* hout = (TH1F*)h->Clone(histName.c_str());
    if(debug>1){
      std::cout << "...folder    : " << hist->substr(0, hist->find("/")).c_str() << std::endl;
      std::cout << "...histogram : " << hout->GetName () << " / " << hist->c_str() << std::endl; 
    }
    for(int i0 = 0; i0 < hout->GetNbinsX()+1; i0++) if(hout->GetXaxis()->GetBinCenter(i0) < iMin) hout->SetBinContent(i0,0);
    if(debug>1){ 
      std::cout << "...new name  : " << hout->GetName() << std::endl; 
      std::cout << "...new title : " << hout->GetTitle() << std::endl; 
    }
    if(armed){
      if(hist->find("/")!=std::string::npos){
	file->cd(hist->substr(0, hist->find("/")).c_str());
      }
      else{
	file->cd();
      }
      std::cout << "writing to file: " << hout->GetName() << " -- " << hist->substr(hist->find("/")+1).c_str() << std::endl;
      hout->Write(hist->substr(hist->find("/")+1).c_str(), TObject::kOverwrite); 
    }
  }
  file->Close();
  return;
}