C++ TDirectory::GetDirectory方法代码示例

/**
 * @brief Write all TObjects from a given TCollection into a certain directory structure in the file
 *
 * @param col pointer to a TCollection-based container
 * @param dirname name of a directory inside the output file to which the objects should be written
 * @param subdirname optional name of a subdirectory inside dirname to which the objects should be written
 *
 * This method whites all TObject-based objects contained in the TCollection-based container (see ROOT documentation)
 * into a directory whose name is given by dirname inside the output file. If dirname does not exist
 * in the output file, it will be created. Otherwise, contents of the col collection will be appended to an existing
 * directory.
 *
 * If the optional subdirectory name is specified (subdirname parameter, defaults to empty string) then the
 * contents of the collection will be written to "dirname/subdirname". If the "subdirname" directory does not
 * exist inside the "dirname" directory, it will be created.
 *
 */
void JPetWriter::writeCollection(const TCollection* col, const char* dirname, const char* subdirname)
{

  TDirectory* current =  fFile->GetDirectory(dirname);

  if (!current) {
    current = fFile->mkdir(dirname);
  }

  assert(current);

  // use a subdirectory if requested by user
  if (!std::string(subdirname).empty()) {

    if (current->GetDirectory(subdirname)) {
      current = current->GetDirectory(subdirname);
    } else {
      current = current->mkdir(subdirname);
    }
  }

  assert(current);

  current->cd();

  TIterator* it = col->MakeIterator();

  TObject* obj;
  while ((obj = it->Next())) {
    obj->Write();
  }

  fFile->cd();
}

void store(){
	 if (!TClass::GetDict("NBD")) {
                gROOT->ProcessLine(".L /afs/cern.ch/user/q/qixu/CMSSW_6_2_5/src/Centrality/NBD_Glauber_fit/NBD/NBDclass.C+");
        }
	TH1::AddDirectory(kFALSE); 
	int Gth=atoi(getenv("GTH"));
	int sth=atoi(getenv("STH"));
	NBD *l;
	struct para1 var;
        if(sth==0){TString dirname = "std"; var = var1; if(Gth<nGlau) l = new NBD(datafile,stdGlaulist[Gth],histoname);else l = new NBD(datafile,stdGlaulist[0],histoname);}
        else if(sth==1){TString dirname ="Gri055"; var = var2; if(Gth<nGlau) l = new NBD(datafile,Gri055Glaulist[Gth],histoname);else l = new NBD(datafile,Gri055Glaulist[0],histoname);}
        else {TString dirname ="Gri101"; var = var3; if(Gth<nGlau) l = new NBD(datafile,Gri101Glaulist[Gth],histoname);else l = new NBD(datafile,Gri101Glaulist[0],histoname);}
	l->initmu(var.mumin,var.mumax,var.mustep);
     	l->initk(var.kmin,var.kmax,var.kstep);
	if(Gth<nGlau)
        	l->initx(var.xmin,var.xmax);
	else if(Gth-nGlau==0)
		l->initx(var.xmin-binshift,var.xmax);
	else
		l->initx(var.xmin+binshift,var.xmax);
//	l->fit();

 	if(sth==0) l->assign(bestlist1[Gth].mubest,bestlist1[Gth].kbest);
	else if(sth==1) l->assign(bestlist2[Gth].mubest,bestlist2[Gth].kbest);
	else l->assign(bestlist3[Gth].mubest,bestlist3[Gth].kbest);

        l->initN(bin,N,method);
        l->calcvar();
	TFile *outfile = new TFile(outG,"Update");
	if(Gth==0 && sth==0){
	l->dataname.Write("dataname",TObject::kOverwrite);
	l->histoname.Write("histoname",TObject::kOverwrite);
	}
	TDirectory *dir = outfile->GetDirectory(dirname);
	if (!dir) {outfile->mkdir(dirname);	TDirectory *dir = outfile->GetDirectory(dirname);}
	dir->cd();
	TString name;
	if(Gth==0)
		name = "G0";
	else if(Gth<nGlau)
		name = Form("Glau_%d",Gth);
	else
		name = Form("bin_%d",Gth-nGlau+1);
	TDirectory *subdir = dir->GetDirectory(name);
	if(!subdir) {dir->mkdir(name);	TDirectory *subdir = dir->GetDirectory(name);}
	subdir->cd();
	l->method.Write("method",TObject::kOverwrite);
	l->Glaubername.Write("Glaubername",TObject::kOverwrite);
	l->xmin.Write("xmin",TObject::kOverwrite);l->xmax.Write("xmax",TObject::kOverwrite);
	l->mubest.Write("mubest",TObject::kOverwrite);l->kbest.Write("kbest",TObject::kOverwrite);
	l->chis.Write("chis",TObject::kOverwrite);l->Ndf.Write("Ndf",TObject::kOverwrite);
	l->NcollAver.Write("NcollAver",TObject::kOverwrite);l->NpartAver.Write("NpartAver",TObject::kOverwrite);l->BAver.Write("BAver",TObject::kOverwrite);
	l->centbin.Write("centbin",TObject::kOverwrite);	l->kpoint.Write("kpoint",TObject::kOverwrite);
	l->centbin_.Write("centbin_",TObject::kOverwrite);	l->kpoint_.Write("kpoint_",TObject::kOverwrite);
	l->Npartdis->Write();
        l->Grgrid->Write("Grgrid",TObject::kOverwrite);
        outfile->Close();

}

