C++ TH1::GetTitle方法代码示例

  void getPlotData() {
    TH1 * h = (TH1*) m_file->Get(m_direc.c_str());

    for (int i=0; i<h->GetXaxis()->GetNbins(); i++) {
      m_xs.push_back(h->GetXaxis()->GetBinCenter(i));
      m_ys.push_back(h->GetBinContent(i));
    }

    m_plot->m_xAxisTitle = std::string(h->GetXaxis()->GetTitle());
    m_plot->m_yAxisTitle = std::string(h->GetYaxis()->GetTitle());
    m_plot->m_title = std::string(h->GetTitle());

    std::stringstream ssN, ssMu, ssSig, ssUF, ssOF;
    ssN << std::setprecision(4) << h->GetEntries();
    ssMu << std::setprecision(4) << h->GetMean();
    ssSig << std::setprecision(4) << h->GetRMS();
    ssUF << std::setprecision(4) << h->GetBinContent(0);
    ssOF << std::setprecision(4) << h->GetBinContent(h->GetNbinsX() + 1);

    m_statsTitles.push_back("N:");
    m_statsTitles.push_back("mu:");
    m_statsTitles.push_back("sig:");
    m_statsTitles.push_back("UF:");
    m_statsTitles.push_back("OF:");

    m_statsValues.push_back(ssN.str());
    m_statsValues.push_back(ssMu.str());
    m_statsValues.push_back(ssSig.str());
    m_statsValues.push_back(ssUF.str());
    m_statsValues.push_back(ssOF.str());
  }

void counts(int run, int lumistart, int lumiend, string type, map<string,vector<tripletI> > &cnt, vector<tripletD> &cntref, string hlttype, bool docnt, bool doref) {
   TString filename = basedir + Form("/DQM_V0001_HLTpb_R000%i.root",run);
   TFile *f = new TFile(filename);
   if (!f->IsOpen()) {
      cout << "Error, could not open " << filename << endl;
      return;
   }
   TString tdirname = Form("DQMData/Run %i/HLT/Run summary/TriggerRates/",run) + TString(type);
   f->cd(tdirname);
   TProfile *hlumi = (TProfile*) f->Get(Form("DQMData/Run %i/HLT/Run summary/LumiMonitoring/lumiVsLS",run));
   if (extrapol) extrapolate(hlumi);

   // if HLT: accept, error, pass L1 seed, pass prescaler, reject

   TIter next(gDirectory->GetListOfKeys());
   TKey *key;
   while ((key = (TKey*)next())) {
      TClass *cl = gROOT->GetClass(key->GetClassName());
      
      // it must be an histogram
      if (!cl->InheritsFrom("TH1")) continue;
      TH1 *h = (TH1*)key->ReadObj();

      // the name must match one of the requested patterns
      bool match=false; TString hname(h->GetName());
      for (vector<TRegexp>::const_iterator it=patterns.begin(); it!=patterns.end(); it++) {
         if (hname(*it).Length()!=0) {
            match=true;
            break;
         }
      }
      if (!match) continue;

      int nlumis = (lumiend+1-lumistart);

      if (extrapol) extrapolate(h);

      if (type != "HLT") fill(cnt[h->GetName()], h, run, lumistart, lumiend, docnt);
      else {
         string htitle(h->GetTitle());
         if (htitle.find(hlttype) == string::npos) continue;
         else {
            TString thepath; Ssiz_t from=0; TString(htitle).Tokenize(thepath,from," ");
            fill(cnt[thepath.Data()], h, run, lumistart, lumiend, docnt);
         }
      }
   }

   if (doref) {
      fill(cntref, hlumi, run, lumistart, lumiend);
   }

   f->Close();
   delete f;
}

  /** 
   * Create ratios to other data 
   * 
   * @param ib      Bin number  
   * @param res     Result
   * @param alice   ALICE result if any
   * @param cms     CMS result if any
   * @param all     Stack to add ratio to 
   */
  void Ratio2Stack(Int_t ib, TH1* res, TGraph* alice, TGraph* cms, THStack* all)
  {
    if (!all || !res || !(alice || cms)) return;

    Int_t        off  = 5*ib;
    TGraph*      gs[] = { (alice ? alice : cms), (alice ? cms : 0), 0 };
    TGraph**     pg   = gs;
    while (*pg) { 
      TGraph*     g = *pg;
      const char* n = (g == alice ? "ALICE" : "CMS");

      TH1*    r = static_cast<TH1*>(res->Clone(Form("ratio%s", n)));
      TString tit(r->GetTitle());
      tit.ReplaceAll("Corrected", Form("Ratio to %s", n));
      r->SetTitle(tit);
      r->SetMarkerColor(g->GetMarkerColor());
      r->SetLineColor(g->GetLineColor());

      TObject* tst = r->FindObject("legend");
      if (tst) r->GetListOfFunctions()->Remove(tst);

      for (Int_t i = 1; i <= r->GetNbinsX(); i++) {
	Double_t c = r->GetBinContent(i);
	Double_t e = r->GetBinError(i);
	Double_t o = g->Eval(r->GetBinCenter(i));
	if (o < 1e-12) { 
	  r->SetBinContent(i, 0);
	  r->SetBinError(i, 0);
	  continue;
	}
	r->SetBinContent(i, (c - o) / o + off);
	r->SetBinError(i, e / o);
      }
      all->Add(r);
      pg++;
    }
    TLegend* leg = StackLegend(all);
    if (!leg) return;
      
    TString   txt      = res->GetTitle();
    txt.ReplaceAll("Corrected P(#it{N}_{ch}) in ", "");
    if      (ib == 0) txt.Append(" "); // (#times1)");
    // else if (ib == 1) txt.Append(" (#times10)");
    else              txt.Append(Form(" (+%d)", off));

    TObject* dummy = 0;
    TLegendEntry* e = leg->AddEntry(dummy, txt, "p");
    e->SetMarkerStyle(res->GetMarkerStyle());
    e->SetMarkerSize(res->GetMarkerSize());
    e->SetMarkerColor(kBlack);
    e->SetFillColor(0);
    e->SetFillStyle(0);
    e->SetLineColor(kBlack);
  }

