C++ TFile::Delete方法代码示例

void combineHistos(std::string type) {

    bool _dt = (type == "DATA");
    bool _mc = !_dt;

    TDirectory *curdir = gDirectory;

    // Input file: normalized pT spectra
    TFile *fin = 
        new TFile(Form("outputs/output-%s-2a.root", type.c_str()), "READ");
    assert(fin && !fin->IsZombie());
    
    // Top-level directory
    _top = gDirectory;

    // Output file: combined spectra 
    TFile *fout =
        new TFile(Form("outputs/output-%s-2b.root", type.c_str()), "RECREATE");
    assert(fout && !fout->IsZombie());


    std::cout << "Calling combineHistos(" << type << ");" << std::endl;
    std::cout << "Input file " << fin->GetName() << std::endl;
    std::cout << "Output file " << fout->GetName() << std::endl;
    std::cout << "Starting recursions. These may take a few seconds" << std::endl;

    // Store pT ranges to a nice map
    for (int itrg = 0; itrg != _jp_ntrigger; ++itrg) {

        std::string name = _jp_triggers[itrg];
        double lower = _jp_trigranges[itrg][0];
        double upper = _jp_trigranges[itrg][1];

        _ptranges[name] = pair<double, double>(lower, upper);
    }


    // Loop over all the directories recursively
    recurseFile(fin, fout, "hpt");
    recurseFile(fin, fout, "hpt_pre");
    if (_dt)
        recurseFile(fin, fout, "hlumi");
    if (_mc)
        recurseFile(fin, fout, "hpt_g0tw");

    curdir->cd();

    std::cout << std::endl << "Output stored in " << fout->GetName() << std::endl;
    
    // Close files
    fout->Close();
    fout->Delete();
    std::cout << "Output file closed" << std::endl;

    fin->Close();
    fin->Delete();
    std::cout << "Input file closed" << std::endl;
}

void normalizeHistos(string type) {

  assert(type==_jp_type);
  _nh_mc = (type=="MC" || type=="HW");
  _nh_dt = (type=="DATA");
  assert((_nh_mc || _nh_dt) && !(_nh_mc && _nh_dt));

  TFile *fin = new TFile(Form("output-%s-1.root",type.c_str()),"READ");
  assert(fin && !fin->IsZombie());

  TFile *fout = new TFile(Form("output-%s-2a.root",type.c_str()),"RECREATE");
  assert(fout && !fout->IsZombie());

  if (_lumiscale!=1 && !_nh_mc)
    cout << "Attention! : Scaling luminosity to the new estimate"
         << " by multiplying with " << _lumiscale << endl;

  if (_jp_usetriglumi) {
    cout << "Reading trigger luminosity from settings.h" << endl;
    for (int i = 0; i != _jp_ntrigger; ++i) {
      double lumi = _jp_triglumi[i]/1e6; // /ub to /pb
      cout << Form(" *%s: %1.3f /pb", _jp_triggers[i].c_str(),lumi) << endl;
      triglumi[_jp_triggers[i]] = lumi;
    }
  }

  cout << "Calling normalizeHistos("<<type<<");" << endl;
  cout << "Input file " << fin->GetName() << endl;
  cout << "Output file " << fout->GetName() << endl;
  cout << "Starting recursive loop. This may take a minute" << endl << flush;

  // Loop over all the directories recursively
  recurseFile(fin, fout);

  cout << endl;
  cout << "Recursive loop done." << endl;
  cout << "Writing output to " << fout->GetName() << endl;
  cout << "This may again take a minute" << endl << flush;
  fout->Write();
  cout << "Output written in " << fout->GetName() << endl;
  fout->Close();
  cout << "Output file closed" << endl;
  fout->Delete();
  cout << "Output file pointer deleted" << endl << flush;


  fin->Close();
  fin->Delete();

} // normalizeHistos

void func(const char * name) {

   TFile *g = new TFile("temp.root","UPDATE");

   // delete any pre-existing trees of the same name

   char nameCycle[100];
   sprintf(nameCycle,"%s;*",name);
   g->Delete(nameCycle);

   TTree *tree1 = new TTree(name,name);

   Float_t flt = 1;
   tree1->Branch("branch",&flt,"temp/F");

   for(int i; i < 500; i++) {
      tree1->Fill();
   }

   g->Write();
   tree1->Reset();
   delete tree1;
   tree1 = 0;

   g->Close();
   delete g;
   g = 0;

}

void make_scan_results()
{
  TFile *f = TFile::Open("scan_results.root", "UPDATE");

  f->Delete("SR;*");

  T = new TTree("SR", "Scanning results");

  TClonesArray* ts = new TClonesArray("IlcESDtrack", 32);
  TBranch * tb = T->Branch("T", &ts);
  delete ts;

  IlcMultiplicity *ms = 0;
  TBranch *mb = T->Branch("M", &ms);

  for (Int_t v = 0; v < 3; ++v)
  {
    vvv[v].vert   = 0;
    vvv[v].branch = T->Branch(vvv[v].bname, &vvv[v].vert);
  }

  for (Int_t i=0; i<=9999; ++i)
  {
    TString name;

    name.Form("Tracks_%04d", i);
    ts = (TClonesArray*) f->Get(name);
    if (ts == 0)
      continue;

    name.Form("Tracklets_%04d", i);
    ms = (IlcMultiplicity*) f->Get(name);
    if (ms == 0)
      Error("make_scan_results", "'%s' not found.", name.Data());

    tb->SetAddress(&ts);
    mb->SetAddress(&ms);

    for (Int_t v = 0; v < 3; ++v)
    {
      name.Form("%s_%04d", vvv[v].oname, i);
      vvv[v].vert = (IlcESDVertex*) f->Get(name);
      if (vvv[v].vert == 0)
        Error("make_scan_results", "'%s' not found.", name.Data());
      vvv[v].branch->SetAddress(&vvv[v].vert);
    }

    T->Fill();

    delete ts;
    delete ms;
    for (Int_t v = 0; v < 3; ++v) delete vvv[v].vert;
  }

  T->Write();

  f->Close();
  delete f;
}