LjmetEventContent::~LjmetEventContent(){

  // save tree
  mpTree->Write();

  // save histograms
  std::map<std::string,std::map<std::string,LjmetEventContent::HistMetadata> >::iterator iMod;
  std::map<std::string,LjmetEventContent::HistMetadata>::iterator iHist;
  TDirectory * rootDir = mpTree->GetDirectory();
  for (iMod=mDoubleHist.begin();iMod!=mDoubleHist.end();++iMod){
    TDirectory * _dir = 0;
    _dir = rootDir->mkdir( iMod->first.c_str() );
    if (!_dir) rootDir->GetDirectory( iMod->first.c_str() );
    for (iHist=iMod->second.begin();iHist!=iMod->second.end();++iHist){
      std::cout << mLegend
		<< "Saving " << iMod->first << "/"
		<< iHist->second.GetName() << std::endl;
      _dir->cd();
      iHist->second.GetHist( )->SetDirectory(_dir);
      iHist->second.GetHist( )->Print();
      iHist->second.GetHist( )->Write();
    }
  }

}

void store(int type, int sth=0, int Gth=0){
	TH1::AddDirectory(kFALSE); 
	NBD *l;
	struct para1 var;
        TString dirname;
        double sss=0.1;
	if(sth==0){
	dirname = "std";
	if(Gth<nGlau)
            l = new NBD(datafile,stdGlaulist[Gth],histoname);
	else
            l = new NBD(datafile,stdGlaulist[0],histoname);
	var = var1[Gth];
	l->initmu(bestlist1[Gth].mubest-sss,bestlist1[Gth].mubest+sss,sss/5);
	l->initk(bestlist1[Gth].kbest-sss,bestlist1[Gth].kbest+sss,sss/5);
	l->initf(bestlist1[Gth].fbest-sss,bestlist1[Gth].fbest+sss,sss/5);
	}
        
	else if(sth==1){
	dirname ="Gri055"; 
	if(Gth<nGlau)
            l = new NBD(datafile,Gri055Glaulist[Gth],histoname);
	else
            l = new NBD(datafile,Gri055Glaulist[0],histoname);
	var = var2[Gth];
	l->initmu(bestlist2[Gth].mubest-sss,bestlist2[Gth].mubest+sss,sss/5);
	l->initk(bestlist2[Gth].kbest-sss,bestlist2[Gth].kbest+sss,sss/5);
	l->initf(bestlist2[Gth].fbest-sss,bestlist2[Gth].fbest+sss,sss/5);
	}
        
	else {
	dirname ="Gri101";
	if(Gth<nGlau)
            l = new NBD(datafile,Gri101Glaulist[Gth],histoname);
	else
            l = new NBD(datafile,Gri101Glaulist[0],histoname);
        var = var3[Gth];
	l->initmu(bestlist3[Gth].mubest-sss,bestlist3[Gth].mubest+sss,sss/5);
	l->initk(bestlist3[Gth].kbest-sss,bestlist3[Gth].kbest+sss,sss/5);
	l->initf(bestlist3[Gth].fbest-sss,bestlist3[Gth].fbest+sss,sss/5);
	}

	//l->initmu(var.mumin,var.mumax,var.mustep);
    	//l->initk(var.kmin,var.kmax,var.kstep);
    	//l->initf(var.fmin,var.fmax,var.fstep);
	if(Gth<nGlau)
        	l->initx(var.xmin,var.xmax);
	else if(Gth-nGlau==0)
		l->initx(var.xmin-binshift,var.xmax);
	else
		l->initx(var.xmin+binshift,var.xmax);
if(type==0)
	l->fit();	
if(type==1){
 	if(sth==0) l->assign(bestlist1[Gth].mubest,bestlist1[Gth].kbest,bestlist1[Gth].fbest);
	else if(sth==1) l->assign(bestlist2[Gth].mubest,bestlist2[Gth].kbest,bestlist2[Gth].fbest);
	else l->assign(bestlist3[Gth].mubest,bestlist3[Gth].kbest,bestlist3[Gth].fbest);

        l->initN(bin,N,method);
        l->calcvar();
	TFile *outfile = new TFile(outG,"Update");
	if(Gth==0 && sth==0){
	l->dataname.Write("dataname",TObject::kOverwrite);
	l->histoname.Write("histoname",TObject::kOverwrite);
	}
	TDirectory *dir = (TDirectory*)outfile->GetDirectory(dirname);
	if (!dir) {outfile->mkdir(dirname);	dir = (TDirectory*)outfile->GetDirectory(dirname);}
	dir->cd();
	TString name;
	if(Gth==0)
		name = "G0";
	else if(Gth<nGlau)
		name = Form("Glau_%d",Gth);
	else
		name = Form("bin_%d",Gth-nGlau+1);
	TDirectory *subdir = (TDirectory*)dir->GetDirectory(name);
	if(!subdir) {dir->mkdir(name);	subdir = (TDirectory*)dir->GetDirectory(name);}
	subdir->cd();
	l->method.Write("method",TObject::kOverwrite);
	l->Glaubername.Write("Glaubername",TObject::kOverwrite);
	l->xmin.Write("xmin",TObject::kOverwrite);l->xmax.Write("xmax",TObject::kOverwrite);
	l->mubest.Write("mubest",TObject::kOverwrite);l->kbest.Write("kbest",TObject::kOverwrite);l->fbest.Write("fbest",TObject::kOverwrite);
	l->chis.Write("chis",TObject::kOverwrite);l->Ndf.Write("Ndf",TObject::kOverwrite);
	l->NcollAver.Write("NcollAver",TObject::kOverwrite);l->NpartAver.Write("NpartAver",TObject::kOverwrite);l->BAver.Write("BAver",TObject::kOverwrite);
	l->centbin.Write("centbin",TObject::kOverwrite);	l->kpoint.Write("kpoint",TObject::kOverwrite);
	l->centbin_.Write("centbin_",TObject::kOverwrite);	l->kpoint_.Write("kpoint_",TObject::kOverwrite);
	l->Npartdis->Write("Npartdis",TObject::kSingleKey | TObject::kOverwrite);
	l->Grgrid->Write("Grgrid",TObject::kOverwrite);
        outfile->Close();
}
}