RooHistN::RooHistN(const TH1 &data1, const TH1 &data2, Double_t nominalBinWidth, Double_t nSigma, Double_t xErrorFrac) :
  TGraphAsymmErrors(), _nominalBinWidth(nominalBinWidth), _nSigma(nSigma), _rawEntries(-1)
{
  // Create a histogram from the asymmetry between the specified TH1 objects
  // which may have fixed or variable bin widths, but which must both have
  // the same binning. The asymmetry is calculated as (1-2)/(1+2). Error bars are
  // calculated using Binomial statistics. Prints a warning and rounds
  // any bins with non-integer contents. Use the optional parameter to
  // specify the confidence level in units of sigma to use for
  // calculating error bars. The nominal bin width specifies the
  // default used by addAsymmetryBin(), and is used to set the relative
  // normalization of bins with different widths. If not set, the
  // nominal bin width is calculated as range/nbins.

  initialize();
  // copy the first input histogram's name and title
  SetName(data1.GetName());
  SetTitle(data1.GetTitle());
  // calculate our nominal bin width if necessary
  if(_nominalBinWidth == 0) {
    const TAxis *axis= ((TH1&)data1).GetXaxis();
    if(axis->GetNbins() > 0) _nominalBinWidth= (axis->GetXmax() - axis->GetXmin())/axis->GetNbins();
  }
  setYAxisLabel(Form("Asymmetry (%s - %s)/(%s + %s)",
		     data1.GetName(),data2.GetName(),data1.GetName(),data2.GetName()));
  // initialize our contents from the input histogram contents
  Int_t nbin= data1.GetNbinsX();
  if(data2.GetNbinsX() != nbin) {
    coutE(InputArguments) << "RooHistN::RooHistN: histograms have different number of bins" << endl;
    return;
  }
  for(Int_t bin= 1; bin <= nbin; bin++) {
    Axis_t x= data1.GetBinCenter(bin);
    if(fabs(data2.GetBinCenter(bin)-x)>1e-10) {
      coutW(InputArguments) << "RooHistN::RooHistN: histograms have different centers for bin " << bin << endl;
    }
    Stat_t y1= data1.GetBinContent(bin);
    Stat_t y2= data2.GetBinContent(bin);
    addAsymmetryBin(x,roundBin(y1),roundBin(y2),data1.GetBinWidth(bin),xErrorFrac);
  }
  // we do not have a meaningful number of entries
  _entries= -1;
}

TString* get_var_names( Int_t nVars )
{
   const TString directories[6] = { "InputVariables_NoTransform",
                                    "InputVariables_DecorrTransform",
                                    "InputVariables_PCATransform",
				    "InputVariables_Id",
				    "InputVariables_Norm",
				    "InputVariables_Deco"};

   TDirectory* dir = 0;
   for (Int_t i=0; i<6; i++) {
      dir = (TDirectory*)Network_GFile->Get( directories[i] );
      if (dir != 0) break;
   }
   if (dir==0) {
      cout << "*** Big troubles in macro \"network.C\": could not find directory for input variables, "
           << "and hence could not determine variable names --> abort" << endl;
      return 0;
   }
   cout << "--> go into directory: " << dir->GetName() << endl;
   dir->cd();

   TString* vars = new TString[nVars];
   Int_t ivar = 0;

   // loop over all objects in directory
   TIter next(dir->GetListOfKeys());
   TKey* key = 0;
   while ((key = (TKey*)next())) {
      if (key->GetCycle() != 1) continue;

      if(!TString(key->GetName()).Contains("__S") &&
	 !TString(key->GetName()).Contains("__r")) continue;

      // make sure, that we only look at histograms
      TClass *cl = gROOT->GetClass(key->GetClassName());
      if (!cl->InheritsFrom("TH1")) continue;
      TH1 *sig = (TH1*)key->ReadObj();
      TString hname = sig->GetTitle();
      
      vars[ivar] = hname; ivar++;

      if (ivar > nVars-1) break;
   }      
   
   if (ivar != nVars-1) { // bias layer is also in nVars counts
      cout << "*** Troubles in \"network.C\": did not reproduce correct number of "
           << "input variables: " << ivar << " != " << nVars << endl;
   }

   return vars;

   // ------------- old way (not good) -------------

   //    TString fname = "weights/TMVAnalysis_MLP.weights.txt";
   //    ifstream fin( fname );
   //    if (!fin.good( )) { // file not found --> Error
   //       cout << "Error opening " << fname << endl;
   //       exit(1);
   //    }
   
   //    Int_t   idummy;
   //    Float_t fdummy;
   //    TString dummy = "";
   
   //    // file header with name
   //    while (!dummy.Contains("#VAR")) fin >> dummy;
   //    fin >> dummy >> dummy >> dummy; // the rest of header line
   
   //    // number of variables
   //    fin >> dummy >> idummy;
   //    // at this point, we should have idummy == nVars
   
   //    // variable mins and maxes
   //    TString* vars = new TString[nVars];
   //    for (Int_t i = 0; i < idummy; i++) fin >> vars[i] >> dummy >> dummy >> dummy;
   
   //    fin.close();
   
   //    return vars;
}

