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

void createInputs(int n = 2) 
{
   for(UInt_t i = 0; i < (UInt_t)n; ++i ) {
      TFile *file = TFile::Open(TString::Format("input%d.root",i),"RECREATE");
      TH1F * h = new TH1F("h1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
 
      Int_t nbins[5];
      Double_t xmin[5];
      Double_t xmax[5];
      for(UInt_t j = 0; j < 5; ++j) {
         nbins[j] = 10; xmin[j] = 0; xmax[j] = 10;
      }
      THnSparseF *sparse = new THnSparseF("sparse", "sparse", 5, nbins, xmin, xmax);
      Double_t coord[5] = {0.5, 1.5, 2.5, 3.5, 4.5};
      sparse->Fill(coord);
      sparse->Write();
      
      THStack *stack = new THStack("stack","");
      h = new TH1F("hs_1","",10,0,100);
      h->Fill(10.5); h->Fill(20.5);
      h->SetDirectory(0);
      stack->Add(h);
      h = new TH1F("hs_2","",10,0,100);
      h->Fill(30.5); h->Fill(40.5);
      h->SetDirectory(0);
      stack->Add(h);
      stack->Write();

      TGraph *gr = new TGraph(3);
      gr->SetName("exgraph");
      gr->SetPoint(0,1,1);
      gr->SetPoint(1,2,2);
      gr->SetPoint(2,3,3);
      
      gr->Write();
      
      TTree *tree = new TTree("tree","simplistic tree");
      Int_t data = 0;
      tree->Branch("data",&data);
      for(Int_t l = 0; l < 2; ++l) {
         data = l;
         tree->Fill();
      }
      
      file->Write();
      delete file;
   }
}

void PlotEnergySpecs(int num){

  cout<<"Run is "<<num<<endl;
  
  stringstream ss;
  
  TH1F * currentHisto;
  
  for (int i=0;i<40;i++){
    ss.str("");
    ss<<"Run"<<num<<"LA"<<i<<".root";
    TFile f(ss.str().c_str());
    
    currentHisto =(TH1F*)gDirectory->Get("ENoOverFlows1");
    if (i==0)
      currentHisto->Draw();
    else if ( i==38)
      currentHisto->Draw("same");

    currentHisto->SetLineColor(i+1);
    currentHisto->SetDirectory(0);
  }



}

TH1F* getHisto(string filename, string directoryname, string histoname) {
  TFile *file = new TFile(filename.c_str(),"READ");
  if (!file) {
    cout << "Could not open file " << filename << endl;
    return 0;
  }

  TDirectory *dir = (TDirectory*)file->FindObjectAny(directoryname.c_str());
  if (!dir) {
    cout << "Could not find directory " << directoryname << endl;
    delete file;
    return 0;
  }

  TH1F *hist = (TH1F*)dir->Get(histoname.c_str());
  if (!hist) {
    cout << "Could not find histogram " <<  histoname << " in directory " << directoryname << endl;
    delete dir;
    delete file;
    return 0;
  }

  hist->SetDirectory(0);
  delete dir;
  delete file;
  return hist;

}

void beffAnalysis(const char* input, const char* output) {
  // input
  TFile* inputFile = TFile::Open(input);
  TTree* btagEff = (TTree*)inputFile->Get("btagEff");

  // output
  TFile* outputFile = TFile::Open(output,"RECREATE");

  // histogram with proper binning... cloned later on 
  Double_t binning[23] = {20,25,30,35,40,45,50,55,60,65,70,75,80,85,90,95,100,120,140,160,180,200,1000};
  TH1F* ptSpectrum = new TH1F("PtSpectrum","PtSpectrum",22,binning);
  ptSpectrum->Sumw2();

  // produce the ratio plot for the 12 combinations of (CSVL,CSVM,CSVT),(Barrel,Endcap),(b,c,l)
  TClonesArray algorithms("TCut",3);
  new(algorithms[0]) TCut("CSVL","csv>0.244");
  new(algorithms[1]) TCut("CSVM","csv>0.679");
  new(algorithms[2]) TCut("CSVT","csv>0.898");

  TClonesArray etaRegions("TCut",2);
  new(etaRegions[0]) TCut("Barrel","abs(eta)<=1.2");
  new(etaRegions[1]) TCut("Endcaps","abs(eta)>1.2");

  TClonesArray flavor("TCut",3);
  new(flavor[0]) TCut("l","abs(flavor)!=4 && abs(flavor)!=5");
  new(flavor[1]) TCut("c","abs(flavor)==4");
  new(flavor[2]) TCut("b","abs(flavor)==5");

  for(int i=0; i< algorithms.GetEntries() ; ++i) {
    outputFile->mkdir(((TCut*)algorithms.At(i))->GetName());
    outputFile->cd(((TCut*)algorithms.At(i))->GetName());
    for(int j=0; j< etaRegions.GetEntries() ; ++j) {
      for(int k=0; k< flavor.GetEntries() ; ++k) {
        // histogram before tagging
        TH1F* pretag = ptSpectrum->Clone("pretag");
        btagEff->Draw("pt>>pretag",((*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // histogram after tagging
        TH1F* posttag = ptSpectrum->Clone("posttag");
        btagEff->Draw("pt>>posttag",((*(TCut*)algorithms.At(i))&&(*(TCut*)etaRegions.At(j))&&(*(TCut*)flavor.At(k)))*"eventWeight");
        // ratio
        TH1F* ratio = posttag->Clone(Form("h_eff_bTagOverGoodJet_pt%s_%s",((TCut*)flavor.At(k))->GetName(),
                                                                          ((TCut*)etaRegions.At(j))->GetName()));
        ratio->Divide(pretag);
        // cleanup
        delete pretag;
        delete posttag;
      }
    }
  }

  // cleanup
  algorithms.Delete();
  etaRegions.Delete();
  flavor.Delete();
  ptSpectrum->SetDirectory(0);
  outputFile->Write();
  outputFile->Close();
  inputFile->Close();
}

void Histogrammer::make_hist(const char* hname, const char* htitle, int nbins, double xlow, double xhigh, const char* xlabel, const char* ylabel){
	TH1F* h = new TH1F(hname, htitle, nbins, xlow, xhigh);
	h->GetXaxis()->SetTitle(xlabel);
	h->GetYaxis()->SetTitle(ylabel);
	h->SetDirectory(0);
	//h->Sumw2();
	hists[hname] = h;
}

  Hists(std::string name)
  {
    title            = name;
    h1McPt           = new TH1F(Form("h1McPt_%s",name.c_str()),";mcP_{T}(GeV/c)",40,1,5);
    h1HftMatchedMcPt = new TH1F(Form("h1HftMatchedMcPt_%s",name.c_str()),";mcP_{T}(GeV/c)",40,1,5);
    h1TpcMcPt        = new TH1F(Form("h1TpcMcPt_%s",name.c_str()),";mcP_{T}(GeV/c)",40,1,5);
    gTpcEff          = new TGraphAsymmErrors; gTpcEff->SetName(Form("gTpcEff_%s",name.c_str()));
    gHftEff          = new TGraphAsymmErrors; gHftEff->SetName(Form("gHftEff_%s",name.c_str()));
    h1DcaZPosEta     = new TH1F(Form("h1DcaZPosEta_%s",name.c_str()),";dcaZ(cm)",400,-1,1);
    h1DcaZNegEta     = new TH1F(Form("h1DcaZNegEta_%s",name.c_str()),";dcaZ(cm)",400,-1,1);
    h1DcaZ           = new TH1F(Form("h1DcaZ_%s",name.c_str()),";dcaZ(cm)",400,-1,1);

    h1McPt->SetDirectory(0);
    h1HftMatchedMcPt->SetDirectory(0);
    h1TpcMcPt->SetDirectory(0);
    h1DcaZ->SetDirectory(0);
  }

TH1F* LikelihoodSpace::get_projection(std::string name) {
  int default_nbins = 100;
  gEnv->GetValue("Hist.Binning.1D.x", default_nbins);
  gEnv->SetValue("Hist.Binning.1D.x", 10000);
  this->samples->Draw((name + ">>_hp").c_str(), "", "goff");
  gEnv->SetValue("Hist.Binning.1D.x", default_nbins);
  TH1F* hp = dynamic_cast<TH1F*>(gDirectory->FindObject("_hp"));
  assert(hp);
  hp->SetDirectory(NULL);

  return hp;
}

TH1F* loadHistogram1D(
    const std::vector<std::string>& fileList, 
    const std::string& histName, 
    const std::string& sys, 
    double scale,
    int rebin
    )
{
    TH1F* result = nullptr;
    for (const std::string& file: fileList)
    {
        TFile input(file.c_str(),"r");
        
        TH1F* hist=nullptr;
        if ((sys=="" or sys=="nominal") and input.FindKey(histName.c_str()))
        {
            hist = static_cast<TH1F*>(input.Get(histName.c_str()));
        }
        else if (input.FindKey((histName+"__"+sys).c_str()))
        {
            hist = static_cast<TH1F*>(input.Get((histName+"__"+sys).c_str()));
        }
        else if (input.FindKey(histName.c_str()))
        {
            log(WARNING,"using fallback '%s' for '%s'\n",histName.c_str(),(histName+"__"+sys).c_str());
            hist = static_cast<TH1F*>(input.Get(histName.c_str()));
        }
        else
        {
            log(ERROR,"neither '%s' nor '%s' found in file '%s'\n",histName.c_str(),(histName+"__"+sys).c_str(),file.c_str());
            continue;
        }
        hist->Scale(scale);
        hist->Rebin(hist->GetNbinsX()/rebin);
        if (!result)
        {
            result = new TH1F(*hist);
            result->Sumw2();
            result->SetDirectory(0);
        }
        else
        {
            result->Add(hist);
        }
    }
    if (!result)
    {
        log(CRITICAL,"Error while finding histogram '"+histName+"' with sys '"+sys+"' in given files");
        throw std::string("Error while finding histogram '"+histName+"' with sys '"+sys+"' in given files");
    }
    return result;
}

//*************************************************************************************************
//Function to create a distribution histogram for the signal for the variable given
//*************************************************************************************************
TH1F* makeSignalDistributionHistogram ( int variableIndex, MitNtupleEvent* event , Int_t nentries, TLegend *legend, int OnlyThisFinalState = -1) {
  assert(variableIndex >= 0 && variableIndex < NVARIABLES);
  string label = "NSignalEventsPass";
  string histName = "h" + fVariableNames[variableIndex] + "_" + label;
  string axisLabel = ";" + fVariableNames[variableIndex] + ";Number of Events";
  TH1F *histogram = new TH1F(histName.c_str(),axisLabel.c_str(),NBINS, 
                             fVariableRangeMin[variableIndex], fVariableRangeMax[variableIndex]);
  histogram->SetDirectory(0);

  for (int n=0;n<nentries;n++) {
    event->GetEntry(n);
    float eventweight = event->H_weight;
    int finalstatetype = (int)event->H_ltype;

    //only look at events with the given finalstatetype
    if (OnlyThisFinalState >= 0 && OnlyThisFinalState != finalstatetype)
      continue;

    Float_t cutAboveValues[NVARIABLES];
    Float_t cutBelowValues[NVARIABLES];
    Float_t variableValues[NVARIABLES];   

    //only look at events with the given finalstatetype
    //convert finalstatetype into the array index
    int fs;
    if (finalstatetype == 10)
      fs = 0;
    else if (finalstatetype == 11)
      fs = 1;
    else if (finalstatetype == 12)
      fs = 2;
    else {
      fs = -1;
      continue;
    }

    for (int j=0;j<NVARIABLES;j++) {
      cutAboveValues[j] = fCutAboveInitialValue[fs][j];
      cutBelowValues[j] = fCutBelowInitialValue[fs][j];
      variableValues[j] = GetVariableValue(j, event);      
    }

    //for N-1 distribution
    //if (passCut(event,variableValues, cutAboveValues, cutBelowValues, variableIndex)) { 
    if (passCut(event,variableValues, cutAboveValues, cutBelowValues)) {
      histogram->Fill(variableValues[variableIndex], eventweight);
    }  
  }
  legend->AddEntry(histogram, "signal", "LP"); 
  return histogram;
}

int ScanChain( TChain* chain, int nEvents = -1, std::string skimFilePrefix="") {

  TObjArray *listOfFiles = chain->GetListOfFiles();

  unsigned int nEventsChain=0;
  if(nEvents==-1) 
    nEvents = chain->GetEntries();
  nEventsChain = nEvents;
  unsigned int nEventsTotal = 0;
  TDirectory *rootdir = gDirectory->GetDirectory("Rint:");

  TH1F *samplehisto = new TH1F("samplehisto", "Example histogram", 200,0,200);
  samplehisto->SetDirectory(rootdir);
  // file loop
  TIter fileIter(listOfFiles);
  TFile *currentFile = 0;
  while ( currentFile = (TFile*)fileIter.Next() ) {
    TFile f(currentFile->GetTitle());
    TTree *tree = (TTree*)f.Get("Events");
    cms2.Init(tree);
    
    //Event Loop
    unsigned int nEvents = tree->GetEntries();
    for( unsigned int event = 0; event < nEvents; ++event) {
      cms2.GetEntry(event);
      ++nEventsTotal;
      // Progress feedback to the user
      if(nEventsTotal%1000 == 0) {
        // xterm magic from L. Vacavant and A. Cerri
        if (isatty(1)) {
          printf("\015\033[32m ---> \033[1m\033[31m%4.1f%%"
          "\033[0m\033[32m <---\033[0m\015", (float)nEventsTotal/(nEventsChain*0.01));
          fflush(stdout);
        }
      }//if(nEventsTotal%20000 == 0) {


    }
    delete tree;
    f.Close();
  }

  if ( nEventsChain != nEventsTotal ) {
    std::cout << "ERROR: number of events from files is not equal to total number of events" << std::endl;
  }

  samplehisto->Draw();
  return 0;
}

void make(TDirectory & out, TObject * o) {
  TDirectory * dir;
  TH1F * th1f;
  TH1D * th1d;
  TH2F * th2f;
  TH2D * th2d;
  out.cd();
  if((dir = dynamic_cast<TDirectory*>(o)) != 0) {
    TDirectory * outDir = out.mkdir(dir->GetName(), dir->GetTitle());
    TIter next(dir->GetListOfKeys());
    TKey *key;
    while( (key = dynamic_cast<TKey*>(next())) ) {
      string className(key->GetClassName());
      string name(key->GetName());
      TObject * obj = dir->Get(name.c_str());
      if(obj == 0) {
	cerr <<"error: key " << name << " not found in directory " << dir->GetName() << endl;
	exit(-1);
      }
      make(*outDir, obj);
    }
  } else if((th1f = dynamic_cast<TH1F*>(o)) != 0) {
    TH1F *h = (TH1F*) th1f->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th1d = dynamic_cast<TH1D*>(o)) != 0) {
    TH1D *h = (TH1D*) th1d->Clone();
    h->Reset();
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2f = dynamic_cast<TH2F*>(o)) != 0) {
    TH2F *h = (TH2F*) th2f->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  } else if((th2d = dynamic_cast<TH2D*>(o)) != 0) {
    TH2D *h = (TH2D*) th2d->Clone();
    h->Reset();   
    h->Sumw2();
    h->SetDirectory(&out);
  }
}

PUReweight::PUReweight(int nFiles, char** fileNames, std::string PUfilename){
	PUweightSum = 0.0;
	events = 0;
	TFile* pileupFile = new TFile(PUfilename.c_str(),"READ");
	PUweightHist = (TH1D*)pileupFile->Get("pileup");
	PUweightHist->SetDirectory(0);
	pileupFile->Close();
	TH1D* PUbackup;
	if(PUweightHist->GetNbinsX() != 100){
		std::cout << "Wrong number of bins in the pileup histogram" << std::endl;
		PUbackup = new TH1D("pileup_new","pileup_new",100,0,100);
		for(int ibin=1; ibin <= PUweightHist->GetNbinsX(); ibin++){
			PUbackup->SetBinContent(ibin, PUweightHist->GetBinContent(ibin));
			// assuming the same scale
		}
		PUweightHist = PUbackup;
	}
	double PUweightInt = PUweightHist->Integral();
	TH1F* mcPU = NULL;
	for(int nmcfile = 0; nmcfile<nFiles; nmcfile++){
		std::cout << "reading file " << std::string(fileNames[nmcfile]) << std::endl;
		TFile* mcFile = new TFile(fileNames[nmcfile],"READ");
		if(!(mcFile->Get("ggNtuplizer/hPU"))) {
			std::cout << "no hPU histogram here!" << std::endl;
			delete PUweightHist;
			PUweightHist = NULL;
			return;
		}
		if( mcPU==NULL) mcPU = (TH1F*)mcFile->Get("ggNtuplizer/hPU");
		else mcPU->Add((TH1F*)mcFile->Get("ggNtuplizer/hPU"));
		mcPU->SetDirectory(0);
		mcFile->Close();
	}
	mcPU->Scale(1.0/mcPU->Integral());
	PUweightHist->Divide(mcPU);
	PUweightHist->Scale(1.0/PUweightInt);
	delete mcPU;
}

void SignfificanceT2tt(bool pval, int exp=0){


    TFile *f = TFile::Open("Significances2DHistograms_T2tt.root");
    TH2F *h;
    if(pval){
        if(exp==0) h = (TH2F*)f->Get("hpObs");
        else if(exp==1) h = (TH2F*)f->Get("hpExpPosteriori");
        else if(exp==2) h = (TH2F*)f->Get("hpExpPriori");
    }
    else {
        if(exp==0) h = (TH2F*)f->Get("hObs");
        else if(exp==1) h = (TH2F*)f->Get("hExpPosteriori");
        else if(exp==2) h = (TH2F*)f->Get("hExpPriori");
    }
    h->GetXaxis()->SetTitle("m_{#tilde{t}} [GeV]");
    h->GetXaxis()->SetLabelFont(42);
    h->GetXaxis()->SetLabelSize(0.035);
    h->GetXaxis()->SetTitleSize(0.05);
    h->GetXaxis()->SetTitleOffset(1.2);
    h->GetXaxis()->SetTitleFont(42);
    h->GetYaxis()->SetTitle("m_{#tilde{#chi}_{1}^{0}} [GeV]");
    h->GetYaxis()->SetLabelFont(42);
    h->GetYaxis()->SetLabelSize(0.035);
    h->GetYaxis()->SetTitleSize(0.05);
    h->GetYaxis()->SetTitleOffset(1.35);
    h->GetYaxis()->SetTitleFont(42);
    double maximum = h->GetMaximum();
    double minimum = h->GetMinimum();
    double sigmin = 99; int sigminx=-1; int sigminy=-1; if(pval) sigmin = -99;
    h->GetZaxis()->SetRangeUser(minimum,maximum);
    for(int x = 1; x<=h->GetNbinsX();++x){
        for(int y = 1; y<=h->GetNbinsX();++y){
            if(!pval&&h->GetBinContent(x,y)<sigmin){ sigmin =h->GetBinContent(x,y); sigminx = x; sigminy = y; }
            if( pval&&h->GetBinContent(x,y)>sigmin){ sigmin =h->GetBinContent(x,y); sigminx = x; sigminy = y; }
            if(!pval&&h->GetXaxis()->GetBinLowEdge(x)<h->GetYaxis()->GetBinLowEdge(y)+75) h->SetBinContent(x,y,-999);
            if(h->GetXaxis()->GetBinLowEdge(x)<374) continue;
            if(h->GetXaxis()->GetBinLowEdge(x)>424) continue;
            if(h->GetYaxis()->GetBinLowEdge(y)<199) continue;
            if(h->GetYaxis()->GetBinLowEdge(y)>249) continue;
            if(!pval&&h->GetBinContent(x,y)>0.3) h->SetBinContent(x,y,0.05);
        }
    }
    h->GetZaxis()->SetRangeUser(minimum,maximum);
    if(!pval) cout << "minimal significance " << sigmin << " at " << h->GetXaxis()->GetBinLowEdge(sigminx) << "-" << h->GetYaxis()->GetBinLowEdge(sigminy) << endl;
    else cout << "maximal p- value " << sigmin << " at " << h->GetXaxis()->GetBinLowEdge(sigminx) << "-" << h->GetYaxis()->GetBinLowEdge(sigminy) << endl;
    
   TCanvas *c1 = new TCanvas("c1", "c1",50,50,600,600);
   gStyle->SetOptFit(1);
   gStyle->SetOptStat(0);
   gStyle->SetOptTitle(0);
   gStyle->SetErrorX(0.5); 
   //c1->Range(-6.311689,-1.891383,28.75325,4.56342);
   c1->SetFillColor(0);
   c1->SetBorderMode(0);
   c1->SetBorderSize(2);
   //c1->SetLogy();
   c1->SetTickx(1);
   c1->SetTicky(1);
   c1->SetLeftMargin(0.15);
//   c1->SetRightMargin(0.05);
    c1->SetRightMargin(0.15);
   c1->SetTopMargin(0.07);
   c1->SetBottomMargin(0.15);
   c1->SetFrameFillStyle(0);
   c1->SetFrameBorderMode(0);
   c1->SetFrameFillStyle(0);
   c1->SetFrameBorderMode(0);
   gStyle->SetHatchesLineWidth(0);

   TH1F *hSum = new TH1F("hSum","",10,100,900);
   hSum->SetMinimum(0.);
   hSum->SetMaximum(550);
   hSum->SetDirectory(0);
   hSum->SetStats(0);
    hSum->Draw();
    hSum->GetYaxis()->SetRangeUser(0,500);
    hSum->GetXaxis()->SetRangeUser(100,900);

   Int_t ci;   // for color index setting
   ci = TColor::GetColor("#000099");
   hSum->SetLineColor(ci);
   hSum->SetLineStyle(0);
   hSum->SetMarkerStyle(20);
   hSum->GetXaxis()->SetTitle("m_{#tilde{t}} [GeV]");
   //hSum->GetXaxis()->SetBit(TAxis::kLabelsVert);
   hSum->GetXaxis()->SetLabelFont(42);
   //hSum->GetXaxis()->SetLabelOffset(0.005);
   hSum->GetXaxis()->SetLabelSize(0.035);
   hSum->GetXaxis()->SetTitleSize(0.06);
   hSum->GetXaxis()->SetTitleOffset(1.2);
   hSum->GetXaxis()->SetTitleFont(42);
   hSum->GetYaxis()->SetTitle("m_{#tilde{#chi}}_{1}^{0} [GeV]");
   hSum->GetYaxis()->SetLabelFont(42);
   //hSum->GetYaxis()->SetLabelOffset(0.007);
   hSum->GetYaxis()->SetLabelSize(0.035);
   hSum->GetYaxis()->SetTitleSize(0.05);
   hSum->GetYaxis()->SetTitleOffset(1.3);
   hSum->GetYaxis()->SetTitleFont(42);
    
//.........这里部分代码省略.........

TObjArray getDistributionFromPlotter(TString plot,TString baseURL="~/scratch0/top-nosyst/plotter.root")
{
  using namespace std;

  setStyle();
  
  //the relevant processes
  TString procs[]={"Di-bosons","Single top", "W+jets", "Z-#gamma^{*}+jets#rightarrow ll", "other t#bar{t}", "t#bar{t} dileptons", "data"};
  const size_t nprocs=sizeof(procs)/sizeof(TString);

  //containers for the histos
  TList *data     = new TList;
  TH1F *totalData=0;
  TList *mc       = new TList;
  TH1F *totalMC=0;
  TH1F *ttbarMC=0;
  TList *spimpose = new TList;
  
  //open file with histograms
  TFile *f=TFile::Open(baseURL);
  for(size_t iproc=0; iproc<nprocs; iproc++)
    {
      TH1F *histo = (TH1F *) f->Get(procs[iproc]+"/"+plot);
      if(histo==0) { cout << "[Warning] " << plot << " not found for " << procs[iproc] << " ...skipping" << endl; continue; }
      histo->SetDirectory(0);
      if(procs[iproc].Contains("data")) 
	{
	  data->Add(histo);
	  if(totalData==0) { totalData = (TH1F *) histo->Clone(plot+"totaldata"); totalData->SetDirectory(0); }
	  else             { totalData->Add(histo); }
	}
      else
	{
	  mc->Add(histo);
	  if(totalMC==0) { totalMC = (TH1F *) histo->Clone(plot+"totalmc"); totalMC->SetDirectory(0); }
	  else           { totalMC->Add(histo); }
	  if(procs[iproc]=="t#bar{t} dileptons") { ttbarMC = (TH1F *)histo->Clone(plot+"ttbar"); ttbarMC->SetDirectory(0); }
	}
    }
  f->Close();

  cout << "Drawing now" << endl;
  
  //draw
  TCanvas *cnv = getNewCanvas(plot+"c",plot+"c",false);
  cnv->Clear();
  cnv->SetWindowSize(600,600);
  cnv->cd(); 
  TLegend *leg=showPlotsAndMCtoDataComparison(cnv,*mc,*spimpose,*data,false);
  formatForCmsPublic(cnv,leg,"CMS preliminary", 4);
  cnv->SaveAs(plot+".C");

  cout << "Stat comparison for " << plot << endl
       << "Sample \t Average \t\t RMS " << endl
       << "----------------------------------------" << endl;
  if(totalData) cout << "Data \t " << totalData->GetMean() << " +/- " << totalData->GetMeanError() << " \t " << totalData->GetRMS() << " +/- " << totalData->GetRMSError() << endl;
  if(totalMC)   cout << "MC \t "   << totalMC->GetMean() << " +/- " << totalMC->GetMeanError() << " \t " << totalMC->GetRMS() << " +/- " << totalMC->GetRMSError() << endl;
  if(ttbarMC)   cout << "Signal \t "   << ttbarMC->GetMean() << " +/- " << ttbarMC->GetMeanError() << " \t " << ttbarMC->GetRMS() << " +/- " << ttbarMC->GetRMSError() << endl;

  TObjArray res;
  res.Add(totalData);
  res.Add(totalMC);
  res.Add(ttbarMC);
  res.Add(mc);
  return res;
}

void compare_wls1(TString filename="../p15m_nwls/wcsim.root",TString histoname="p15m_nwls", Int_t *flag, TString rootfile = "temp.root") {

    TFile *file1;
    if (*flag!=0) {
        //file1 = new TFile(rootfile,"Update");
        file1 = new TFile(rootfile,"RECREATE");
    } else {
        file1 = new TFile(rootfile,"RECREATE");
    }
    TString filename1;
    filename1 = histoname + "_digi";
    TTree *T = new TTree(filename1,filename1);
    T->SetDirectory(file1);
    Double_t diginpe,digitime,cor_digitime,digitheta,dis_digihit;
    Int_t neve;

    Double_t mom;
    Double_t pos_x,pos_y,pos_z;
    Double_t dir_x,dir_y,dir_z;
    Double_t tube_x,tube_y,tube_z;
    Double_t totankdis;
    Double_t vertex[3],dir[3];

    Double_t tube_pos[3];

    T->Branch("digi_eve",&neve,"data/I");
    T->Branch("diginpe",&diginpe,"data/D");
    T->Branch("digitime",&digitime,"data/D");
    T->Branch("cor_digitime",&cor_digitime,"data/D");
    T->Branch("digitheta",&digitheta,"data/D");
    T->Branch("dis_dighit",&dis_digihit,"data/D");

    T->Branch("mom",&mom,"data/D");
    T->Branch("totankdis",&totankdis,"data/D");
    T->Branch("pos_x",&vertex[0],"data/D");
    T->Branch("pos_y",&vertex[1],"data/D");
    T->Branch("pos_z",&vertex[2],"data/D");

    T->Branch("dir_x",&dir[0],"data/D");
    T->Branch("dir_y",&dir[1],"data/D");
    T->Branch("dir_z",&dir[2],"data/D");

    T->Branch("tube_x",&tube_pos[0],"data/D");
    T->Branch("tube_y",&tube_pos[1],"data/D");
    T->Branch("tube_z",&tube_pos[2],"data/D");

    filename1 = histoname + "_hit";
    TTree *t1 = new TTree(filename1,filename1);
    t1->SetDirectory(file1);

    Double_t wavelength, truetime, corr_time,theta,distance,index;
    Int_t qe_flag,parentid,tubeid,totalpe;

    Int_t ntracks;

    t1->Branch("ntracks",&ntracks,"data/I");
    t1->Branch("neve",&neve,"data/I");
    t1->Branch("wavelength",&wavelength,"data/D");
    t1->Branch("truetime",&truetime,"data/D");
    t1->Branch("corr_time",&corr_time,"data/D");
    t1->Branch("theta",&theta,"data/D");
    t1->Branch("distance",&distance,"data/D");
    t1->Branch("index",&index,"data/D");

    t1->Branch("mom",&mom,"data/D");
    t1->Branch("totankdis",&totankdis,"data/D");
    t1->Branch("pos_x",&vertex[0],"data/D");
    t1->Branch("pos_y",&vertex[1],"data/D");
    t1->Branch("pos_z",&vertex[2],"data/D");

    t1->Branch("dir_x",&dir[0],"data/D");
    t1->Branch("dir_y",&dir[1],"data/D");
    t1->Branch("dir_z",&dir[2],"data/D");

    t1->Branch("tube_x",&tube_pos[0],"data/D");
    t1->Branch("tube_y",&tube_pos[1],"data/D");
    t1->Branch("tube_z",&tube_pos[2],"data/D");

    // t1->Branch("pos_x",&pos_x,"data/D");
//   t1->Branch("pos_y",&pos_y,"data/D");
//   t1->Branch("pos_z",&pos_z,"data/D");

//   t1->Branch("dir_x",&dir_x,"data/D");
//   t1->Branch("dir_y",&dir_y,"data/D");
//   t1->Branch("dir_z",&dir_z,"data/D");

//   t1->Branch("tube_x",&tube_x,"data/D");
//   t1->Branch("tube_y",&tube_y,"data/D");
//   t1->Branch("tube_z",&tube_z,"data/D");

    t1->Branch("qe_flag",&qe_flag,"data/I");
    t1->Branch("parentid",&parentid,"data/I");
    t1->Branch("tubeid",&tubeid,"data/I");
    t1->Branch("totalpe",&totalpe,"data/I");


    TFile *file = new TFile(filename);
    TTree  *wcsimT = file->Get("wcsimT");
    WCSimRootEvent *wcsimrootsuperevent = new WCSimRootEvent();
    wcsimT->SetBranchAddress("wcsimrootevent",&wcsimrootsuperevent);
//.........这里部分代码省略.........