void SubtractBgndRuns(TDirectory *fbg,
                      TDirectory *fprod,
                      const Char_t *foutname="",
                      const std::string& bl_list="")
{
    std::set<UInt_t> blacklist = GetBlackedChannels(bl_list);

    string basename = FindBaseName(fprod->GetName());

    vector<TH1*> histbg   = ROOTUtils::GetAllTH1InDirectory(fbg);
    ROOTUtils::EnforceProperOrdering(histbg);
    //GetAllHists(*fbg, false);
    FilterOutTDC(histbg);
    vector<TH1*> histprod = ROOTUtils::GetAllTH1InDirectory(fprod,
                                                            Form("^((%s[0-9]+_shifted_norm_vclk)|(%s[0-9]+_norm_vclk_shifted))$",basename.data(),basename.data()));
    ROOTUtils::EnforceProperOrdering(histprod);
    //GetAllHists(*fprod, true);
    //  for_each(histprod.begin(), histprod.end(), InspectSingleHistErrors);
    FilterOutTDC(histprod);

    EnforceAxisRangeConsistency(histbg, histprod);

    TDirectory *fout=0;
    string foname = foutname;
    if (foname.length()>0)
    {

        fout = TFile::Open(foname.data(),"UPDATE");
    }
    else
    {
        TDirectory *basedir = GetBaseDir(fprod);
        if (dynamic_cast<TDirectoryFile*>(basedir)==0)
        {
            std::cerr << "Cannot identify the file" << std::endl;
            return;
        }

        string dname = FormOutputSubDirName(fprod->GetName());
        fout = basedir->GetDirectory(dname.data());
        if (fout==0)
            fout = basedir->mkdir(dname.data());

    }

    TH1 *hbg=0, *h=0, *h_cl=0;
    for (UInt_t i=0; i<min(histbg.size(), histprod.size()); i++)
    {
        hbg = histbg.at(i);
        h   = histprod.at(i);


        h_cl = static_cast<TH1*>(h->Clone(Form("normed_bgsubbed_%s%i",basename.data(), i)));

        if (blacklist.find(i)==blacklist.end())
        {
            cout << histbg.at(i)->GetName() << " subbed from " << histprod.at(i)->GetName() << endl;
            h_cl->SetTitle(Form("Overnight Bgnd Subtracted ADC%i",i));
            h_cl->Add(hbg,-1.0);
            // The next line is only valid if Sumw2 is not turned on
            //  CalcErrors(h_cl, hbg, h);
        }
        else
        {
            cout << histbg.at(i)->GetName() << " blacklisted ... no subtraction" << endl;
            h_cl->SetTitle(Form("ADC%i Without Subtraction",i));
        }

        fout->cd();
        h_cl->Write("",TObject::kOverwrite);
    }
    
    std::cout << "Results saved in " << fout->GetName() << std::endl;

    if (foname.length()>0)
        fout->Close();

}

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;
}