//.........这里部分代码省略.........
      TString hname(sig->GetName());

      //normalize to 1
      NormalizeHist(sig);      

      // create new canvas
      if (countPad%noPadPerCanv==0) {
         ++countCanvas;
         canv = new TCanvas( Form("canvas%d", countCanvas), title,
                             countCanvas*50+50, countCanvas*20, width, height );
         canv->Divide(xPad,yPad);
         canv->SetFillColor(kWhite);
         canv->Draw();
      }

      TPad* cPad = (TPad*)canv->cd(countPad++%noPadPerCanv+1);
      cPad->SetFillColor(kWhite);

      // find the corredponding backgrouns histo
      TString bgname = hname;
      bgname.ReplaceAll("__Signal","__Background");
      TH1 *bgd = (TH1*)dir->Get(bgname);
      if (bgd == NULL) {
         cout << "ERROR!!! couldn't find background histo for" << hname << endl;
         exit;
      }
      //normalize to 1
      NormalizeHist(bgd);


      // this is set but not stored during plot creation in MVA_Factory
      TMVAGlob::SetSignalAndBackgroundStyle( sig, (isRegression ? 0 : bgd) );            

      sig->SetTitle( TString( htitle ) + ": " + sig->GetTitle() );
      TMVAGlob::SetFrameStyle( sig, 1.2 );

      // normalise both signal and background
//       if (!isRegression) TMVAGlob::NormalizeHists( sig, bgd );
//       else {
//          // change histogram title for target
//          TString nme = sig->GetName();
//          if (nme.Contains( "_target" )) {
//             TString tit = sig->GetTitle();
//             sig->SetTitle( tit.ReplaceAll("Input variable", "Regression target" ) );
//          }
//       }
      sig->SetTitle( "" );            
      

      // finally plot and overlay
      Float_t sc = 1.1;
      if (countPad == 1) sc = 1.3;
      sig->SetMaximum( TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*sc );
      sig->Draw( "hist" );
      cPad->SetLeftMargin( 0.17 );

      sig->GetYaxis()->SetTitleOffset( 1.50 );
      if (!isRegression) {
         bgd->Draw("histsame");
         TString ytit = TString("(1/N) ") + sig->GetYaxis()->GetTitle();
         ytit = TString("Fraction of Events");
         sig->GetYaxis()->SetTitle( ytit ); // histograms are normalised
      }

      if (countPad == 1) sig->GetXaxis()->SetTitle("Leading Lepton p_{T} [GeV/c]");
      if (countPad == 2) sig->GetXaxis()->SetTitle("Trailing Lepton p_{T} [GeV/c]");

TH1* GetCentK(TDirectory* top, Double_t c1, Double_t c2, Int_t s,
	      TLegend* l)
{
  TString dname; dname.Form("cent%06.2f_%06.2f", c1, c2);
  dname.ReplaceAll(".", "d");
  TDirectory* d = top->GetDirectory(dname);
  if (!d) {
    Warning("GetCetnK", "Directory %s not found in %s",
	    dname.Data(), top->GetName());
    return;
  }

  TDirectory* det = d->GetDirectory("details");
  if (!det) {
    Warning("GetCetnK", "Directory details not found in %s",
	    d->GetName());
    d->ls();
    return;
  }

  TObject* o = det->Get("scalar");
  if (!o) {
    Warning("GetCetnK", "Object scalar not found in %s",
	    det->GetName());
    return;
  }

  if (!o->IsA()->InheritsFrom(TH1::Class())) {
    Warning("GetCetnK", "Object %s is not a TH1, but a %s",
	    o->GetName(), o->ClassName());
    return;
  }
  TH1* h = static_cast<TH1*>(o->Clone());
  Color_t col = cc[(s-1)%10];
  h->SetLineColor(col);
  h->SetMarkerColor(col);
  h->SetFillColor(col);
  h->SetFillStyle(1001);
  // h->SetTitle(Form("%5.2f-%5.2f%% #times %d", c1, c2, s));
  h->SetTitle(Form("%2.0f-%2.0f%% + %d", c1, c2, s-1));
  TF1* f = new TF1("", "[0]",-2.2,2.2);
  f->SetParameter(0,s-1);
  f->SetLineColor(col);
  f->SetLineStyle(7);
  f->SetLineWidth(1);
  // h->Scale(s);
  h->Add(f);
  h->GetListOfFunctions()->Add(f);
  f->SetParameter(0,s);
  for (Int_t i = 1; i <= h->GetNbinsX(); i++) {
    if (TMath::Abs(h->GetBinCenter(i)) > 2) {
      h->SetBinContent(i,0);
      h->SetBinError(i,0);
    }
  }
  
  TLegendEntry* e = l->AddEntry(h, h->GetTitle(), "f");
  e->SetFillColor(col);
  e->SetFillStyle(1001);
  e->SetLineColor(col);

  return h;
}