void ChangeDataTreeName(std::string suffix){

  // Define names
  TString filename = "", treename = "", histname = "";
  filename.Form("central_%s.root",suffix.c_str());
  treename.Form("id_%s",suffix.c_str());
  histname.Form("sumOfMcWeights_%s",suffix.c_str());

  TString _pathFile(string(getenv("WORKAREA")) + "/histoAna" + "/SusyAna/histos_110312_13fb_n0111_DD_MMtrial9_SYS_HFT/HFTOutputs_merge/");
  filename = _pathFile+filename;
  cout << "open " << filename << endl;

  // Open file
  TFile* file         = new TFile(filename,"UPDATE");
  if(file->IsOpen()==kFALSE){
    cerr << "Cannot open file " << filename << endl;
    abort();
  }
  file->ls();
  TTree* treeOriginal = (TTree*) file->Get(treename);
  TH1D*  histOriginal = (TH1D*)  file->Get(histname);
  TTree* treeNew      = (TTree*) treeOriginal->Clone();
  TH1D*  histNew      = (TH1D*)  histOriginal->Clone();

  // Set-Write-Delete
  histNew->SetName ("sumOfMcWeights_Data" );
  histNew->SetTitle("sumOfMcWeights_Data" );
  histNew->Write   ("",TObject::kOverwrite);
  treeNew->SetName ("id_Data"             );
  treeNew->SetTitle("id_Data"             );
  treeNew->Write   ("",TObject::kOverwrite);
  file->Delete     (histname+";1"         );
  file->Delete     (treename+";1"         );

  // Close file
  file->Close();
}

void testweight::Loop(bool cs)
{
//   In a ROOT session, you can do:
//      Root > .L testweight.C
//      Root > testweight t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   TFile * hOutputFile   = new TFile("spectrumrew.root", "RECREATE" ) ;
   TH1D* pt_rew =  new TH1D("pt_rew","pt_rew",50,20,100);
   TH1D* pt_norew =   new TH1D("pt_norew","pt_norew",50,20,100);
   TH1D* mass_rew =  new TH1D("mass_rew","mass_rew",40,100,180);
   TH1D* mass_norew =  new TH1D("mass_norew","mass_norew",40,100,180);
   TH2D* pt2d_norew = new TH2D("pt2d_norew","pt2d_norew",15,25,100,15,55,160);
   TH2D* pt2d_rew = new TH2D("pt2d_rew","pt2d_rew",15,25,100,15,55,160);

   getweights();

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      
      if(massggnewvtx<90 || massggnewvtx>190) continue;
      //if(massggnewvtx<100 || massggnewvtx>180) continue;
      
      if((TMath::Abs(etascphot1)>1.4442&&TMath::Abs(etascphot1)<1.566)||(TMath::Abs(etascphot2)>1.4442&&TMath::Abs(etascphot2)<1.566)
	 || TMath::Abs(etascphot1)>2.5 || TMath::Abs(etascphot2)>2.5) continue;  // acceptance
      
      if(ptphot1<55) continue; //pt first photon
      if(ptphot2<25) continue; //pt second photon
      
      //      if(ptcorrjet1<30 || TMath::Abs(etajet1)>4.7) continue; //pt first jet

      bool idphot1 = (idcicphot1 >= 4);
      bool idphot2 = (idcicphot2 >= 4);

      if(!cs){ // photon id no control sample

	if(!(idphot1)) continue;
	if(!(idphot2)) continue;
  
      }else{ // photon id for control sample
	
	if( !( (idphot1 && !idphot2 && !pid_hasMatchedPromptElephot2) || (idphot2 && !idphot1 && !pid_hasMatchedPromptElephot1) ) ) continue; 

      }

      pt_norew->Fill(ptphot2);
      mass_norew->Fill(massggnewvtx);

      double minptsublead(25), maxptsublead(100);
      double minptlead(55), maxptlead(160);
      double sizex = (maxptsublead - minptsublead)/15.;
      double sizey = (maxptlead - minptlead)/15.;
      int i = int((ptphot2-25)/sizex);
      int j = int((ptphot1-55)/sizey);
      if(i<0) i=0;
      if(j>14) j=14;
      if(i<0) i=0;
      if(j>14) j=14;
      if(i>-1 && i<15 && j>-1 && j<15){
	pt_rew->Fill(ptphot2,weights_[i][j]);
	mass_rew->Fill(massggnewvtx,weights_[i][j]);
	pt2d_norew->Fill(ptphot2,ptphot1,1);
	pt2d_rew->Fill(ptphot2,ptphot1,weights_[i][j]);
      }
      //      else  pt_rew->Fill(ptphot2,1.);
   }
   
   hOutputFile->Write() ;
   hOutputFile->Close() ;
   hOutputFile->Delete();
     
}