int InputForLimits(){
  TH1::SetDefaultSumw2(true);
  if(pcp)cout<<"going to set inputs"<<endl;

  Int_t NBR = 3;
  Float_t BR[] = { 1., 0.75, 0.5};

  TFile* bkgFile = new TFile( "../../BkgPrediction/BkgPrediction.root", "READ"); 
  TTree* bkgTree;
  bkgFile->GetObject( "ElAndMu", bkgTree);
  
  Float_t bkg = 0.;
  Float_t bkgTotUnc = 0.;

  Float_t obs = 0.;

  bkgTree->SetBranchAddress( "srData", &obs);
  bkgTree->SetBranchAddress( "srAllBkgCorr", &bkg);
  bkgTree->SetBranchAddress( "TotUnc", &bkgTotUnc);

  TFile* sigFile = new TFile( "../../SignalSystematics/SignalSys.root", "READ"); 

  std::vector<std::vector<TString> > sysColl;
  std::vector<TString> sys;
  sys.push_back(TString("JES_Up"));
  sys.push_back(TString("JES_Down"));
  sysColl.push_back(sys);

  sys.clear();
  sys.push_back(TString("BTagReweight_UpBC"));
  sys.push_back(TString("BTagReweight_DownBC"));
  sysColl.push_back(sys);

  sys.clear();
  sys.push_back(TString("BTagReweight_UpLight"));
  sys.push_back(TString("BTagReweight_DownLight"));
  sysColl.push_back(sys);

  std::vector<TString> decayMode;
  decayMode.push_back(TString("tt"));
  decayMode.push_back(TString("tb"));
  decayMode.push_back(TString("bb"));
		      
  Systematics* systematics[3];
  TString dirname;
  TString histoname;

  TDirectory* srDir;
  TDirectory* histoDir;

  TFile* outFile = new TFile( "InputForLimits.root", "RECREATE");
  TDirectory* outBRDir;
  TDirectory* outSRDir;

  TH1F* datah = new TH1F( "data", "data", 1, 0., 1.);
  TH1F* bkgh = new TH1F( "bkg", "bkg", 1, 0., 1.);

  TH2F* sigh;
  TH2F* sig_toth;
  TH2F* effh;
  
  TH2F* sigLh;
  TH2F* effLh;

  TH2F* sigRh;
  TH2F* effRh;
  
  TH2F* jesh;
  TH2F* btagBCh;
  TH2F* btagLighth;
  TH2F* btagh;
  
  TH2F* sysh;

  TH2F* unch;

  TH2F* jesPercenth;
  TH2F* btagBCPercenth;
  TH2F* btagLightPercenth;
  TH2F* btagPercenth;
  
  TH2F* sysPercenth;

  TH2F* uncPercenth;

  Float_t sig = 0.;
  Float_t stat = 0.;

  Float_t jes = 0.;
  Float_t bc = 0.;
  Float_t light = 0.;
  
  Float_t unc = 0.;

  int N = bkgTree->GetEntries();

  for ( int ibr = 0; ibr < NBR; ibr++){
    dirname = ""; dirname += BR[ibr];
    outFile->mkdir(dirname);
    outBRDir = outFile->GetDirectory(dirname);
//.........这里部分代码省略.........

void recurseFile(TDirectory *indir, TDirectory *outdir,
                 double etawid, double etamid) {

  TDirectory *curdir = gDirectory;

  // Automatically go through the list of keys (directories)
  TList *keys = indir->GetListOfKeys();
  TListIter itkey(keys);
  TObject *key, *obj;
  TDirectory *dir;

  while ( (key = itkey.Next()) ) {

    if (_debug) cout << key->GetName() << endl << flush;
    obj = ((TKey*)key)->ReadObj(); assert(obj);
    dir = indir;

    // Found a subdirectory: copy it to output and go deeper
    if (obj->InheritsFrom("TDirectory")) {

      //assert(outdir->mkdir(obj->GetName()));
      outdir->mkdir(obj->GetName());
      assert(outdir->cd(obj->GetName()));
      TDirectory *outdir2 = outdir->GetDirectory(obj->GetName()); assert(outdir2);
      outdir2->cd();

      assert(indir->cd(obj->GetName()));
      TDirectory *indir2 = indir->GetDirectory(obj->GetName());
      indir2->cd();

      // Check if directory name contains information on eta bin width
      float etamin, etamax;
      if ( (sscanf(indir->GetName(),"Eta_%f-%f",&etamin,&etamax)==2)
           && (etamax>etamin) ) {
        etawid = 2.*(etamax-etamin);
        etamid = 0.5*(etamax+etamin);
        //cout << "Eta bin width: " << etawid << flush << endl;
      }

      recurseFile(indir2, outdir2, etawid, etamid);
      //outdir2->Write(); // does this speedup or slow down?
    } // inherits from TDirectory

    // Found a plot: normalize if hpt, then copy to output
    if (obj->InheritsFrom("TH1")) {

      outdir->cd();
      TObject *obj2 = obj->Clone(obj->GetName());

      // Normalize hpt and hselpt histograms
      // Same for hbpt
      if (string(obj2->GetName())=="hpt" ||
          string(obj2->GetName())=="hpt_evt" ||
          string(obj2->GetName())=="hpt_jet" ||
          string(obj2->GetName())=="hpt_pre" ||
          string(obj2->GetName())=="hpt0" ||
          string(obj2->GetName())=="hpt1" ||
          string(obj2->GetName())=="hpt2" ||
          string(obj2->GetName())=="hpt3" ||
          string(obj2->GetName())=="hpt_jk1" ||
          string(obj2->GetName())=="hpt_jk2" ||
          string(obj2->GetName())=="hpt_jk3" ||
          string(obj2->GetName())=="hpt_jk4" ||
          string(obj2->GetName())=="hpt_jk5" ||
          string(obj2->GetName())=="hpt_jk6" ||
          string(obj2->GetName())=="hpt_jk7" ||
          string(obj2->GetName())=="hpt_jk8" ||
          string(obj2->GetName())=="hpt_jk9" ||
          string(obj2->GetName())=="hpt_jk10" ||
          string(obj2->GetName())=="hpt_l1off" ||
          string(obj2->GetName())=="hpt_l1fast" ||
          string(obj2->GetName())=="hpt_plus" ||
          string(obj2->GetName())=="hpt_minus" ||
          string(obj2->GetName())=="hpt0_plus" ||
          string(obj2->GetName())=="hpt0_minus" ||
          string(obj2->GetName())=="hpt_noid" ||
          string(obj2->GetName())=="hpt_noevtid" ||
          string(obj2->GetName())=="hpt_nojetid" ||
          string(obj2->GetName())=="hpt_ak5calo" ||
          string(obj2->GetName())=="hpt_ak5pf" ||
          string(obj2->GetName())=="hpt_evt_ak5pf" ||
          string(obj2->GetName())=="hpt_jet_ak5pf" ||
          string(obj2->GetName())=="hselpt" ||
          string(obj2->GetName())=="hpt_r" ||
          string(obj2->GetName())=="hpt_g" ||
          string(obj2->GetName())=="hpt_gg" ||
          string(obj2->GetName())=="hpt_g0" ||
          string(obj2->GetName())=="hpt_g0tw" ||
          string(obj2->GetName())=="hdjmass" ||
          string(obj2->GetName())=="hdjmass0" ||
          string(obj2->GetName())=="hdjmass0_hgg") {

          cout << "." << flush;

        TH1D *hpt = (TH1D*)obj2;
        bool isgen = TString(obj2->GetName()).Contains("pt_g");
        bool isoth = (TString(obj2->GetName()).Contains("pt_no") ||
                      TString(obj2->GetName()).Contains("djmass") ||
                      TString(obj2->GetName()).Contains("hpt0") ||
                      TString(obj2->GetName()).Contains("l1off") ||
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  TH1D* jet2h[NDirs];
  TH1D* jet3h[NDirs];
  TH1D* jet4h[NDirs];

  TH1D* nbjetsh[NDirs];
  TH1D* bjet1h[NDirs];
  TH1D* bjetHighDh[NDirs];
  TH1D* bdiscHh[NDirs];
  
  TH1D* hth[NDirs];
  TH1D* ht3h[NDirs];
  TH1D* ht4h[NDirs];
  TH1D* ht5h[NDirs];
  TH1D* htratioh[NDirs];
  TH1D* meth[NDirs];
  TH1D* meffh[NDirs];
  TH1D* yh[NDirs];

  TH1D* mth[NDirs];
  TH1D* mlb1h[NDirs];
  TH1D* mlbh[NDirs];
  TH1D* m3bh[NDirs];
  TH1D* m3h[NDirs];
  TH1D* centralityh[NDirs];
  TH1D* mt2wh[NDirs];
  TH1D* hadchi2h[NDirs];
  TH1D* topnessh[NDirs];

  TH1D* dphiminh[NDirs];
  TH1D* drlb1h[NDirs];
  TH1D* drlbminh[NDirs];

  TString dirName = ""; dirName += iSR;
  if (outFile->GetDirectory(dirName)) outFile->Delete(dirName + ";*");
  outFile->mkdir( dirName);
  TDirectory* SRDir = outFile->GetDirectory( dirName);
  SRDir->cd();
  
  for ( int iDir = 0; iDir < NDirs; iDir++){    
    SRDir->mkdir( controlDirName[iDir]);
    controlDir[iDir] = SRDir->GetDirectory( controlDirName[iDir]);
    controlDir[iDir]->cd();

    npvh[iDir] = new TH1D( "npv", "NPV", 51, -0.5, 50.5);
    ngoodpvh[iDir] = new TH1D( "ngoodpv", "NgoodPV", 51, -0.5, 50.5);

    lpth[iDir]     = new TH1D( "lpt", "lep p_{T} [GeV]",  12, 25., 500.);
    letah[iDir]    = new TH1D( "leta", "lep #Eta", 30, -3., 3.);
    lrelisoh[iDir] = new TH1D( "lRelIso", "lep RelIso", 30, 0., 1.);
    
    njetsh[iDir] = new TH1D(  "njets", "jets multiplicity",    10, -0.5, 9.5);
    jet1h[iDir]  = new TH1D(  "jet1", "1st jet p_{T} [GeV]",   25, 0., 500.);
    jet2h[iDir]  = new TH1D(  "jet2", "2nd jet p_{T} [GeV]",   25, 0., 500.);
    jet3h[iDir]  = new TH1D(  "jet3", "3rd jet p_{T} [GeV]",   25, 0., 500.);
    jet4h[iDir]  = new TH1D(  "jet4", "4th jet p_{T} [GeV]",   25, 0., 500.);
    
    nbjetsh[iDir]    = new TH1D( "nbjets", "b jets multiplicity",  6, -0.5, 5.5);  
    bjet1h[iDir]     = new TH1D(     "bjet1", "Leading b jet p_{T} [GeV]", 25, 0., 500.);
    bjetHighDh[iDir] = new TH1D( "bjetHighD", "p_{T} of the highest b disc jet [GeV]",  12, 0., 500.);
    bdiscHh[iDir]    = new TH1D( "bdisc", "bdisc",  20, 0., 1.);  

    meth[iDir] = new TH1D( "MET", "MET [GeV]", 12, 0., 400.);
    
    hth[iDir]      = new TH1D(   "Ht",   "Ht [GeV]",  20, 0., 1000.);
    ht3h[iDir]     = new TH1D(   "Ht3",   "Ht3 [GeV]",  20, 0., 1000.);
    ht4h[iDir]     = new TH1D(   "Ht4",   "Ht4 [GeV]",  20, 0., 1000.);

int main(int ac, char *av[]) { 
  try {
    double fMin, fMax;
    string ext;
    po::options_description desc("Allowed options");
    desc.add_options()
      ("help", "produce help message")
      ("include-path,I", po::value< vector<string> >(), 
       "include path")
      ("input-file", po::value< vector<string> >(), "input file")
      ("min,m", po::value<double>(&fMin)->default_value(60), "minimum value for fit range")
      ("max,M", po::value<double>(&fMax)->default_value(120), "maximum value for fit range")
      ("convbwintgamg", 
       "fit to the convolution of a breit-wigner plus interference term and gamma propagator and a gaussian")
      ("convexpbwintgamg", 
       "fit to the convolution of the product between an exponential and a breit-wigner plus interference term and gamma propagator with a gaussian")
      ("convbwintgam2gf", 
       "fit to the convolution of the breit-wigner plus interference term and gamma propagator with a linear combination of fixed gaussians")
      ("output-file,O", po::value<string>(&ext)->default_value(".ps"), 
       "output file format")
      ;
    
    po::positional_options_description p;
    p.add("input-file", -1);
    
    po::variables_map vm;
    po::store(po::command_line_parser(ac, av).
	    options(desc).positional(p).run(), vm);
    po::notify(vm);
    
    if (vm.count("help")) {
      cout << "Usage: options_description [options]\n";
      cout << desc;
      return 0;
      }
    
    if (vm.count("include-path")) {
      cout << "Include paths are: " 
	   << vm["include-path"].as< vector<string> >() << "\n";
    }
    
    vector<string> v_file;
    vector<TH1D*> v_ZMassHistos;
    vector<string> v_eps;
    
    if (vm.count("input-file"))	{
      cout << "Input files are: " 
	   << vm["input-file"].as< vector<string> >() << "\n";
      v_file = vm["input-file"].as< vector<string> >();
      for(vector<string>::const_iterator it = v_file.begin(); 
	  it != v_file.end(); ++it) { 
	TFile * root_file = new TFile(it->c_str(),"read");
	TDirectory *Histos = (TDirectory*) root_file->GetDirectory("ZHisto");
	TDirectory *RecoHistos = (TDirectory*) Histos->GetDirectory("ZRecoHisto");
	TH1D * zMass = (TH1D*) RecoHistos->Get("ZMass");
	zMass->Rebin(4); //remember...
	zMass->GetXaxis()->SetTitle("#mu #mu invariant mass (GeV/c^{2})");
	v_ZMassHistos.push_back(zMass);
	gROOT->SetStyle("Plain");
	//gStyle->SetOptFit(1111);
	string f_string = *it;
	replace(f_string.begin(), f_string.end(), '.', '_');
	string eps_string = f_string + ext;
	v_eps.push_back(eps_string);
	cout << ">>> histogram loaded\n";
      }
      cout << v_file.size() << ", " << v_ZMassHistos.size() << ", " << v_eps.size() << endl;
      cout <<">>> Input files loaded\n";
    }

    IntegratorConv integratorConv(1.e-5);
    //Values for Z mass and width
    funct::Parameter mass("Mass", 91.364);
    funct::Parameter gamma("Gamma", 4.11);
    //Parameters for Z Line Shape
    funct::Parameter f_gamma("Photon factor", 0.838);
    funct::Parameter f_int("Interference factor", -0.00197);
    //Parameters for fits with gaussians
    funct::Parameter yield("Yield", 283000);
    funct::Parameter alpha("Alpha", 0.771); //the first gaussian is narrow
    funct::Parameter mean("Mean", 0); //0.229
    funct::Parameter sigma1("Sigma 1", 0.76); 
    funct::Parameter sigma2("Sigma 2", 2.94);
    //Parameter for exponential
    funct::Parameter lambda("Lambda", 0);
    
    if (vm.count("convbwintgamg")) {
      cout << "Fitting histograms in input files to the convolution of the Breit-Wigner plus Z/photon interference and photon propagator with a Gaussian\n";
      cout << ">>> set pars: " << endl;
      cout << yield << endl; 
      cout << mass  << endl; 
      cout << gamma << endl; 
      cout << f_gamma << endl; 
      cout << f_int << endl; 
      cout << mean  << endl; 
      cout << sigma1 << endl;
      for(unsigned int i = 0; i < v_ZMassHistos.size(); ++i) { 
	TH1D * zMass = v_ZMassHistos[i]; 
	funct::ZLineShape zls(mass, gamma, f_gamma, f_int);
	funct::Gaussian gauss(mean, sigma1);
//.........这里部分代码省略.........