void DrawTwoInPad(TVirtualPad* p,
                  Int_t        sub,
                  TH1*         h1,
                  TH1*         h2,
		  Bool_t       ratio,
		  Bool_t       logy=false,
		  Bool_t       legend=false)
{
  TVirtualPad* pp = p->cd(sub);
  pp->SetRightMargin(0.02);
  pp->SetLeftMargin(0.10);
  TVirtualPad* ppp = pp;
  if (ratio) {
    pp->Divide(1,2,0,0);
    ppp = pp->cd(1);				
    ppp->SetRightMargin(0.02);
  }
  if (logy) ppp->SetLogy();
  TH1* hs[] = { h1, h2, 0 };
  if (h1->GetMaximum() < h2->GetMaximum()) {
    hs[0] = h2;
    hs[1] = h1;
  }
  TH1** ph = hs;
  Double_t size = (ratio ? 0.1 : 0.05);
  Double_t off  = (ratio ? 0.6 : 0.5);
  h1->SetFillStyle(3004);
  h2->SetFillStyle(3005);
  while (*ph) {
    TString opt("hist");
    if (ph != hs) opt.Append(" same");

    TH1* copy = (*ph)->DrawCopy(opt);
    copy->GetXaxis()->SetLabelSize(2*size);
    copy->GetYaxis()->SetLabelSize(size);
    copy->GetYaxis()->SetTitleSize(size);
    copy->GetYaxis()->SetTitleOffset(off);
    copy->SetYTitle(copy->GetTitle());
    copy->SetTitle("");
    copy->SetDirectory(0);
    ph++;
  }
  TString s1 = h1->GetYaxis()->GetTitle();
  TString s2 = h2->GetYaxis()->GetTitle();
  
  if (legend) { 
    TLegend* l = new TLegend(0.6, 0.1, 0.9, 0.9);
    l->SetBorderSize(0);
    TLegendEntry* e = l->AddEntry("dummy", s1, "lf");
    l->SetFillColor(kWhite);
    e->SetFillColor(kBlack);
    e->SetFillStyle(h1->GetFillStyle());
    e = l->AddEntry("dummy", s2, "lf");
    e->SetFillColor(kBlack);
    e->SetFillStyle(h2->GetFillStyle());
    l->Draw();
  }
  if (!ratio) return;
  ppp = pp->cd(2);
  ppp->SetRightMargin(0.02);
  TH1* r = static_cast<TH1*>(h1->Clone(Form("ratio%s", h1->GetName())));
  r->SetDirectory(0);
  r->SetTitle("");
  r->GetXaxis()->SetLabelSize(size);
  r->GetYaxis()->SetLabelSize(size);
  r->GetYaxis()->SetTitleSize(0.9*size);
  r->GetYaxis()->SetTitleOffset(0.9*off);
  r->SetMarkerStyle(20);
  r->SetMarkerColor(h1->GetFillColor()+1);
  r->SetFillStyle(3007);
  r->SetYTitle(Form("#frac{%s}{%s}", s1.Data(), s2.Data()));

  // r->Add(h2, -1);
  // r->Divide(h1);
  if (!r->IsA()->InheritsFrom(TProfile::Class())) {
    r->GetSumw2()->Set(0); // r->Sumw2(false);
    h2->GetSumw2()->Set(0); // h2->Sumw2(false);
  }
  r->Divide(h2);
  Printf("%s", r->GetName());
  for (UShort_t bin = 1; bin <= r->GetNbinsX(); bin++) {
    Printf("  bin # %2d: Diff=%g+/-%g", bin, r->GetBinContent(bin),
	   r->GetBinError(bin));
    r->SetBinError(bin, 0);
  }
  r->GetSumw2()->Set(0); //r->Sumw2(false);
  r->SetMarkerSize(4);
  r->SetMaximum(r->GetMaximum()*1.2);
  r->SetMinimum(r->GetMinimum()*0.8);
  r->Draw("hist text30");  
  p->Modified();
  p->Update();
  p->cd();
}