//.........这里部分代码省略.........

  TFile *fout = new TFile(Form("output-%s-2c.root",type.c_str()),"RECREATE");
  assert(fout && !fout->IsZombie());

  // Select top category
  assert(fin->cd("Standard"));
  fin->cd("Standard");
  TDirectory *din0 = gDirectory;

  assert(fmc->cd("Standard"));
  fmc->cd("Standard");
  TDirectory *dmc0 = gDirectory;

  fout->mkdir("Standard");
  assert(fout->cd("Standard"));
  fout->cd("Standard");
  TDirectory *dout0 = gDirectory;

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

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

    obj = ((TKey*)key)->ReadObj(); assert(obj);
    
    // Found a subdirectory
    if (obj->InheritsFrom("TDirectory")
	&& string(obj->GetName())!="Eta_0.0-1.3"
	&& string(obj->GetName())!="Eta_3.0-3.2"
	&& string(obj->GetName())!="Eta_3.2-4.7") {

      assert(din0->cd(obj->GetName()));
      din0->cd(obj->GetName());
      TDirectory *din = gDirectory;

      assert(dmc0->cd(obj->GetName()));
      dmc0->cd(obj->GetName());
      TDirectory *dmc = gDirectory;
      /*
      assert(fnlo_cteq->cd());
      TDirectory *dnlo_cteq = gDirectory;
      assert(fnlo_ct10->cd());
      TDirectory *dnlo_ct10 = gDirectory;
      assert(fnlo_mstw->cd());
      TDirectory *dnlo_mstw = gDirectory;
      assert(fnlo_nnpdf->cd());
      TDirectory *dnlo_nnpdf = gDirectory;
      assert(fnlo_hera->cd());
      TDirectory *dnlo_hera = gDirectory;
      //
      assert(fnlo_as->cd());
      TDirectory *das = gDirectory;

      assert(fnp->cd());
      TDirectory *dnp = gDirectory;
      */

      dout0->mkdir(obj->GetName());
      assert(dout0->cd(obj->GetName()));
      dout0->cd(obj->GetName());
      TDirectory *dout = gDirectory;
      
      // Process subdirectory
      //theoryBin(din, dnlo_cteq, dnp, dnlo_ct10, dnlo_mstw, dnlo_nnpdf,
      //	dnlo_hera,
      //	das, dout);
      theoryBin(din, dmc, dout);
      //dataBin(din, dout);
    } // inherits TDirectory
  } // while

  cout << "Output stored in " << fout->GetName() << endl;
  fout->Write();
  fout->Close();
  fout->Delete();

  fin->Close();
  fin->Delete();

  fmc->Close();
  /*
  fnlo_cteq->Close();
  fnlo_cteq->Delete();
  fnlo_ct10->Close();
  fnlo_ct10->Delete();
  fnlo_mstw->Close();
  fnlo_mstw->Delete();
  fnlo_nnpdf->Close();
  fnlo_nnpdf->Delete();
  fnlo_as->Close();
  fnlo_as->Delete();
  */

  // Empty work space so that next steps won't get messed up
  //if (gROOT->IsBatch()) {
  //gROOT->Clear();
  //}
} // theory

//.........这里部分代码省略.........
 
                    
          int label = cl->GetLabel(0);
          TParticle* part = 0;
          if (label>=0 && (part=stack->Particle(label)) ) {
            cSum.pdg = part->GetPdgCode();
            cSum.eta = part->Eta();
            cSum.pt  = part->Pt();
            cSum.phi = part->Phi();
            cSum.prim = stack->IsPhysicalPrimary(label);
          } 
          cSum.ntr = 0;
          for (int ilb=0;ilb<3;ilb++) if (cl->GetLabel(ilb)>=0) cSum.ntr++;
          for (int i=0;i<3;i++) cSum.xyz[i] = xyzClGloF[i];
          //
          trOut->Fill();
          /*
          if (clsize==5) {
            printf("\nL%d(%c) Mod%d, Cl:%d | %+5.1f %+5.1f (%d/%d)|H:%e %e %e | C:%e %e %e\n",ilr,cl->TestBit(kSplit) ? 'S':'N',
             modID,icl,(txyzH[0]-xyzClTr[0])*1e4,(txyzH[2]-xyzClTr[2])*1e4, row,col,
             gxyzH[0],gxyzH[1],gxyzH[2],xyzClGlo[0],xyzClGlo[1],xyzClGlo[2]);
            cl->Print();
            pHit->Print();
            //
            double a0,b0,c0,a1,b1,c1,e0;
            pHit->GetPositionL0(a0,b0,c0,e0);
            pHit->GetPositionL(a1,b1,c1);
            float cloc[3];
            cl->GetLocalXYZ(cloc);
            printf("LocH: %e %e %e | %e %e %e\n",a0,b0,c0,a1,b1,c1);
            printf("LocC: %e %e %e | %e %e %e\n",cloc[0],cloc[1],cloc[2],xyzClTr[0],xyzClTr[1],xyzClTr[2]);
          }
          */
          //
        }
      }
    }
    
    //    layerClus.Clear();
    //
    arrMCTracks.Delete();
  }//event loop
  //
  DB.EndAndSort();
  DB.SetThresholdCumulative(0.95);
  cout << "Over threshold: : "<< DB.GetOverThr()<<endl;
  DB.Grouping(10,10);
  DB.PrintDB("Database1.txt"); 
  flOut->cd();
  trOut->Write();
  delete trOut;
  flOut->Close();
  flOut->Delete();
  DrawReport("clinfo.ps",&histoArr);
  TFile* flDB = TFile::Open("TopologyDatabase.root", "recreate");
  flDB->WriteObject(&DB,"DB","kSingleKey");
  flDB->Close();
  delete flDB;

  TCanvas* cnv123 = new TCanvas("cnv123","cnv123");
  cnv123->Divide(1,2);
  cnv123->Print("anglesdistr.pdf[");
  cnv123->cd(1);
  hL0A->Draw();
  cnv123->cd(2);
  hL0B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL1A->Draw();
  cnv123->cd(2);
  hL1B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL2A->Draw();
  cnv123->cd(2);
  hL2B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL3A->Draw();
  cnv123->cd(2);
  hL3B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL4A->Draw();
  cnv123->cd(2);
  hL4B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL5A->Draw();
  cnv123->cd(2);
  hL5B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->cd(1);
  hL6A->Draw();
  cnv123->cd(2);
  hL6B->Draw();
  cnv123->Print("anglesdistr.pdf");
  cnv123->Print("anglesdistr.pdf]");
  //
}