TObjArray* ProcessSummary(TObjArray* arrs, int icl, const char* pref)
{
  // Process TObjArray (e.g. for set of pt bins) of TObjArray of KMCTrackSummary objects:
  // pick the KMCTrackSummary for "summary class" icl (definition of acceptable track) and create
  // graphs vs bin.
  // These graphs are returned in new TObjArray
  //
  TString prefs = pref;
  if (!gs) gs = new TF1("gs","gaus",-1,1);
  //
  int nb = arrs->GetEntriesFast();
  TObjArray* sums = (TObjArray*) arrs->At(0);
  int nclass = sums->GetEntriesFast();
  if (icl>=nclass) {printf("summary set has %d classes only, %d requested\n",nclass,icl);return 0;}
  //
  KMCTrackSummary* sm = (KMCTrackSummary*)sums->At(icl);
  //
  TH1* h;
  //
  h = sm->GetHMCSigDCARPhi();  // MC resolution in transverse DCA 
  TGraphErrors * grSigD = 0;
  if (h) {
    grSigD = new TGraphErrors(nb);
    grSigD->SetName(Form("%s%s",prefs.Data(),h->GetName()));
    grSigD->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
  }
  //
  TGraphErrors * grSigZ = 0;
  h = sm->GetHMCSigDCAZ();    // MC resolution in Z DCA
  if (h) {
    grSigZ = new TGraphErrors(nb);
    grSigZ->SetName(Form("%s%s",prefs.Data(),h->GetName()));
    grSigZ->SetTitle(Form("%s%s",prefs.Data(),h->GetTitle()));
  }
  //
  TGraphErrors * grSigAD = 0; // anaitical estimate for resolution in transverse DCA 
  {
    grSigAD = new TGraphErrors(nb);
    grSigAD->SetName(Form("%s%s",prefs.Data(),"sigmaDan"));
    grSigAD->SetTitle(Form("%s%s",prefs.Data(),"#sigmaD an"));
  }
  //
  TGraphErrors * grSigAZ = 0; // anaitical estimate for resolution in Z DCA 
  {
    grSigAZ = new TGraphErrors(nb);
    grSigAZ->SetName(Form("%s%s",prefs.Data(),"sigmaZan"));
    grSigAZ->SetTitle(Form("%s%s",prefs.Data(),"#sigmaZ an"));
  }
  //
  //
  TGraphErrors * grSigPt = 0; // MC res. in pt
  {
    grSigPt = new TGraphErrors(nb);
    grSigPt->SetName(Form("%s%s",prefs.Data(),"sigmaPt"));
    grSigPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt"));
  }
  //
  TGraphErrors * grSigAPt = 0; // analitycal res. in pt
  {
    grSigAPt = new TGraphErrors(nb);
    grSigAPt->SetName(Form("%s%s",prefs.Data(),"sigmaPtan"));
    grSigAPt->SetTitle(Form("%s%s",prefs.Data(),"#sigmaPt an"));
  }

  //
  TGraphErrors * grEff = 0; // MC efficiency
  {
    grEff = new TGraphErrors(nb);
    grEff->SetName(Form("%s_rate",prefs.Data()));
    grEff->SetTitle(Form("%s Rate",prefs.Data()));
  }
  //
  TGraphErrors * grUpd = 0; // number of Kalman track updates
  {
    grUpd = new TGraphErrors(nb);
    grUpd->SetName(Form("%s_updCalls",prefs.Data()));
    grUpd->SetTitle(Form("%s Updates",prefs.Data()));
  }
  //
  for (int ib=0;ib<nb;ib++) {
    sums = (TObjArray*) arrs->At(ib);
    sm = (KMCTrackSummary*)sums->At(icl);
    KMCProbe& prbRef = sm->GetRefProbe();
    KMCProbe& prbAn  = sm->GetAnProbe();
  
    double pt = prbRef.Pt();
    //
    if (grSigAD) {
      grSigAD->SetPoint(ib, pt,prbAn.GetSigmaY2()>0 ? TMath::Sqrt(prbAn.GetSigmaY2()) : 0.);
    }
    //
    if (grSigAZ) {
      grSigAZ->SetPoint(ib, pt,prbAn.GetSigmaZ2()>0 ? TMath::Sqrt(prbAn.GetSigmaZ2()) : 0.);
    }
    //
    if (grSigAPt) {
      double pts = TMath::Sqrt(prbAn.GetSigma1Pt2());
      grSigAPt->SetPoint(ib, pt,pts>0 ? pts*pt : 0.);
    }
    //
//.........这里部分代码省略.........

void comparisonJetMCData(string plot,int rebin){
  string tmp;

  string dir="/gpfs/cms/data/2011/Observables/Approval/";
	
  if (isAngularAnalysis){
    mcfile=dir+"MC_zjets"+version;
    back_w=dir+"MC_wjets"+version;
    back_ttbar=dir+"MC_ttbar"+version;
    WW=dir+"MC_diW"+version;
    ZZ=dir+"MC_siZ"+version;
    WZ=dir+"MC_diWZ"+version;
    datafile=dir+"DATA"+version;
    mcfiletau=dir+"MC_zjetstau"+version;
  }
  // List of files

  TFile *dataf = TFile::Open(datafile.c_str()); //data file
  TFile *mcf = TFile::Open(mcfile.c_str()); //MC file
  TFile *mcftau = TFile::Open(mcfiletau.c_str()); //MC file
  TFile *ttbarf = TFile::Open(back_ttbar.c_str()); //MC background file
  TFile *wf = TFile::Open(back_w.c_str());


  TFile *qcd23emf = TFile::Open(qcd23em.c_str());
  TFile *qcd38emf = TFile::Open(qcd38em.c_str());
  TFile *qcd817emf = TFile::Open(qcd817em.c_str());

  TFile *qcd23bcf = TFile::Open(qcd23bc.c_str());
  TFile *qcd38bcf = TFile::Open(qcd38bc.c_str());
  TFile *qcd817bcf = TFile::Open(qcd817bc.c_str());

  TFile *WZf = TFile::Open(WZ.c_str());
  TFile *ZZf = TFile::Open(ZZ.c_str());
  TFile *WWf = TFile::Open(WW.c_str());


  // Canvas
  if (CanvPlot) delete CanvPlot;
  CanvPlot = new TCanvas("CanvPlot","CanvPlot",0,0,800,600);

  // Getting, defining ...
  dataf->cd("validationJEC");
  if (isMu && isAngularAnalysis) dataf->cd("validationJECmu");
	
  TObject * obj;
  gDirectory->GetObject(plot.c_str(),obj);

  TH1 *data;
  TH2F *data2; 
  TH1D *data3; 

  THStack *hs = new THStack("hs","Total MC");


  int flag=-1;
  if ((data = dynamic_cast<TH1F *>(obj)) ){
    flag=1;
    gROOT->Reset();
    gROOT->ForceStyle();
    gStyle->SetPadRightMargin(0.03);
    gPad->SetLogy(1);
    gPad->Modified();
    gPad->Update();
  }
  if ((data2 = dynamic_cast<TH2F *>(obj)) ){
    flag=2;
    gStyle->SetPalette(1);
    gStyle->SetPadRightMargin(0.15);
    gPad->Modified();
  }



  //===================
  // Dirty jobs :)
  if (flag==1){
    CanvPlot->cd();
    TPad *pad1 = new TPad("pad1","pad1",0.01,0.33,0.99,0.99);
    pad1->Draw();
    pad1->cd();
    pad1->SetTopMargin(0.1);
    pad1->SetBottomMargin(0.01);
    pad1->SetRightMargin(0.1);
    pad1->SetFillStyle(0);
    pad1->SetLogy(1);
    TString str=data->GetTitle();
    if (str.Contains("jet") && !str.Contains("zMass") && !str.Contains("Num") && !str.Contains("Eta") && !str.Contains("Phi") && !str.Contains("eld") && !str.Contains("meanPtZVsNjet")) {
      if (!isAngularAnalysis) rebin=1;
    }

    //======================
    // DATA
    Double_t dataint = data->Integral();
    data->SetLineColor(kBlack);
    data->Rebin(rebin);
    if(str.Contains("nJetVtx")) data->GetXaxis()->SetRangeUser(0,10);	
    if(str.Contains("zMass")) data->GetXaxis()->SetRangeUser(70,110);	
    data->SetMinimum(1.);
    data->Sumw2();
//.........这里部分代码省略.........

//.........这里部分代码省略.........
for (Int_t iev=0; iev< nevent; iev++) {
  // get entry	
  t1->GetEntry(iev);
  t2->GetEntry(iev);
		
  nMCTracks = MCTracks->GetEntriesFast();
  nPoints = TutorialDetPoints->GetEntriesFast();
  nHits = TutorialDetHits->GetEntriesFast();
    		
  cout << " Event"  << iev << ":";
  cout << nMCTracks << " MC tracks ";
  cout << nPoints << "  points ";
  cout << nHits << " Hits "<<endl;
    
  // Hit loop					
  for (Int_t j =0; j<nHits; j++) {
    Hit   = (FairTutorialDet4Hit*) TutorialDetHits->At(j);
    Int_t l = Hit->GetRefIndex();
    Point = (FairTutorialDet4Point*) TutorialDetPoints->At(l);

    // Point info
    x_poi	= Point -> GetX();
    y_poi	= Point -> GetY();
    z_poi	= Point -> GetZ();

    // Hit	info		
    x_hit	= Hit->GetX();
    y_hit	= Hit->GetY();
    z_hit	= Hit->GetZ();
    dy_hit	= Hit->GetDy();
//    Int_t flg_hit = Hit->GetFlag();

    Float_t delta_x	= x_poi - x_hit;
    Float_t delta_y	= y_poi - y_hit;
    Float_t delta_z	= z_poi - z_hit;
			
    dxx ->Fill(x_poi,delta_x);
    dyy ->Fill(y_poi,delta_y);
    pullx ->Fill(delta_x);
    pully ->Fill(delta_y);
    pullz ->Fill(delta_z);
    pointx ->Fill(x_hit);
    pointy ->Fill(y_hit);
		
  } //  Hit loop end

} // event loop end 


// save histos to file
// TFile *fHist = TFile::Open("data/auaumbias.hst.root","RECREATE");
 cout << "Processing done, outflag =" <<iout << endl;
 if (iout==1){
   fHist->Write();
   if(0){ // explicit writing 
     TIter next(gDirectory->GetList());
     TH1 *h;
     TObject* obj;
     while(obj= (TObject*)next()){
       if(obj->InheritsFrom(TH1::Class())){
         h = (TH1*)obj;
         cout << "Write histo " << h->GetTitle() << endl;
         h->Write();
       }
     }
   }
   fHist->ls();
   fHist->Close();
  }

  // -----   Finish   -------------------------------------------------------

  cout << endl << endl;

  // Extract the maximal used memory an add is as Dart measurement
  // This line is filtered by CTest and the value send to CDash
  FairSystemInfo sysInfo;
  Float_t maxMemory=sysInfo.GetMaxMemory();
  cout << "<DartMeasurement name=\"MaxMemory\" type=\"numeric/double\">";
  cout << maxMemory;
  cout << "</DartMeasurement>" << endl;

  timer.Stop();
  Double_t rtime = timer.RealTime();
  Double_t ctime = timer.CpuTime();

  Float_t cpuUsage=ctime/rtime;
  cout << "<DartMeasurement name=\"CpuLoad\" type=\"numeric/double\">";
  cout << cpuUsage;
  cout << "</DartMeasurement>" << endl;

  cout << endl << endl;
  cout << "Output file is "    << outFile << endl;
  cout << "Parameter file is " << parFile << endl;
  cout << "Real time " << rtime << " s, CPU time " << ctime
       << "s" << endl << endl;
  cout << "Macro finished successfully." << endl;

  // ------------------------------------------------------------------------
}

string FullTitle(const TH1 &h) {
    return string(h.GetTitle())
           +";"+h.GetXaxis()->GetTitle()
           +";"+h.GetYaxis()->GetTitle();
}

  /** 
   * Process a single eta bin 
   * 
   * @param measured     Input collection of measured data
   * @param corrections  Input collection of correction data
   * @param method       Unfolding method to use 
   * @param regParam     Regularisation parameter
   * @param out          Output directory. 
   *
   * @return Stack of histograms or null 
   */
  THStack* ProcessBin(TCollection* measured, 
		      TCollection* corrections, 
		      UInt_t       method,
		      Double_t     regParam, 
		      TDirectory*  out)
  {
    Printf("   Processing %s ...", measured->GetName());
    // Try to get the data 
    TH1* inRaw    = GetH1(measured,    "rawDist");
    TH1* inTruth  = GetH1(corrections, "truth");
    TH1* inTruthA = GetH1(corrections, "truthAccepted");
    TH1* inTrgVtx = GetH1(corrections, "triggerVertex");
    TH2* inResp   = GetH2(corrections, "response");
    if (!inRaw || !inTruth || !inTruthA || !inTrgVtx || !inResp) 
      return 0;
    
    // Make output directory
    TDirectory* dir = out->mkdir(measured->GetName());
    dir->cd();

    // Copy the input to the output 
    TH1* outRaw    = static_cast<TH1*>(inRaw    ->Clone("measured"));
    TH1* outTruth  = static_cast<TH1*>(inTruth  ->Clone("truth"));
    TH1* outTruthA = static_cast<TH1*>(inTruthA ->Clone("truthAccepted"));
    TH1* outTrgVtx = static_cast<TH1*>(inTrgVtx ->Clone("triggerVertex"));
    TH2* outResp   = static_cast<TH2*>(inResp   ->Clone("response"));

    // Make our response matrix 
    RooUnfoldResponse matrix(0, 0, inResp);
    
    // Store regularization parameter 
    Double_t             r        = regParam;
    RooUnfold::Algorithm algo     = (RooUnfold::Algorithm)method;
    RooUnfold*           unfolder = RooUnfold::New(algo, &matrix, inRaw, r);
    unfolder->SetVerbose(0);

    // Do the unfolding and get the result
    TH1* res = unfolder->Hreco();
    res->SetDirectory(0);

    // Make a copy to store on the output 
    TH1* outUnfold = static_cast<TH1*>(res->Clone("unfolded"));
    TString tit(outUnfold->GetTitle());
    tit.ReplaceAll("Unfold Reponse matrix", "Unfolded P(#it{N}_{ch})");
    outUnfold->SetTitle(tit);

    // Clone the unfolded results and divide by the trigger/vertex
    // bias correction
    TH1* outCorr   = static_cast<TH1*>(outUnfold->Clone("corrected"));
    outCorr->Divide(inTrgVtx);
    tit.ReplaceAll("Unfolded", "Corrected");
    outCorr->SetTitle(tit);

    // Now normalize the output to integral=1 
    TH1*  hists[] = { outRaw, outUnfold, outCorr, 0 };
    TH1** phist   = hists;
    while (*phist) { 
      TH1* h = *phist;
      if (h) { 
	Double_t intg = h->Integral(1, h->GetXaxis()->GetXmax());
	h->Scale(1. / intg, "width");
      }
      phist++;
    }
    
    // And make ratios
    TH1* ratioTrue = static_cast<TH1*>(outCorr->Clone("ratioCorrTruth"));
    tit = ratioTrue->GetTitle();
    tit.ReplaceAll("Corrected", "Corrected/MC 'truth'");
    ratioTrue->SetTitle(tit);
    ratioTrue->Divide(outTruth);
    ratioTrue->SetYTitle("P_{corrected}(#it{N}_{ch})/P_{truth}(#it{N}_{ch})");

    TH1* ratioAcc  = static_cast<TH1*>(outUnfold->Clone("ratioUnfAcc"));
    tit = ratioAcc->GetTitle();
    tit.ReplaceAll("Unfolded", "Unfolded/MC selected");
    ratioAcc->SetTitle(tit);
    ratioAcc->Divide(outTruthA);
    ratioAcc->SetYTitle("P_{unfolded}(#it{N}_{ch})/P_{MC}(#it{N}_{ch})");
    

    // Make a stack 
    tit = measured->GetName();
    tit.ReplaceAll("m", "-");
    tit.ReplaceAll("p", "+");
    tit.ReplaceAll("d", ".");
    tit.ReplaceAll("_", "<#it{#eta}<");
    THStack* stack = new THStack("all", tit);
    stack->Add(outTruth,  "E2");
//.........这里部分代码省略.........

//.........这里部分代码省略.........
                             Form("probability_%s",methodTitle.Data()) :
                             (htype == CompareType) ? 
                             Form("comparison_%s",methodTitle.Data()) :
                             Form("rarity_%s",methodTitle.Data()));

            c = new TCanvas( Form("canvas%d", countCanvas+1), ctitle, 
                             countCanvas*50+200, countCanvas*20, width, (Int_t)width*0.78 ); 
    
            // set the histogram style
            TMVAGlob::SetSignalAndBackgroundStyle( sig, bgd );
   
            // normalise both signal and background
            TMVAGlob::NormalizeHists( sig, bgd );
   
            // frame limits (choose judicuous x range)
            Float_t nrms = 4;
            cout << "--- Mean and RMS (S): " << sig->GetMean() << ", " << sig->GetRMS() << endl;
            cout << "--- Mean and RMS (B): " << bgd->GetMean() << ", " << bgd->GetRMS() << endl;
            Float_t xmin = TMath::Max( TMath::Min(sig->GetMean() - nrms*sig->GetRMS(), 
                                                  bgd->GetMean() - nrms*bgd->GetRMS() ),
                                       sig->GetXaxis()->GetXmin() );
            Float_t xmax = TMath::Min( TMath::Max(sig->GetMean() + nrms*sig->GetRMS(), 
                                                  bgd->GetMean() + nrms*bgd->GetRMS() ),
                                       sig->GetXaxis()->GetXmax() );
            Float_t ymin = 0;
            Float_t maxMult = (htype == CompareType) ? 1.3 : 1.2;
            Float_t ymax = TMath::Max( sig->GetMaximum(), bgd->GetMaximum() )*maxMult;
   
            // build a frame
            Int_t nb = 500;
            TString hFrameName(TString("frame") + methodTitle);
            TObject *o = gROOT->FindObject(hFrameName);
            if(o) delete o;
            TH2F* frame = new TH2F( hFrameName, sig->GetTitle(), 
                                    nb, xmin, xmax, nb, ymin, ymax );
            frame->GetXaxis()->SetTitle( methodTitle + ((htype == MVAType || htype == CompareType) ? " response" : "") );
            if      (htype == ProbaType  ) frame->GetXaxis()->SetTitle( "Signal probability" );
            else if (htype == RarityType ) frame->GetXaxis()->SetTitle( "Signal rarity" );
            frame->GetYaxis()->SetTitle("Normalized");
            TMVAGlob::SetFrameStyle( frame );
   
            // eventually: draw the frame
            frame->Draw();  
    
            c->GetPad(0)->SetLeftMargin( 0.105 );
            frame->GetYaxis()->SetTitleOffset( 1.2 );

            // Draw legend               
            TLegend *legend= new TLegend( c->GetLeftMargin(), 1 - c->GetTopMargin() - 0.12, 
                                          c->GetLeftMargin() + (htype == CompareType ? 0.40 : 0.3), 1 - c->GetTopMargin() );
            legend->SetFillStyle( 1 );
            legend->AddEntry(sig,TString("Signal")     + ((htype == CompareType) ? " (test sample)" : ""), "F");
            legend->AddEntry(bgd,TString("Background") + ((htype == CompareType) ? " (test sample)" : ""), "F");
            legend->SetBorderSize(1);
            legend->SetMargin( (htype == CompareType ? 0.2 : 0.3) );
            legend->Draw("same");

            // overlay signal and background histograms
            sig->Draw("samehist");
            bgd->Draw("samehist");
   
            if (htype == CompareType) {
               // if overtraining check, load additional histograms
               TH1* sigOv = 0;
               TH1* bgdOv = 0;

  //____________________________________________________________________
  void DrawResCollection(TCollection* top, const TString& name)
  {
    TCollection* c = GetCollection(top, name, false);
    if (!c) return;

    THStack* s = GetStack(c, "all");
    s->SetTitle("");
    DrawInPad(fBody, 0, s, "nostack", kLogy);
    TLegend* l = new TLegend(.5, .75, .98, .98, "P(#it{N}_{ch})");
    l->SetBorderSize(0);
    // l->SetBorderMode(0);
    l->SetFillColor(0);
    l->SetFillStyle(0);
    TIter next(s->GetHists());
    TH1*  h = 0;
    Bool_t hasTrue = false;
    while ((h = static_cast<TH1*>(next()))) { 
      TString n(h->GetTitle());
      if (n.BeginsWith("True")) { hasTrue = true; continue; }
      n.ReplaceAll("Raw P(#it{N}_{ch}) in ", "");
      TLegendEntry* e = l->AddEntry("dummy", n, "p");
      e->SetMarkerStyle(h->GetMarkerStyle());
    }
    if (hasTrue) {
      TLegendEntry* e = l->AddEntry("dummy", "Raw", "p");
      e->SetMarkerStyle(20);
      e->SetMarkerColor(kRed+1);
      e = l->AddEntry("dummy", "MC truth", "p");
      e->SetMarkerStyle(24);
      e->SetMarkerColor(kBlue+1);
      e = l->AddEntry("dummy", "MC truth selected", "p");
      e->SetMarkerStyle(24);
      e->SetMarkerColor(kOrange+1);
    }
    fBody->cd();
    l->Draw();
    
    PrintCanvas(Form("%s results", name.Data()));

    // return;

    TIter       nextO(c);
    TObject*    o = 0;
    while ((o = nextO())) { 
      Double_t etaMin = 999;
      Double_t etaMax = 999;
      TCollection* bin = GetEtaBin(o, etaMin, etaMax);
      if (!bin) continue;

      fBody->Divide(2,3);
      DrawInPad(fBody, 1, GetH1(bin, "rawDist"),      "",     kLogy);
      DrawInPad(fBody, 1, GetH1(bin, "truthAccepted", false),
		"same", kSilent);
      DrawInPad(fBody, 1, GetH1(bin, "truth", false),"same", kSilent|kLegend);
      DrawInPad(fBody, 2, GetH1(bin, "coverage"));
      DrawInPad(fBody, 3, GetH2(bin, "corr"),            "colz");
      DrawInPad(fBody, 4, GetH2(bin, "response", false), "colz", kLogz|kSilent);
      DrawInPad(fBody, 5, GetH1(bin, "triggerVertex", false), "", kSilent);
      
      PrintCanvas(Form("%+5.1f < #eta < %+5.1f", etaMin, etaMax));
    }
  }