//.........这里部分代码省略.........
    if (faketype == 11) {
      fakeTypeString = "Electron";
      maxY = 0.1;
    } else if (faketype == 13) {
      fakeTypeString = "Muon";
      maxY = 0.01;
    } else {
      cerr << "Error: faketype = " << faketype << " not recognized.\n";
      assert(false);
    }
    
    //*****************************************************************************************
    //Load input histogram files
    //*****************************************************************************************
    madgraphDatasetFiles.clear(); 
    madgraphDatasetNames.clear();
    pythiaDatasetFiles.clear();
    pythiaDatasetNames.clear();
    
    //Madgraph WJets
    madgraphDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_f8-wjets-mg-id9.root");
    madgraphDatasetNames.push_back("f8-wjets-mg-id9");
    //Pythia WJets
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-we-id9.root");
    pythiaDatasetNames.push_back("s8-we-id9");
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-wm-id9.root");
    pythiaDatasetNames.push_back("s8-wm-id9");
    pythiaDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-wtau-id9.root");
    pythiaDatasetNames.push_back("s8-wtau-id9");

    //*****************************************************************************************
    //Construct the fake rate Histgrams. Fit with function.
    //*****************************************************************************************
    vector<TH1F*> fakeRateHistograms; fakeRateHistograms.clear();
    vector<TF1*> fakeRateFits; fakeRateFits.clear();
    vector<string> fakeRateLegendNames; fakeRateLegendNames.clear();

    //pythia version
    TH1F *FakeRatePt_pythia = createFakeRatePtHist(pythiaDatasetFiles, pythiaDatasetNames, 1, faketype, chargetype, eventType,
                                                   (chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominatorFakeRate_Pt_PythiaWJets").c_str());
    FakeRatePt_pythia->GetYaxis()->SetTitle((chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominator Fake Rate").c_str());    
    TF1 *FakeRateFitFunction_pythia = NULL;
    if (fitFakeRate) {
      FakeRateFitFunction_pythia = fitFakeRatePtHist(FakeRatePt_pythia);
      FakeRateFitFunction_pythia->SetName((chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominatorFakeRateFunction_Pt_PythiaWJets").c_str());
      fakeRateFits.push_back(FakeRateFitFunction_pythia);
    }
    fakeRateHistograms.push_back(FakeRatePt_pythia);
    fakeRateLegendNames.push_back((chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominator Pythia W+Jets").c_str());

    //madgraph version
//     TH1F *FakeRatePt_madgraph = createFakeRatePtHist(
//       madgraphDatasetFiles, madgraphDatasetNames,1,faketype, chargetype,eventType,
//       (chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominator FakeRate_Pt_MadgraphWJets").c_str());
//     FakeRatePt_madgraph->GetYaxis()->SetTitle((chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominator Fake Rate").c_str());    

//     TF1 *FakeRateFitFunction_madgraph = NULL;
//     if (fitFakeRate) {
//       FakeRateFitFunction_madgraph = fitFakeRatePtHist(FakeRatePt_madgraph);
//       FakeRateFitFunction_madgraph->SetName(
//       (chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominatorFakeRateFunction_Pt_MadgraphWJets").c_str());
//     }
//     if (fitFakeRate) fakeRateFits.push_back(FakeRateFitFunction_madgraph);
//     fakeRateHistograms.push_back(FakeRatePt_madgraph);
//     fakeRateLegendNames.push_back(chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominator Madgraph W+Jets");
    
    //*****************************************************************************************
    //Draw the plots
    //*****************************************************************************************
    drawFakeRatePlots(fakeRateHistograms,fakeRateFits,fakeRateLegendNames,
                      (chargeTypeName + fakeTypeString + eventTypeName + "TrackDenominatorFakeRatePt").c_str(),maxY);

    //*****************************************************************************************
    //save them to the output file
    //*****************************************************************************************
    for (UInt_t i=0; i<fakeRateHistograms.size();i++) {
      //write the histograms. delete any existing one with the same name
      if (file->FindObjectAny(fakeRateHistograms[i]->GetName())) {
        string tmpname = fakeRateHistograms[i]->GetName();        
        file->Delete((tmpname+";*").c_str());
      }

      fakeRateHistograms[i]->SetDirectory(file);
      file->WriteTObject(fakeRateHistograms[i]);

      //write the fit functions. delete any existing one with the same name
      if (fitFakeRate) {
        if (file->FindObjectAny(fakeRateFits[i]->GetName())) {
          string tmpname = fakeRateFits[i]->GetName();        
          file->Delete((tmpname+";*").c_str());
        }
        file->WriteTObject(fakeRateFits[i]);
      }
    }
  }

  file->Close();
  return;

}

/*************************************************************************************
 * createMFOutfile: moves the MLB distributions into an output file for use in 
 *                  R. Nally's MassFit.C code.
 *  input:  the main() arguments array
 *  output: writes to an output file the histograms, in a MassFit.C-readable format 
 *
 *  Structure-wise: can be implemented into class easily.
 ***********************************/
void createMFOutfile(const char* argv[]) {
  TFile *f = new TFile(argv[2]);

  //create the output file we'd like to write histograms to
  TFile *output = new TFile(outfileName, "RECREATE");
  
  //get the filesystem information from the file
  f->cd();
  TList *alokDirs = (TList*) f->GetListOfKeys();

  //create a list of "All" histograms and initialize (TODO)
  TList *allHists = new TList;

  //loop through the directories in the input file
  for(int idir=0; alokDirs->At(idir-1) != alokDirs->Last(); idir++) {
    //if it's not mlb, we don't care
    if(!TString(alokDirs->At(idir)->GetName()).Contains("_Mlb")) continue;
    
    //get the file directory information, its histograms
    TDirectory *cDir  = (TDirectory*) f->Get(alokDirs->At(idir)->GetName());
    TList *alokHistos = (TList*)      cDir->GetListOfKeys();

    //loop through the histograms in the current directory 
    for(int ihisto=0; alokHistos->At(ihisto-1) != alokHistos->Last(); ihisto++) {
      // don't consider the graph objects
      if(TString(alokHistos->At(ihisto)->GetName()).Contains("Graph_")) continue; 

      // clone the histogram, give it a new name
      TString cloneName = formatName(alokHistos->At(ihisto)->GetName(),nominalWidth);
      TH1F *thisto = (TH1F*) cDir->Get(alokHistos->At(ihisto)->GetName());
      TH1F *tclone = (TH1F*) thisto->Clone(cloneName);

      // open the outfile and write the histo in
      output->cd();
      TList *outkeys = (TList*) output->GetListOfKeys();

      // if the histogram already exists, add thisto to the existing one
      // messy, but removes the need for magic
      for(int iout=0; outkeys->At(iout-1) != outkeys->Last(); iout++) {
        if(outkeys->At(iout)->GetName() == cloneName) {
          cout<<" - found another histogram in output with name "<<cloneName<<endl;
          TH1F *theHisto = (TH1F*) output->Get(cloneName);
          cout<<" - got the same histogram from the output file"<<endl;
          TH1F *tHclone  = (TH1F*) theHisto->Clone(cloneName);
          cout<<" - cloned the histogram"<<endl;

          cout<<" - adding in clone"<<endl;
          tclone->Add(tHclone);

          cout<<" - deleting the original histogram from the output file"<<endl;
          output->Delete(cloneName + TString(";1"));
          cout<<" - deleted thing from output file"<<endl;
          cout<<" - tclone looks like "<<tclone<<endl;
        }
      }
      cout<<" - writing the tclone to file"<<endl;
      tclone->Write();

      // reopen the input root file and clean up
      f->cd();
      delete thisto;
      delete tclone;
    }
  }

  f->Close();
  output->cd();
  output->Close();

  // if we want to interpolate, start making more histograms
  if(interpolate) {
    // Have to reopen the outfile because ROOT is odd
    TFile *output2 = new TFile(outfileName, "UPDATE");
    output2->cd();

    std::pair<double,TString> maxPair = *moreFiles.begin();
    double maxWidth = maxPair.first;

    // check that it makes sense to interpolate with our settings
    if(maxWidth <= nominalWidth) {
      cout<<"\n\nERROR: Max width is less than or equal to the nominal width! Exiting."<<endl;
      exit(EXIT_FAILURE);
    }

    // open the input file and loop through the relevant directories
    TFile *maxFile = new TFile(maxPair.second);
    for(int i=0; i<lepsSize; i++) {
      // change directory and get the name of the folder we want to access
      maxFile->cd();
      char a[128];
      sprintf(a, "mlbwa_%s_Mlb", leps[i]);
      TDirectory *tDir = (TDirectory*) maxFile->Get(a);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
  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.);

//.........这里部分代码省略.........
                // double Z_integrals[2*Nsamples];
                double W_integrals[2*Nsamples];

                //LOOP OVER THE SAMPLES
                for(int isample=0; isample<Nsamples;isample++){
                  for(int k=0;k<WMass::NFitVar;k++){

                    if(!TemplatesW_NonScaled[c][m][n][h][k][isample][ieta][jmass]) continue;
                    histcounter++;

                    // DECLARE NEW HISTOS
                    TH1D* Wtempl_NonScaled;

                    // To BE CHECKED
                    // W_histoname[isample] = samples_str[isample] == "DATA" ? Form("data_obs_W%s%s_%s",Wlike.Data(),WCharge_str[c].Data(),WMass::FitVar_str[k].Data()) : Form("W%s%s_%s_%s_pdf%d-%d%s",Wlike.Data(),WCharge_str[c].Data(),samples_str[isample].Data(),WMass::FitVar_str[k].Data(),WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets,h, RecoilCorrVarDiagoParU1orU2fromDATAorMC>0?Form("_RecoilCorrVar%d",m):"");
                    // Wtempl=new TH1D(W_histoname[isample],W_histoname[isample],TemplatesW[c][m][n][h][k][isample][ieta][jmass]->GetXaxis()->FindBin(xmax*fitrange_Scaling)-TemplatesW[c][m][n][h][k][isample][ieta][jmass]->GetXaxis()->FindBin(xmin*fitrange_Scaling),xmin*fitrange_Scaling,xmax*fitrange_Scaling);
                    W_histoname_NonScaled[isample] = samples_str[isample] == "DATA" ? Form("data_obs_W%s%s_%sNonScaled",Wlike.Data(),WCharge_str[c].Data(),WMass::FitVar_str[k].Data()) : Form("W%s%s_%s_%sNonScaled_pdf%d-%d%s%s",Wlike.Data(),WCharge_str[c].Data(),samples_str[isample].Data(),WMass::FitVar_str[k].Data(),WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets,h, RecoilCorrVarDiagoParU1orU2fromDATAorMC>0?Form("_RecoilCorrVar%d",m):"",WMass::KalmanNvariations>1?Form("_KalmanVar%d",n):"");
                    // Wtempl_NonScaled=new TH1D(W_histoname_NonScaled[isample],W_histoname_NonScaled[isample],50, WMass::fit_xmin[k],WMass::fit_xmax[k]);
                    Wtempl_NonScaled=(TH1D*)TemplatesW_NonScaled[c][m][n][h][k][isample][ieta][jmass]->Clone(W_histoname_NonScaled[isample]);
                    Wtempl_NonScaled->SetName(W_histoname_NonScaled[isample]);
                    cout << "W_histoname_NonScaled[isample]="<<W_histoname_NonScaled[isample]<<endl;
                    Wtempl_NonScaled->SetTitle(W_histoname_NonScaled[isample]);

                    Wtempl_NonScaled->Write();

                    int nspaces1 = 50 - W_histoname_NonScaled[isample].Length();
                    // outTXTfile << Wtempl->GetName();
                    outTXTfile << Wtempl_NonScaled->GetName();
                    for(int ispace=0;ispace<nspaces1;ispace++) outTXTfile << " ";
                    // W_integrals[isample] = Wtempl->Integral();
                    W_integrals[isample] = Wtempl_NonScaled->Integral();
                    outTXTfile << W_integrals[isample] << endl;

                    Wtempl_NonScaled->Delete();

                    // 2D
                    for(int k2=k+1;k2<WMass::NFitVar-1;k2++){
                      if(k==k2 || k==3) continue;
                      if(TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]){
                        
                        // DECLARE NEW HISTOS

                        W_histoname_NonScaled_2d[isample] = samples_str[isample] == "DATA" ? Form("data_obs_W%s%s_%svs%s",Wlike.Data(),WCharge_str[c].Data(),WMass::FitVar_str[k].Data(),WMass::FitVar_str[k2].Data()) : Form("W%s%s_%s_%svs%s_pdf%d-%d%s%s",Wlike.Data(),WCharge_str[c].Data(),samples_str[isample].Data(),WMass::FitVar_str[k].Data(),WMass::FitVar_str[k2].Data(),WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets,h, RecoilCorrVarDiagoParU1orU2fromDATAorMC>0?Form("_RecoilCorrVar%d",m):"",WMass::KalmanNvariations>1?Form("_KalmanVar%d",n):"");
                        cout << "W_histoname_NonScaled_2d[isample]= " << W_histoname_NonScaled_2d[isample] << endl;
                        // Wtempl_NonScaled_2d_unrolled=(TH2D*)TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->Clone(W_histoname_NonScaled_2d[isample]);
                        int binsx = TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsX();
                        int binsy = TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsY();
                        cout << "binsx= " << binsx << " binsy= " << binsy << " binsx*binsy= " << binsx*binsy << endl;
                        TH1D* Wtempl_NonScaled_2d_unrolled = new TH1D(W_histoname_NonScaled_2d[isample],W_histoname_NonScaled_2d[isample],(binsx*binsy), 0, (binsx*binsy) );
                        cout << "Wtempl_NonScaled_2d_unrolled->GetNbinsX()= " << Wtempl_NonScaled_2d_unrolled->GetNbinsX() << endl;
                        cout << "Wtempl_NonScaled_2d_unrolled->GetNbinsY()= " << Wtempl_NonScaled_2d_unrolled->GetNbinsX() << endl;
                        cout << "TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsX()= " << TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsX() << endl;
                        cout << "TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsY()= " << TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetNbinsX() << endl;
                        int new_bin=1;
                        for(int biny=1; biny<binsy+1; biny++){
                          for(int binx=1; binx<binsx+1; binx++){
                            Wtempl_NonScaled_2d_unrolled->SetBinContent(new_bin,TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->GetBinContent(binx,biny));
                            new_bin++;
                          }
                        }
                        // cout << "W_histoname_NonScaled_2d[isample]="<<W_histoname_NonScaled_2d[isample]<<endl;

                        Wtempl_NonScaled_2d_unrolled->Print();
                        Wtempl_NonScaled_2d_unrolled->Write();
                        TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->SetName(Form("%s_2d",W_histoname_NonScaled_2d[isample].Data()));
                        TemplatesW_NonScaled_2d[c][m][n][h][k][isample][ieta][jmass]->Write();

//.........这里部分代码省略.........
              h_deltaM_PDF[k][c]->Fill(deltaM[m][h][k][c]);
            }
          }
          leg1->Draw("same");
          c_summary->Write();

        }
      }
    }
    
    for(int k=0;k<3;k++){
      if(WMass::PDF_members>1){

        cout << endl;
        int usedpdf = WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets;
        cout << Form("PDF %d with %s: #DeltaM_{W, min}= %.0f MeV, #DeltaM_{W,max}= %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMmin[k][c], deltaMmax[k][c]) << endl;
        double denominator = 1;

        TCanvas*c1= new TCanvas(Form("c_deltaM_PDF_W%s_%s",WCharge_str[c].Data(),WMass::FitVar_str[k].Data()),Form("c_deltaM_PDF_W%s_%s",WCharge_str[c].Data(),WMass::FitVar_str[k].Data()));
        
        g_deltaM_PDF[k][c]->Draw("ape");
        TLatex *text;
        text = new TLatex(0.2,0.85,Form("W%s, PDF %d with %s: #DeltaM_{W, min}= %.0f MeV, #DeltaM_{W,max}= %.0f MeV", WCharge_str[c].Data(), WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMmin[k][c], deltaMmax[k][c]));
        text->SetTextSize(0.035);
        text->SetNDC();
        text->Draw();
        
        if(usedpdf==232000){ // NNPDF, see ref. http://arxiv.org/pdf/1301.6754v1.pdf formulae 2.23 and 2.25
          // cout << Form("PDF %d with %s: #DeltaM_{W} square summed = -%.0f MeV, #DeltaM_{W} square summed = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMnegSummed[k]/denominator, deltaMposSummed[k]/denominator) << endl;
          denominator = TMath::Sqrt(WMass::PDF_members/2-1);
          cout << Form("PDF %d with %s: #DeltaM_{W} square summed = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMSummed[k][c]/denominator) << endl;
          text = new TLatex(0.2,0.8,Form("PDF %d with %s: #DeltaM_{W} square summed = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMSummed[k][c]/denominator));
          text->SetTextSize(0.035);
          text->SetNDC();
          text->Draw();
        }else{
          cout << Form("PDF %d with %s: #DeltaM_{W} as of J. Rojo = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), (deltaMJuan[k][c]/2/denominator)) << endl;
          text = new TLatex(0.2,0.8,Form("PDF %d with %s: #DeltaM_{W} as of J. Rojo = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), (deltaMJuan[k][c]/2/denominator)));
          text->SetTextSize(0.035);
          text->SetNDC();
          text->Draw();
        }

        g_deltaM_PDF[k][c]->Write();
        h_deltaM_PDF[k][c]->Write();
        c1->Write();
        
        // mean = TMath::Mean(n, &v[0]);
        // vector<double> v;
        // // std::generate(v.begin(), v.end(), 1);
        // if(k==1)
          // for(int h=0; h<WMass::PDF_members; h++){
            // // v.push_back(deltaM[h][k]);
            // cout << deltaM[h][k] << endl;
          // }
        // double meanWmass= TMath::Mean(v.begin(), v.end());
        // cout << "meanWmass= " << meanWmass << endl;
      }
      if(WMass::NtoysMomCorr>1){
        cout << endl;
        // int usedpdf = WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets;
        cout << Form("MomCorr toys with %s: #DeltaM_{W, min}= %.0f MeV, #DeltaM_{W,max}= %.0f MeV", WMass::FitVar_str[k].Data(), deltaMmin[k][c], deltaMmax[k][c]) << endl;
        double denominator = 1;

        TCanvas*c1= new TCanvas(Form("c_deltaM_PDF_%s",WMass::FitVar_str[k].Data()),Form("c_deltaM_PDF_%s",WMass::FitVar_str[k].Data()));
        g_deltaM_PDF[k][c]->Draw("ape");
        TLatex *text;
        text = new TLatex(0.2,0.85,Form("MomCorr toys with %s: #DeltaM_{W, min}= %.0f MeV, #DeltaM_{W,max}= %.0f MeV", WMass::FitVar_str[k].Data(), deltaMmin[k][c], deltaMmax[k][c]));
        text->SetTextSize(0.035);
        text->SetNDC();
        text->Draw();
        
        // if(usedpdf==232000){ // NNPDF, see ref. http://arxiv.org/pdf/1301.6754v1.pdf formulae 2.23 and 2.25
          // cout << Form("PDF %d with %s: #DeltaM_{W} square summed = -%.0f MeV, #DeltaM_{W} square summed = %.0f MeV", WMass::PDF_sets<0?generated_PDF_set:WMass::PDF_sets, WMass::FitVar_str[k].Data(), deltaMnegSummed[k]/denominator, deltaMposSummed[k]/denominator) << endl;
          denominator = TMath::Sqrt(WMass::NtoysMomCorr/2-1);
          cout << Form("MomCorr toys with %s: #DeltaM_{W} square summed = %.0f MeV", WMass::FitVar_str[k].Data(), deltaMSummed[k][c]/denominator) << endl;
          text = new TLatex(0.2,0.8,Form("MomCorr toys with %s: #DeltaM_{W} square summed = %.0f MeV", WMass::FitVar_str[k].Data(), deltaMSummed[k][c]/denominator));
          text->SetTextSize(0.035);
          text->SetNDC();
          text->Draw();
        // }else{
          // cout << Form("MomCorr toys with %s: #DeltaM_{W} as of J. Rojo = %.0f MeV", WMass::FitVar_str[k].Data(), (deltaMJuan[k]/2/denominator)) << endl;
          // text = new TLatex(0.2,0.8,Form("MomCorr toys with %s: #DeltaM_{W} as of J. Rojo = %.0f MeV", WMass::FitVar_str[k].Data(), (deltaMJuan[k]/2/denominator)));
          // text->SetTextSize(0.035);
          // text->SetNDC();
          // text->Draw();
        // }

        g_deltaM_PDF[k][c]->Write();
        h_deltaM_PDF[k][c]->Write();
        c1->Write();

      }
    }
  }
    
    fout->Write();
    
    fout->Delete();
}

//.........这里部分代码省略.........
    pythiaInclMuonDatasetFiles.push_back("/home/sixie/hist/Compute" + fakeTypeString + "FakeRate/filler/006/Compute" + fakeTypeString + "FakeRate_s8-incmu_5_50-id9.root");
    pythiaInclMuonDatasetNames.push_back("s8-incmu_5_50-id9");


    //*****************************************************************************************
    //Construct the fake rate Histgrams. Fit with function.
    //*****************************************************************************************
    vector<TH1F*> fakeRateHistograms; fakeRateHistograms.clear();
    vector<TF1*> fakeRateFits; fakeRateFits.clear();
    vector<string> fakeRateLegendNames; fakeRateLegendNames.clear();

    //pythia version
    TH1F *FakeRatePt_pythia = NULL;
    if (faketype == 11)
      FakeRatePt_pythia = createFakeRatePtHist(pythiaQCDemDatasetFiles, pythiaQCDemDatasetNames,
                                               0,faketype, chargetype,eventType,
                                               (chargeTypeName + fakeTypeString + eventTypeName + "FakeRate_Pt_PythiaQCDEM").c_str());
    else if (faketype == 13) {
      //FakeRatePt_pythia = createFakeRatePtHist(pythiaInclMuonDatasetFiles, pythiaInclMuonDatasetNames, 
      //                                         0,faketype,chargetype,eventType,
      //                                         (chargeTypeName + fakeTypeString + eventTypeName + "FakeRate_Pt_PythiaInclMuon").c_str());
      FakeRatePt_pythia = createFakeRatePtHist(pythiaQCDemDatasetFiles, pythiaQCDemDatasetNames, 
                                               0,faketype, chargetype,eventType,
                                               (chargeTypeName + fakeTypeString + eventTypeName + "FakeRate_Pt_PythiaQCDEM").c_str());
    } else {
      cerr << "faketype out of range\n";
      assert(false);
    }
    FakeRatePt_pythia->GetYaxis()->SetTitle((chargeTypeName + fakeTypeString + eventTypeName + " Fake Rate").c_str());    
      
    TF1 *FakeRateFitFunction_pythia = NULL;
    if (fitFakeRate) {
      TF1 *FakeRateFitFunction_pythia = fitFakeRatePtHist(FakeRatePt_pythia);
      if (faketype == 11)
        FakeRateFitFunction_pythia->SetName((chargeTypeName + fakeTypeString + eventTypeName + "FakeRateFunction_Pt_PythiaQCDEM").c_str());
      else if (faketype == 13) {
      FakeRateFitFunction_pythia->SetName((chargeTypeName + fakeTypeString + eventTypeName + "FakeRateFunction_Pt_PythiaQCDEM").c_str());    
      }
    }

    fakeRateHistograms.push_back(FakeRatePt_pythia);
    if (fitFakeRate) fakeRateFits.push_back(FakeRateFitFunction_pythia);
    if (faketype == 11)
      fakeRateLegendNames.push_back(chargeTypeName + fakeTypeString + eventTypeName + " Pythia QCDEM");
    else if (faketype == 13)
      fakeRateLegendNames.push_back(chargeTypeName + fakeTypeString + eventTypeName + " Pythia QCDEM");
    else {
      cerr << "faketype out of range\n";
      assert(false);
    }

    //madgraph version
    TH1F *FakeRatePt_madgraph = createFakeRatePtHist(
      madgraphQCDDatasetFiles, madgraphQCDDatasetNames,0,faketype,chargetype,eventType,
      (chargeTypeName + fakeTypeString + "FakeRate_Pt_MadgraphQCD").c_str());
    FakeRatePt_madgraph->GetYaxis()->SetTitle((chargeTypeName +fakeTypeString + eventTypeName + " Fake Rate").c_str());    

    TF1 *FakeRateFitFunction_madgraph = NULL;
    if (fitFakeRate) {
      FakeRateFitFunction_madgraph = fitFakeRatePtHist(FakeRatePt_madgraph);
      FakeRateFitFunction_madgraph->SetName(
        (chargeTypeName + fakeTypeString + eventTypeName + "FakeRateFunction_Pt_MadgraphQCD").c_str());      
    }
    
    fakeRateHistograms.push_back(FakeRatePt_madgraph);
    if (fitFakeRate) fakeRateFits.push_back(FakeRateFitFunction_madgraph);
    fakeRateLegendNames.push_back(chargeTypeName + fakeTypeString + eventTypeName + " Madgraph QCD");
    
    //*****************************************************************************************
    //Draw the plots
    //*****************************************************************************************
    drawFakeRatePlots(fakeRateHistograms,fakeRateFits,fakeRateLegendNames,
                      (chargeTypeName + fakeTypeString + eventTypeName + "FakeRatePt").c_str(),maxY);

    //*****************************************************************************************
    //save them to the output file
    //*****************************************************************************************
    for (UInt_t i=0; i<fakeRateHistograms.size();i++) {
      //write the histograms. delete any existing one with the same name
      if (file->FindObjectAny(fakeRateHistograms[i]->GetName())) {
        string tmpname = fakeRateHistograms[i]->GetName();        
        file->Delete((tmpname+";*").c_str());
      }

      fakeRateHistograms[i]->SetDirectory(file);
      file->WriteTObject(fakeRateHistograms[i]);

      //write the fit functions. delete any existing one with the same name
      if (fitFakeRate){
        if (file->FindObjectAny(fakeRateFits[i]->GetName())) {
          string tmpname = fakeRateFits[i]->GetName();        
          file->Delete((tmpname+";*").c_str());
        }
        file->WriteTObject(fakeRateFits[i]);
      }
    }
  }
  file->Close();
  return;
}

//.........这里部分代码省略.........
		/*
		// make counter hists pretty (and distinguishable)
		tof1->SetTitle(data_cut);
		tof1->SetTitleSize(0.5);
		tof1->SetLineColor(kBlue);
		tof2->SetLineColor(kRed);
		trig->SetLineColor(kGreen);

		// draw histograms and save canvas
		char *cname = Form("%s/reco_%s_%d_%d_ccuts.png",
					  savepath,runtype,lens,(int)angle);
		TCanvas *ccanv = new TCanvas();
		gPad->SetLogy();
		tof1->Draw();
		tof2->Draw("same");
		trig->Draw("same");
		ccanv->Print(cname);
		*/

		t->Project("theta","theta",data_cut);
	}
	else
		t->Project("theta","theta",sim_cut);

	// preliminary paramters for fit of theta
	double height = theta->GetMaximum();
	double center = theta->GetBinCenter(theta->GetMaximumBin());
	double sigma  = 0.015;
	double slope  = 1.0;
	double shift  = height/3;

	// define custom fit, set parameters, and fit to data
	TF1 *gaus0 = new TF1("gaus0", "gaus+[3]*x+[4]", 
						 center-0.05, center+0.05);
	gaus0->SetParameters(height,center,sigma,slope,shift);
	gaus0->SetParNames("height","mean","sigma","slope","const");
	theta->Fit("gaus0","QR");
	double truesig = 1000*gaus0->GetParameter("sigma");
	
	// define strings for picture and root file names
	char *namt = Form("theta_{c} lens %d angle %d diff %1.2f",
					  lens,(int)angle,tdiff);
	char *pict = Form("%s/reco_%s_%d_%d_%1.2f.png",
					  savepath,runtype,lens,(int)angle,tdiff);

	/*
	// format histogram for saving
	TCanvas *c1 = new TCanvas();
	c1->cd();
	double max = 1.1*theta->GetMaximum();
	theta->SetTitle("");
	theta->GetXaxis()->SetLabelSize(0.06);
	theta->GetYaxis()->SetLabelSize(0.06);
	theta->GetXaxis()->SetRangeUser(0.61,1);
	theta->GetYaxis()->SetRangeUser(0,max);
	theta->Draw();
	c1->Print(pict);

	// save histogram again without fit parameters
   	gStyle->SetOptFit(0);
	c1->cd();
	theta->Draw();
	pict = Form("%s/reco_%s_%d_%d_%1.2f_nofit.png",
				savepath,runtype,lens,(int)angle,tdiff);
	c1->Print(pict);
	*/

	// print relevant info to screen and file
	char *outnm = Form("../macro/timecut_%s.tsv",runtype);
	ofstream out;
	out.open(outnm, fstream::in | fstream::out | fstream::app);
	out << angle 
		<< "\t" << lens  
		<< "\t" << tdiff
		<< "\t" << truesig
		<< "\t" << yield 
		<< "\t" << nph << endl;
	out.close();
	cout << "\033[1;31m" 
		 << "\nangle:\t" << angle
		 << "\nlens:\t"  << lens
		 << "\ndiff:\t"  << tdiff
		 << "\npeak:\t"  << diffmax
		 << "\nsigma:\t" << truesig
		 << "\nhits/e:\t" << nph
		 << "\nphts/e:\t" << yield 
		 << "\033[0m" <<  endl;
	TFile *tfile = new TFile(Form("timecut_%1.1f_%d_%1.2f.root",
								  angle,lens,tdiff),"create");
	TTree *ttree = new TTree("tcut","tcut");
	ttree->Branch("angle",&angle,"angle/D");
	ttree->Branch("lens",&lens,"lens/I");
	ttree->Branch("tcut",&tdiff,"tcut/D");
	ttree->Branch("sigma",&truesig,"truesig/D");
	ttree->Branch("nph",&yield,"yield/D");

	ttree->Fill();
	ttree->Write();
	tfile->Delete();
}