C++ TChain::GetBranch方法代码示例

void readmcsample(TH1D &nom, TH1D &rw, TH1D &nomratio, bool doGen){
 cout<<"reading data for "<<File_Signal_reco<<endl;
  TChain* t;
  if (!doGen) t= new TChain(reco_name.c_str(),reco_name.c_str());
  else t= new TChain(gen_name.c_str(),gen_name.c_str());
  int nfiles;
  if (!doGen) nfiles=t->Add(File_Signal_reco.c_str());
  else nfiles=t->Add(File_Signal_gen.c_str());

  TBranch *b_reco=t->GetBranch("reco");
  TBranch *b_truth=t->GetBranch("truth");
  TLeaf *l_mass=b_reco->GetLeaf("z_y");
  if (doMass) l_mass=b_reco->GetLeaf("z_m");
  if (doGen) {
    l_mass=b_truth->GetLeaf("z_y");
    if (doMass) l_mass=b_truth->GetLeaf("z_m");
  }
  TLeaf *l_phistar=b_reco->GetLeaf("z_phistar_dressed");
  TLeaf *l_phistar_true=b_truth->GetLeaf("z_phistar_dressed");
  int nweights;
  t->SetBranchAddress("weight_size",&nweights);
  t->GetEntry(0);
  cout<<"The sample has nweights: "<<nweights<<endl;
  double weights[nweights];
  int weightid[nweights];
  t->SetBranchAddress("weights",&weights);
  t->SetBranchAddress("weight_ids",&weightid);
  
  // TFile f_zy("ratio_y.root");
  // ratio = (TH1F*)f_zy.Get("ratio");
  // f_zy.Close();
  TFile f_zm("ratio_zmass.root");
  TH1F *ratio = (TH1F*)f_zm.Get("ratio_histo");
  cout<<"Entries: "<<t->GetEntries()<<endl;
  for (int i=0; i<t->GetEntries();i++){
    t->GetEntry(i);
    double phistar=l_phistar->GetValue();
    double mass=l_mass->GetValue();
    double phistar_true=l_phistar_true->GetValue();
    double weight =1;
    for (int w=0; w<nweights;w++){
      if (weightid[w]==1 || weightid[w]==2 || weightid[w]==12 || weightid[w]==13 || weightid[w]==20 || weightid[w]==30) {weight=weight*weights[w];}
    }
    int idx=ratio->GetXaxis()->FindBin(mass);
    if (ratio->GetBinContent(idx)!=0){
     if (!doGen){
	nom.Fill(phistar,weight);
	nomratio.Fill(phistar,weight);
	weight=weight*ratio->GetBinContent(idx);
	rw.Fill(phistar,weight);
     }
      else{
	nom.Fill(phistar_true,weight);
	nomratio.Fill(phistar_true,weight);
	weight=weight*ratio->GetBinContent(idx);
	rw.Fill(phistar_true,weight);
      }
    }
  }
}

TH1D* GetDataPhiStar() {
    cout << "reading data" << endl;
    vector<std::string> File_Data = Get_File_Data();
    TH1D* h_phistar = new TH1D("phistar", "phistar", nphistar, phistarBins);
    //TH1D *h_phistar= new TH1D("phistar","phistar",160,50,130);
    h_phistar->Sumw2();
    TChain* t = new TChain(reco_name.c_str(), reco_name.c_str());
    int nfiles;
    for (int i = 0; i < File_Data.size(); i++) {
        nfiles = t->Add(File_Data[i].c_str());
    }
    TBranch* b_reco = t->GetBranch("reco");
    TLeaf* l_m = b_reco->GetLeaf("z_m");
    TLeaf* l_phistar = b_reco->GetLeaf("z_phistar_dressed");
    TLeaf* l_e0_q = b_reco->GetLeaf("e_charge0");
    TLeaf* l_e1_q = b_reco->GetLeaf("e_charge1");
    cout << "Entries: " << t->GetEntries() << endl;
    for (int i = 0; i < t->GetEntries(); i++) {
        t->GetEntry(i);
        double E0_q = l_e0_q->GetValue();
        double E1_q = l_e1_q->GetValue();
        //if (E0_q!=E1_q) continue;
        h_phistar->Fill(l_phistar->GetValue(), 1);
        //h_phistar->Fill(l_m->GetValue(),1);
    }
    cout << "filled data phistar histogram" << endl;
    return h_phistar;
}

int Fitter::Save_sWeights()
{
	struct NameValTuple {string name; double val;};
	TFile* pFile = new TFile(outputSweights.c_str(), "RECREATE");
	TChain* pChain = new TChain;
	TTree* pTree = NULL;
	Long64_t nEntries = 0, iEntry=0;
	BranchProxy<double> bp;
	vector<NameValTuple> nameValTuples;
	
	if(!pFile){ cerr << "Error creating file" << endl; return 1; }
	if(InitChain(pChain)) return 1;	
	
	// Deactivate branches that would have the same name as the added branches
	for(auto& pVec_pdfs : {&sigPdfs, &bkgPdfs})
	for(auto& pdf : *pVec_pdfs)
	{
		NameValTuple nv = {pdf.GetYld().GetName()};
		string swBranchName = nv.name+"_sw";
		if(pChain->GetBranch(swBranchName.c_str()))
			pChain->SetBranchStatus(swBranchName.c_str(), 0);
		nameValTuples.push_back(nv);
	}
	
	pTree = pChain->CloneTree(0);
	bp.Connect(pChain, branchName.c_str());
	nEntries = pChain->GetEntries();
	
	for(auto& nv : nameValTuples)
	{
		string swBranchName = nv.name + "_sw";
		string branchTitle = swBranchName + "/D";
		
		pTree->Branch(swBranchName.c_str(), &nv.val, branchTitle.c_str());	
	}
	
	for(Long64_t i=0; i<nEntries; ++i)
	{
		pChain->GetEntry(i);
		if(!range.In(bp)) continue;
		
		for(auto& nv : nameValTuples) 
			nv.val = pSPlot->GetSWeight(iEntry, (nv.name+"_sw").c_str());
		
		++iEntry;
		pTree->Fill();
	}
	pTree->AutoSave();
	
	delete pFile;
	delete pChain;
	
	return 0;
}

TH1D* GetBGPhiStar(std::string FileName, double sampleweight) {
    TH1D* phistartemp = new TH1D("phistar", "phistar", nphistar, phistarBins);
    //TH1D *phistartemp= new TH1D("phistar","phistar",160,50,130);
    phistartemp->Sumw2();

    gErrorIgnoreLevel = kError;
    cout << "reading data for " << FileName << endl;
    TChain* t = new TChain(reco_name.c_str(), reco_name.c_str());
    int nfiles;
    nfiles = t->Add(FileName.c_str());

    TBranch* b_reco = t->GetBranch("reco");
    TLeaf* l_phistar = b_reco->GetLeaf("z_phistar_dressed");
    TLeaf* l_m = b_reco->GetLeaf("z_m");
    TLeaf* l_e0_q = b_reco->GetLeaf("e_charge0");
    TLeaf* l_e1_q = b_reco->GetLeaf("e_charge1");
    int nweights;
    t->SetBranchAddress("weight_size", &nweights);
    t->GetEntry(0);
    cout << "The sample has nweights: " << nweights << endl;
    double weights[nweights];
    int weightid[nweights];
    t->SetBranchAddress("weights", &weights);
    t->SetBranchAddress("weight_ids", &weightid);

    cout << "reading data for " << FileName << endl;
    cout << "Entries: " << t->GetEntries() << endl;
    for (int i = 0; i < t->GetEntries(); i++) {
        t->GetEntry(i);
        double E0_q = l_e0_q->GetValue();
        double E1_q = l_e1_q->GetValue();
        //cout<<E0_q<<" "<<E1_q<<endl;
        //if (E0_q!=E1_q) continue;
        double phistar = l_phistar->GetValue();
        double weight = sampleweight;
        for (int w = 0; w < nweights; w++) {
            if (weightid[w] == 1 || weightid[w] == 2 || weightid[w] == 12
                    || weightid[w] == 13 || weightid[w] == 20 || weightid[w] == 30) {
                weight = weight * weights[w];
            }
        }
        phistartemp->Fill(phistar, weight);
        //phistartemp->Fill(l_m->GetValue(),weight);
    }
    cout << "done reading data for " << FileName << endl;
    return phistartemp;
}

void GetGen(double sampleweight, TH1D* &h_phistar, vector<TH1D*> &h_cteq, vector<TH1D*> &h_fsr_pileup) {
  gErrorIgnoreLevel = kError;
  cout << "reading signal gen" << endl;
  h_phistar = new TH1D("phistar", "phistar", nbins, 0, nbins);
  h_phistar->Sumw2();
  TChain* t = new TChain(gen_name.c_str(), gen_name.c_str());
  if (doMG) {
    if (elec == 0) t->Add(File_Signal_gen.c_str());
    else if (elec == 1) t->Add(File_Signal_gen_born.c_str());
    else t->Add(File_Signal_gen_bare.c_str());
  } else {
    if (elec == 0) t->Add(File_Powheg_gen.c_str());
    else if (elec == 1) t->Add(File_Powheg_gen_born.c_str());
    else t->Add(File_Powheg_gen_bare.c_str());
  }
  ///t->SetBranchStatus("event_info", 0); //to disable all branches
  //  t->SetBranchStatus("reco", 0); //to disable all branches
  TBranch *b_truth = t->GetBranch("truth");
  TBranch *b_reco = t->GetBranch("reco");
  TLeaf * l_ZY = b_truth->GetLeaf("z_y");
  if (elec == 1) l_ZY = b_reco->GetLeaf("z_yBorn"); //TO BE CHANGED
  if (elec == 2) l_ZY = b_reco->GetLeaf("z_yNaked");
  TLeaf *l_phistar = b_truth->GetLeaf("z_phistar_dressed");
  if (elec == 1) l_phistar = b_truth->GetLeaf("z_phistar_born");
  if (elec == 2) l_phistar = b_truth->GetLeaf("z_phistar_naked");
  int nweights;
  int nwcteq;
  t->SetBranchAddress("weight_size", &nweights);
  t->SetBranchAddress("weight_cteq_size", &nwcteq);
  t->GetEntry(0);
  double weights[nweights];
  int weightid[nweights];
  double weights_cteq[nwcteq];
  t->SetBranchAddress("weights", &weights);
  t->SetBranchAddress("weight_ids", &weightid);
  t->SetBranchAddress("weights_cteq", &weights_cteq);
  double weight_fsr;
  t->SetBranchAddress("weight_fsr", &weight_fsr);
  cout << "Entries: " << t->GetEntries() << endl;
  h_cteq = EmptyHVec(nwcteq);
  h_fsr_pileup = EmptyHVec(3);

  TH1D* h_PDFWeight = new TH1D("h_PDFWeight", "h_PDFWeight", 10000, 0, 100);
  h_PDFWeight->Sumw2();
  // cout << "Nweightcteq=" << nwcteq << endl;

  for (int i = 0; i < t->GetEntries(); i++) {
    // for (int i=0; i<50000;i++){
    t->GetEntry(i);
    double weight = sampleweight;
    double pdfnorm = weights_cteq[0];
    double weightpu_0 = 0;
    double weightpu_p = 0;
    double weightpu_m = 0;
    double phistar = l_phistar->GetValue();
    double Z_Y = l_ZY->GetValue();
    int bin=GetBin(phistar,Z_Y);
    //cout<<pdfnorm<<" "<<weights_cteq[0]<<endl;
    for (int w = 0; w < nweights; w++) {
      if (weightid[w] == 1 || weightid[w] == 2) {
        weight = weight * weights[w];
      }
      //if (weightid[w]==1) {weight=weight*weights[w];}
      if (weightid[w] == 2) weightpu_0 = weights[w];
      if (weightid[w] == 3) weightpu_p = weights[w];
      if (weightid[w] == 4) weightpu_m = weights[w];
    }
    // if (weightpu_0 == 0 || weightpu_p == 0 || weightpu_m == 0) cout << "pile-up weights not there" << endl;
    h_phistar->Fill(bin, weight);
    if (pdfnorm != 0) {
      for (int w = 0; w < nwcteq; w++) {
        h_cteq[w] ->Fill(bin, weight * weights_cteq[w] / pdfnorm);
        if (w != 0) {
          h_PDFWeight->Fill(weights_cteq[w] / pdfnorm, 1);
          // cout<<weights_cteq[w] / pdfnorm<<endl;
        }
      }
    }
    if (doMG) weight_fsr = 1;
    h_fsr_pileup[0]->Fill(bin, weight * weight_fsr);
    if (weightpu_0 != 0) {
      h_fsr_pileup[1]->Fill(bin, weight * weightpu_p / weightpu_0);
      h_fsr_pileup[2]->Fill(bin, weight * weightpu_m / weightpu_0);
    }
  }
  cout << "done reading signal gen" << endl;
  delete t;
  return;
}

//.........这里部分代码省略.........
    std::cout << " -- output file " << outputFile->GetName() << " successfully opened." << std::endl;
  }
  outputFile->cd();


  ///initialise PU density object

  HGCSSPUenergy puDensity("data/EnergyDensity.dat");


    //////////////////////////////////////////////////
    //////////////////////////////////////////////////
    ///////// positionFit /////////////////////////////
    //////////////////////////////////////////////////
    //////////////////////////////////////////////////
  
  const unsigned nEvts = ((pNevts > lSimTree->GetEntries() || pNevts==0) ? static_cast<unsigned>(lSimTree->GetEntries()) : pNevts) ;
  

  PositionFit lChi2Fit(nSR,residualMax,nLayers,nSiLayers,applyPuMixFix,debug);
  lChi2Fit.initialise(outputFile,"PositionFit",outFolder,geomConv,puDensity);

  //try getting z position from input file, if doesn't exit,
  //perform first loop over simhits to find z positions of layers
  if ((redoStep<2 && !lChi2Fit.getZpositions(versionNumber)) || redoStep>1)
    lChi2Fit.getZpositions(versionNumber,lSimTree,nEvts);
  
  //perform second loop over events to find positions to fit and get energies
  SignalRegion SignalEnergy(outFolder, nLayers, nEvts, geomConv, puDensity,applyPuMixFix,versionNumber);
  SignalEnergy.initialise(outputFile,"Energies");

  //initialise
  bool dofit = redoStep>0 || !SignalEnergy.initialiseFitPositions();
  if (!dofit && redoStep==0) std::cout << " -- Info: fit positions taken from file on disk." << std::endl;
  else std::cout << " -- Info: redoing least square fit." << std::endl;

  if (redoStep>0 || (redoStep==0 && dofit)){
    lChi2Fit.getInitialPositions(lSimTree,lRecTree,nEvts);
    lChi2Fit.finaliseErrorMatrix();
    lChi2Fit.initialiseLeastSquareFit();
  }
  
  //loop on events
  HGCSSEvent * event = 0;
  std::vector<HGCSSSamplingSection> * ssvec = 0;
  std::vector<HGCSSSimHit> * simhitvec = 0;
  std::vector<HGCSSRecoHit> * rechitvec = 0;
  std::vector<HGCSSGenParticle> * genvec = 0;
  unsigned nPuVtx = 0;

  lSimTree->SetBranchAddress("HGCSSEvent",&event);
  lSimTree->SetBranchAddress("HGCSSSamplingSectionVec",&ssvec);
  lSimTree->SetBranchAddress("HGCSSSimHitVec",&simhitvec);
  lSimTree->SetBranchAddress("HGCSSGenParticleVec",&genvec);

  lRecTree->SetBranchAddress("HGCSSRecoHitVec",&rechitvec);
  if (lRecTree->GetBranch("nPuVtx")) lRecTree->SetBranchAddress("nPuVtx",&nPuVtx);
  const unsigned nRemove = 12;
  std::vector<unsigned> lToRemove;
  unsigned list[nRemove] = {25,27,15,1,10,3,18,5,12,7,23,20};

  for (unsigned ievt(0); ievt<nEvts; ++ievt){//loop on entries
    if (debug) std::cout << "... Processing entry: " << ievt << std::endl;
    else if (ievt%50 == 0) std::cout << "... Processing entry: " << ievt << std::endl;

    lSimTree->GetEntry(ievt);
    lRecTree->GetEntry(ievt);
    if (dofit) {
      bool found = lChi2Fit.setTruthInfo(genvec,1);
      if (!found) continue;
      //mask layers in turn
      lToRemove.clear();
      for (unsigned r(0); r<nRemove+1;++r){
	FitResult fit;
	if ( lChi2Fit.performLeastSquareFit(ievt,fit,lToRemove)==0 ){
	  //SignalEnergy.fillEnergies(ievt,(*ssvec),(*simhitvec),(*rechitvec),nPuVtx,fit);
	}
	else std::cout << " -- remove " << r << " Fit failed." << std::endl;
	if (r<nRemove) lToRemove.push_back(list[r]);
      }

    }
    else SignalEnergy.fillEnergies(ievt,(*ssvec),(*simhitvec),(*rechitvec),nPuVtx);

  }//loop on entries

  //finalise

  if (dofit) lChi2Fit.finaliseFit();
  SignalEnergy.finalise();

  outputFile->Write();
  //outputFile->Close();
  
  std::cout << " - End of egammaResolution program." << std::endl;

  return 0;
  

}//main

int main(int argc, char *argv[])
{
	TString runNumber, runType, filename;
	Int_t barrelCode, retCode;
	Int_t iModule;
	TCanvas *chi, *clo;
	struct stat statBuf;

	//cout << getenv("PATH") << endl;
	//cout << getenv("ROOTSYS") << endl;

	if (argc < 3)
	{
		//cout << "Usage:"<< endl << "    ./checkPulseHeight <runNumber> <runType> [<filename>]" << endl;
		results.addError("Usage:\n    ./checkPulseHeight <runNumber> <runType> [<filename>]");
		return 0;
	}

	runNumber = argv[1];
	runType = argv[2];

	if (argc > 3)
		filename = argv[3];
	else
	{
		filename  = "tiletb_";
		filename += runNumber;
		filename += "_";
		filename += runType;
		filename += ".0.root";
		TString dir = "/work/tiletransfer/wis/";
		TString aux = dir + filename;
		retCode = stat((char*)aux.Data(), &statBuf);
		if (retCode < 0)
		{
			dir = "/castor/cern.ch/user/t/tilebeam/commissioning/";
			aux = dir + filename;
			retCode = rfio_stat((char*)aux.Data(), &statBuf);
			if (retCode < 0)
			{
				//cout << "<h3>Ntuple does not exist in CASTOR default directory.</h3>" << endl;
				results.addError("Ntuple does not exist in CASTOR default directory.");
				return 0;
			}
			else
				filename = "rfio:" + aux;
		}
		else
		{
			filename = aux;
		}
	}
	//cout<<filename<<endl;
	TChain *t = new TChain("TileRec/h1000");
	t->Add(filename);

	Int_t nevt = t->GetEntries();

	results.setMacroName("Pulse Height");
	results.setParameter("Run Number", runNumber);
	results.setParameter("Number of Events", nevt);

/**********************************************************/

	Int_t i, j, k;

    /******************* Set names ********************/
	TString moduleName[Nmodules];
	TString barrelName[Nbarrels];

	// barrel names
	barrelName[0] = "A";
	barrelName[1] = "C";
	barrelName[2] = "EBA";
	barrelName[3] = "EBC";

	// module names
	for (i = 0; i < 2; i++)
	{
		for (j = 0; j < 64; j++)
		{
			moduleName[i*64 + j] = barrelName[i];
			moduleName[i*64 + j] += (j + 1);
		}
	}
	/**************************************************/

	Int_t SampleHi[Nmodules][48][7];
	Int_t SampleLo[Nmodules][48][7];
	TBranch *b_SampleHi[Nmodules];
	TBranch *b_SampleLo[Nmodules];
	bool moduleInNtuple[Nmodules];

	TString branchName;
	for (iModule = 0; iModule < Nmodules; iModule++)
	{
		branchName = "Sample" + moduleName[iModule] + "hi";
		t->SetBranchAddress(branchName, SampleHi[iModule]);
		b_SampleHi[iModule] = t->GetBranch(branchName.Data());

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

void GetBinM(vector<RooUnfoldResponse*> &BinM, vector<TH1D*> &h_gen, vector<TH1D*> &h_reco, bool madgraph=1, int elec=0){

  TH2D* BinMigration=new TH2D("BinMigration","BinMigration",nphistar,phistarBins,nphistar,phistarBins);
  BinMigration->Sumw2();
  TH1D* Gen=new TH1D("Gen","gen",nphistar,phistarBins);
  Gen->Sumw2();
  TH1D* Reco=new TH1D("Reco","Reco",nphistar,phistarBins);
  Reco->Sumw2();

  TChain* t = new TChain(reco_name.c_str(),reco_name.c_str());
  int nfiles;
  if (madgraph)  nfiles=t->Add(File_Signal_reco.c_str());
  else nfiles=t->Add(File_Powheg_reco.c_str());

  TBranch *b_reco=t->GetBranch("reco");
  TBranch *b_truth=t->GetBranch("truth");
  TBranch *b_event=t->GetBranch("event_info");
  TLeaf *l_phistar=b_reco->GetLeaf("z_phistar_dressed");
  TLeaf *l_phistar_true=b_truth->GetLeaf("z_phistar_dressed");
  TLeaf *l_en=b_event->GetLeaf("event_number");
  if (elec==1) l_phistar_true=b_truth->GetLeaf("z_phistar_born");
  if (elec==2) l_phistar_true=b_truth->GetLeaf("z_phistar_naked");
  int nweights;
  t->SetBranchAddress("weight_size",&nweights);
  t->GetEntry(0);
  cout<<"The sample has nweights: "<<nweights<<endl;
  double weights[nweights];
  int weightid[nweights];
  t->SetBranchAddress("weights",&weights);
  t->SetBranchAddress("weight_ids",&weightid);

  cout<<"Entries: "<<t->GetEntries()<<endl;
  cout<<"reading signal reco "<<endl;
  for (int i=0; i<t->GetEntries();i++){
    // if (!madgraph && (i>N_MC+50000 || i<50000) &&N_MC!=-1) continue;
    if (!madgraph && i>N_MC &&N_MC!=-1) continue;
    //for (int i=0; i<50000;i++){
    t->GetEntry(i);
    double phistar=l_phistar->GetValue();
    double phistar_true=l_phistar_true->GetValue();
    double weight =1;
    double eventn=l_en->GetValue();
    // else en=eventn;
    for (int w=0; w<nweights;w++){
      if (weightid[w]==1 || weightid[w]==2 || weightid[w]==12 || weightid[w]==13 || weightid[w]==20 || weightid[w]==30) {weight=weight*weights[w];}
    }
    Gen ->Fill(phistar_true,weight);
    Reco->Fill(phistar,weight);
    BinMigration->Fill(phistar,phistar_true,weight); 
  }
  if (!madgraph){
    TCanvas* BM = new TCanvas("BM","BM",800,900);
    BinMigration->GetXaxis()->SetRangeUser(0.001,3.2);
    BinMigration->GetXaxis()->SetTitle("#phi^{*}(reconstructed)");
    BinMigration->GetXaxis()->SetTitleOffset(0.8);
    BinMigration->GetXaxis()->SetTitleSize(0.04);
    BinMigration->GetXaxis()->SetLabelOffset(-0.01);
    BinMigration->GetXaxis()->SetLabelSize(0.04);
    BinMigration->GetYaxis()->SetTitleOffset(1.05);
    BinMigration->GetYaxis()->SetTitleSize(0.04);
    BinMigration->GetYaxis()->SetLabelSize(0.04);
    BinMigration->GetYaxis()->SetRangeUser(0.001,3.2);
    BinMigration->GetYaxis()->SetTitle("#phi^{*}(generated)");
    BinMigration->GetZaxis()->SetTitleOffset(-0.004);
    BinMigration->SetStats(0);
    BinMigration->SetBit( TH2::kNoTitle, true );
    BM->SetLogx();
    BM->SetLogy();
    BM->SetLogz();
    BinMigration->Draw("COLZ");
    TLatex mark;
    mark.SetTextSize(0.04);
    mark.SetTextFont(42);
    mark.SetNDC(true);
    mark.DrawLatex(0.19,0.80,"Powheg");
  }
  RooUnfoldResponse* BinM1  =new RooUnfoldResponse (Reco,Gen,BinMigration);
  BinM.push_back(BinM1);
  h_gen.push_back(Gen);
  h_reco.push_back(Reco);
  for (int i=0; i<N_TOY; i++){
    TH1D* Recotemp=new TH1D("Reco","Reco",nphistar,phistarBins);
    Recotemp->Sumw2();
    TH2D* BinMigrationtemp=new TH2D("BinMigration","BinMigration",nphistar,phistarBins,nphistar,phistarBins);
    BinMigrationtemp->Sumw2();
    for (int j=0; j<nphistar; j++){
      double x=gRandom->Gaus(Reco->GetBinContent(j+1),Reco->GetBinError(j+1));
      Recotemp->SetBinContent(j+1,x);
      Recotemp->SetBinError(j+1,Reco->GetBinError(j+1));
      for (int k=0; k<nphistar; k++){
	double mean=BinMigration->GetBinContent(j+1,k+1);
	if (mean==0) continue;
	double error=BinMigration->GetBinError(j+1,k+1);
	if (mean/error<5){
	  x=gRandom->Poisson(mean);
	}
	else{
	  x=gRandom->Gaus(mean,error);
	}
	BinMigrationtemp->SetBinContent(j+1,k+1,x);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
        }
    }
/*
    HiForestInfoController hfic(treeHiForestInfo);
    std::cout<<"### HiForestInfo Tree ###"<< std::endl;
    hfic.printHiForestInfo();
    std::cout<<"###"<< std::endl;
*/
    treeHLT->SetBranchStatus("*",0);     // disable all branches
    treeHLT->SetBranchStatus("HLT_HI*SinglePhoton*Eta*",1);     // enable photon branches
    treeHLT->SetBranchStatus("HLT_HI*DoublePhoton*Eta*",1);     // enable photon branches

    float vz;
    Int_t hiBin;
    UInt_t run, lumis;
    ULong64_t event;
    float pthat, pthatWeight;
    treeHiEvt->SetBranchAddress("vz",&vz);
    treeHiEvt->SetBranchAddress("hiBin",&hiBin);
    treeHiEvt->SetBranchAddress("run", &run);
    treeHiEvt->SetBranchAddress("evt", &event);
    treeHiEvt->SetBranchAddress("lumi", &lumis);
    if(isMC) {
        //treeHiEvt->Branch("pthatWeight", &pthatWeight, "pthatWeight/F");
        treeHiEvt->SetBranchAddress("pthat", &pthat);
    }   

    // specify explicitly which branches to store, do not use wildcard
    treeSkim->SetBranchStatus("*",0);

    Int_t pcollisionEventSelection;  // this filter is used for HI.
    if (isHI) {
        treeSkim->SetBranchStatus("pcollisionEventSelection",1);
        if (treeSkim->GetBranch("pcollisionEventSelection")) {
            treeSkim->SetBranchAddress("pcollisionEventSelection",&pcollisionEventSelection);
        }
        else {   // overwrite to default
            pcollisionEventSelection = 1;
            std::cout<<"could not get branch : pcollisionEventSelection"<<std::endl;
            std::cout<<"set to default value : pcollisionEventSelection = "<<pcollisionEventSelection<<std::endl;
        }
    }
    else {
        pcollisionEventSelection = 0;    // default value if the collision is not HI, will not be used anyway.
    }
    Int_t pPAprimaryVertexFilter;    // this filter is used for PP.
    if (isPP) {
        treeSkim->SetBranchStatus("pPAprimaryVertexFilter",1);
        if (treeSkim->GetBranch("pPAprimaryVertexFilter")) {
            treeSkim->SetBranchAddress("pPAprimaryVertexFilter",&pPAprimaryVertexFilter);
        }
        else {   // overwrite to default
            pPAprimaryVertexFilter = 1;
            std::cout<<"could not get branch : pPAprimaryVertexFilter"<<std::endl;
            std::cout<<"set to default value : pPAprimaryVertexFilter = "<<pPAprimaryVertexFilter<<std::endl;
        }
    }
    else {
        pPAprimaryVertexFilter = 0;      // default value if the collision is not PP, will not be used anyway.
    }
    Int_t pBeamScrapingFilter;   // this filter is used for PP.
    if (isPP) {
        treeSkim->SetBranchStatus("pBeamScrapingFilter",1);
        if (treeSkim->GetBranch("pBeamScrapingFilter")) {
            treeSkim->SetBranchAddress("pBeamScrapingFilter",&pBeamScrapingFilter);
        }

void readTree(TString filelist="/global/project/projectdirs/star/pwg/starhf/zamiller/minBiasTemplates101015/Trees.list",TString outFile="/global/project/projectdirs/star/pwg/starhf/zamiller/minBiasTemplates101015/readTreeOutSubmitTEST.root")
{

  //TString filelist=Form("%s/%s",localDir,fileIn);
  //TString outFile=Form("%s/%s",outDir,fileOut);

  TFile fout(outFile,"RECREATE");
  TH1F::SetDefaultSumw2(); 
  TH1D::SetDefaultSumw2(); 

  for(Int_t c=0; c< numPtBins; c++)
  {
    lowpt[c] = anaConst::lpt[c];
    highpt[c] = anaConst::hpt[c];
  }
  for(Int_t ptbin = 0; ptbin < numPtBins; ptbin++)
  {
    hdEtaRawce[ptbin] = new TH1D(Form("hdEtaRawce_%i",ptbin),"hdEtaRawce",100, -5,5);
    hdEtaRawbe[ptbin] = new TH1D(Form("hdEtaRawbe_%i",ptbin),"hdEtaRawbe",100, -5,5);
    hEventTallyce[ptbin] = new TH1D(Form("ceEventTally_%i",ptbin),"ceEvent Tally",1,0,1);
    hEventTallybe[ptbin] = new TH1D(Form("beEventTally_%i",ptbin),"beEvent Tally",1,0,1);
    hEventTallybce[ptbin] = new TH1D(Form("bceEventTally_%i",ptbin),"bceEvent Tally",1,0,1);
    hdPhiRawce[ptbin] = new TH1D(Form("hdPhiRawce_%i",ptbin),"hdPhiRawce",Phibin, -10,10);
    hdPhiRawbe[ptbin] = new TH1D(Form("hdPhiRawbe_%i",ptbin),"hdPhiRawbe",Phibin, -10,10);
    hdPhiRawbce[ptbin] = new TH1D(Form("hdPhiRawbce_%i",ptbin),"hdPhiRawbce",Phibin, -10,10);
    hept[ptbin] = new TH1D(Form("hept_%i",ptbin),"hept",200,0.,20.);
  }

  // Make Chain
  TChain* chain = new TChain("tree");
  int nfile = 0;
  /*  nfile += chain->Add("Oct30_set13/pythia_tree_Oct30_13.root");
      nfile += chain->Add("Oct30_set14/pythia_tree_Oct30_14.root");
      nfile += chain->Add("Oct30_set15/pythia_tree_Oct30_15.root");
      nfile += chain->Add("Oct30_set16/pythia_tree_Oct30_16.root");
      nfile += chain->Add("Oct30_set17/pythia_tree_Oct30_17.root");
  //  nfile += chain->Add("Oct30_set12/pythia_tree_Oct30_12.root");
  //  nfile += chain->Add("Oct26_set3/pythia_tree_Oct26_3.root");
  // nfile += chain->Add("Oct18_set4/pythia_tree_Oct18_4.root");
  //nfile += chain->Add("liweiTemplate_part2.root");*/
  char filename[1000];
  ifstream fstream(filelist.Data());
  int ifile = 0;
  while (fstream >> filename)
  {
    nfile+= chain->Add(filename);
    std::cout << "Got File: " << filename << std::endl;
  }

  std::cout <<"Added "<<nfile<<" files"<<std::endl;
  std::cout<<"# entries in chain: "<<chain->GetEntries()<<std::endl;
  if (chain == 0) return;

  int ceNtrigger=0;
  int beNtrigger=0;
  int bceNtrigger=0;
  int ptbin,maxptbin;

  //define variables
  Int_t   Event, numberofcElectrons, numberofbElectrons,numberofbcElectrons, numberofHadrons, noTracks;   //
  Float_t celectron_id,celectron_status,celectron_pt,celectron_pz,celectron_phi,celectron_eta,celectron_y;                //track info
  Float_t belectron_id,belectron_status,belectron_pt,belectron_pz,belectron_phi,belectron_eta,belectron_y;
  Float_t bcelectron_id,bcelectron_status,bcelectron_pt,bcelectron_pz,bcelectron_phi,bcelectron_eta,bcelectron_y;
  Float_t assoh_id,assoh_status,assoh_pt,assoh_pz,assoh_phi,assoh_eta,assoh_y;

  int goodEvent = 0;


  Long64_t nentries = chain->GetEntriesFast();
  int x = 0; int n = nentries; int w = 25;

  for(int i = 0; i < nentries; i++){
    Long64_t ientry = chain->LoadTree(i);
    if (ientry < 0) break;
    TBranch *b_destep = chain->GetBranch("hf2eDecay");

    TLeaf* leaf_Event_id            = b_destep->GetLeaf("Event_id");
    TLeaf* leaf_refMult             = b_destep->GetLeaf("refMult");
    TLeaf* leaf_numberofcElectrons  = b_destep->GetLeaf("numberofcElectrons");
    TLeaf* leaf_numberofbElectrons  = b_destep->GetLeaf("numberofbElectrons");
    TLeaf* leaf_numberofbcElectrons = b_destep->GetLeaf("numberofbcElectrons");
    TLeaf* leaf_numberofHadrons     = b_destep->GetLeaf("numberofHadrons");
    TLeaf* leaf_noTracks            = b_destep->GetLeaf("noTracks");

    TLeaf* leaf_ce_id               = b_destep->GetLeaf("ce_id");
    TLeaf* leaf_ce_status           = b_destep->GetLeaf("ce_status");
    TLeaf* leaf_ce_pt               = b_destep->GetLeaf("ce_pt");
    TLeaf* leaf_ce_pz               = b_destep->GetLeaf("ce_pz");
    TLeaf* leaf_ce_phi              = b_destep->GetLeaf("ce_phi");
    TLeaf* leaf_ce_eta              = b_destep->GetLeaf("ce_eta");
    TLeaf* leaf_ce_y                = b_destep->GetLeaf("ce_y");

    TLeaf* leaf_be_id               = b_destep->GetLeaf("be_id");
    TLeaf* leaf_be_status           = b_destep->GetLeaf("be_status");
    TLeaf* leaf_be_pt               = b_destep->GetLeaf("be_pt");
    TLeaf* leaf_be_pz               = b_destep->GetLeaf("be_pz");
    TLeaf* leaf_be_phi              = b_destep->GetLeaf("be_phi");
    TLeaf* leaf_be_eta              = b_destep->GetLeaf("be_eta");
    TLeaf* leaf_be_y                = b_destep->GetLeaf("be_y");

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

void GetBinM(RooUnfoldResponse* &BinM, TH1D* &h_gen, TH1D* &h_reco, RooUnfoldResponse* &BinM_1, TH1D* &h_gen_1, TH1D* &h_reco_1, RooUnfoldResponse* &BinM_2, TH1D* &h_gen_2, TH1D* &h_reco_2, bool madgraph=1, int elec=0){

  TH2D* BinMigration=new TH2D("BinMigration","BinMigration",nbins,0,nbins,nbins,0,nbins);
  BinMigration->Sumw2();
  TH1D* Dphistar=new TH1D("Dphistar","Dphistar",2000,-1,1);
  Dphistar->Sumw2();
  TH1D* Dy=new TH1D("Dy","Dy",2000,-1,1);
  Dy->Sumw2();
  TChain* t = new TChain(reco_name.c_str(),reco_name.c_str());
  int nfiles;
  if (madgraph){
    if (elec==0) nfiles=t->Add(File_Signal_reco_d.c_str());
    if (elec==1) nfiles=t->Add(File_Signal_reco_b.c_str());
    if (elec==2) nfiles=t->Add(File_Signal_reco_n.c_str());
  }
  else {
    if (elec==0) nfiles=t->Add(File_Powheg_reco_d.c_str());
    if (elec==1) nfiles=t->Add(File_Powheg_reco_b.c_str());
    if (elec==2) nfiles=t->Add(File_Powheg_reco_n.c_str());
  }

  TBranch *b_reco=t->GetBranch("reco");
  TBranch *b_truth=t->GetBranch("truth");
  TLeaf *l_phistar=b_reco->GetLeaf("z_phistar_dressed");
  TLeaf *l_phistar_true=b_truth->GetLeaf("z_phistar_dressed");
  if (elec==1) l_phistar_true=b_truth->GetLeaf("z_phistar_born");
  if (elec==2) l_phistar_true=b_truth->GetLeaf("z_phistar_naked");
  TLeaf *l_y=b_reco->GetLeaf("z_y");
  TLeaf *l_y_true=b_truth->GetLeaf("z_y");
  if (elec==1) l_y_true=b_reco->GetLeaf("z_yBorn");
  if (elec==2) l_y_true=b_reco->GetLeaf("z_yNaked");
  int nweights;
  t->SetBranchAddress("weight_size",&nweights);
  t->GetEntry(0);
  cout<<"The sample has nweights: "<<nweights<<endl;
  double weights[nweights];
  int weightid[nweights];
  t->SetBranchAddress("weights",&weights);
  t->SetBranchAddress("weight_ids",&weightid);

  cout<<"Entries: "<<t->GetEntries()<<endl;
  h_gen   = new TH1D("phistar","phistar",nbins,0,nbins);
  h_gen->Sumw2();
  h_reco  = new TH1D("phistar","phistar",nbins,0,nbins);
  h_reco->Sumw2();
  h_gen_1 = new TH1D("phistar","phistar",nbins,0,nbins);
  h_gen_1->Sumw2();
  h_reco_1= new TH1D("phistar","phistar",nbins,0,nbins);
  h_reco_1->Sumw2();
  h_gen_2 = new TH1D("phistar","phistar",nbins,0,nbins);
  h_gen_2->Sumw2();
  h_reco_2= new TH1D("phistar","phistar",nbins,0,nbins);
  h_reco_2->Sumw2();
  BinM  =new RooUnfoldResponse (h_reco,h_gen);
  BinM_1=new RooUnfoldResponse (h_reco,h_gen);
  BinM_2=new RooUnfoldResponse (h_reco,h_gen);

  cout<<"reading signal reco "<<endl;
  for (int i=0; i<t->GetEntries();i++){
    // if (i>20000) continue;
    //for (int i=0; i<50000;i++){
    t->GetEntry(i);
    double weight =1;
    for (int w=0; w<nweights;w++){
      if (weightid[w]==1 || weightid[w]==2 || weightid[w]==12 || weightid[w]==13 || weightid[w]==20 || weightid[w]==30) {weight=weight*weights[w];}
    }
    double phistar=l_phistar->GetValue();
    double phistar_true=l_phistar_true->GetValue();
    Dphistar->Fill((phistar_true-phistar)/phistar_true,weight);
    double y=l_y->GetValue();
    double y_true=l_y_true->GetValue();
    Dy->Fill((y_true-y)/y_true,weight);

    int bin=GetBin(phistar,y);
    int bin_true=GetBin(phistar_true,y_true);

    h_gen ->Fill(bin_true,weight);
    h_reco->Fill(bin,weight);
    if (N_MC==-1 || i<N_MC) BinM  ->Fill(bin,bin_true,weight); 
    BinMigration->Fill(bin,bin_true,weight);
    if (i%2==0){
      h_gen_1 ->Fill(bin_true,weight);
      h_reco_1->Fill(bin,weight);
      if (N_MC==-1 || i<N_MC) BinM_1  ->Fill(bin,bin_true,weight); 
    }
    else{
      h_gen_2 ->Fill(bin_true,weight);
      h_reco_2->Fill(bin,weight);
      if (N_MC==-1 || i<N_MC) BinM_2  ->Fill(bin,bin_true,weight); 
    }
  }
////////////////////////////here is where I left.
  TLatex mark;
  mark.SetTextSize(0.04);
  mark.SetTextFont(42);
  mark.SetNDC(true);

//  BinMigration->GetXaxis()->SetRangeUser(0.001,3.2);
  BinMigration->GetXaxis()->SetTitle("(#phi*,y) bin (reconstructed)");
  BinMigration->GetXaxis()->SetTitleOffset(0.8);
//.........这里部分代码省略.........

//.........这里部分代码省略.........
       TChain* treeSkim  = new TChain("skimanalysis/HltTree");
       TChain* treeHiForestInfo = new TChain("HiForest/HiForestInfo");

       for (std::vector<std::string>::iterator it = inputFiles.begin() ; it != inputFiles.end(); ++it) {
          treeHLT->Add((*it).c_str());
          treeggHiNtuplizer->Add((*it).c_str());
          treeHiEvt->Add((*it).c_str());
          treeSkim->Add((*it).c_str());
          treeHiForestInfo->Add((*it).c_str());
       }

       HiForestInfoController hfic(treeHiForestInfo);
       std::cout<<"### HiForestInfo Tree ###"<< std::endl;
       hfic.printHiForestInfo();
       std::cout<<"###"<< std::endl;

       treeHLT->SetBranchStatus("*",0);     // disable all branches
       treeHLT->SetBranchStatus("HLT_HI*SinglePhoton*Eta*",1);     // enable photon branches
       treeHLT->SetBranchStatus("HLT_HI*DoublePhoton*Eta*",1);     // enable photon branches
       treeHLT->SetBranchStatus("*DoubleMu*",1);                      // enable muon branches
       treeHLT->SetBranchStatus("HLT_HIL1Mu*",1);                     // enable muon branches
       treeHLT->SetBranchStatus("HLT_HIL2Mu*",1);                     // enable muon branches
       treeHLT->SetBranchStatus("HLT_HIL3Mu*",1);                     // enable muon branches
       
       // specify explicitly which branches to store, do not use wildcard
       treeHiEvt->SetBranchStatus("*",1);

       // specify explicitly which branches to store, do not use wildcard
       treeSkim->SetBranchStatus("*",0);

       Int_t pcollisionEventSelection;  // this filter is used for HI.
       if (isHI) {
           treeSkim->SetBranchStatus("pcollisionEventSelection",1);
           if (treeSkim->GetBranch("pcollisionEventSelection")) {
               treeSkim->SetBranchAddress("pcollisionEventSelection",&pcollisionEventSelection);
           }
           else {   // overwrite to default
               pcollisionEventSelection = 1;
               std::cout<<"could not get branch : pcollisionEventSelection"<<std::endl;
               std::cout<<"set to default value : pcollisionEventSelection = "<<pcollisionEventSelection<<std::endl;
           }
       }
       else {
           pcollisionEventSelection = 0;    // default value if the collision is not HI, will not be used anyway.
       }
       Int_t pPAprimaryVertexFilter;    // this filter is used for PP.
       if (isPP) {
           treeSkim->SetBranchStatus("pPAprimaryVertexFilter",1);
           if (treeSkim->GetBranch("pPAprimaryVertexFilter")) {
               treeSkim->SetBranchAddress("pPAprimaryVertexFilter",&pPAprimaryVertexFilter);
           }
           else {   // overwrite to default
               pPAprimaryVertexFilter = 1;
               std::cout<<"could not get branch : pPAprimaryVertexFilter"<<std::endl;
               std::cout<<"set to default value : pPAprimaryVertexFilter = "<<pPAprimaryVertexFilter<<std::endl;
           }
       }
       else {
           pPAprimaryVertexFilter = 0;      // default value if the collision is not PP, will not be used anyway.
       }
       Int_t pBeamScrapingFilter;   // this filter is used for PP.
       if (isPP) {
           treeSkim->SetBranchStatus("pBeamScrapingFilter",1);
           if (treeSkim->GetBranch("pBeamScrapingFilter")) {
               treeSkim->SetBranchAddress("pBeamScrapingFilter",&pBeamScrapingFilter);
           }

void plotPFRandomCone(TString str = "forest.root", Int_t nRCPerEvent = 1, int maxEvents = -1) {

  // gStyle->SetOptStat(0000);
  // gStyle->SetOptTitle(0);

  double minEta = -1.5;
  double maxEta = 1.5;
  double minPhi = -TMath::Pi();
  double maxPhi = TMath::Pi();
  double radiusRC = 0.4;
  
  TChain *fChain = new TChain("hiEvtAnalyzer/HiTree");
  TChain *pfTree = new TChain("pfcandAnalyzer/pfTree");
  TChain *skimTree = new TChain("skimanalysis/HltTree");
  TChain *hltTree = new TChain("hltanalysis/HltTree");
  //  TFile *f = TFile::Open(str.Data());
  fChain->Add(str.Data());
  pfTree->Add(str.Data());
  skimTree->Add(str.Data());
  hltTree->Add(str.Data());
  fChain->AddFriend(pfTree);
  fChain->AddFriend(skimTree);
  fChain->AddFriend(hltTree);
  
  if(!fChain) {
    Printf("Couldn't find pfTree. Aborting!");
    return;
  }

  Int_t MinBiasTriggerBit;
  Int_t phfCoincFilter;
  
  Int_t           hiBin;
  TBranch        *b_hiBin;   //!
  fChain->SetBranchAddress("hiBin", &hiBin, &b_hiBin);
  fChain->SetBranchAddress("HLT_HIL1MinimumBiasHF1ANDExpress_v1",&MinBiasTriggerBit);
  fChain->SetBranchAddress("phfCoincFilter3",&phfCoincFilter);
  
  ForestPFs                    fPFs;              //!PFs in tree
  if (fChain->GetBranch("nPFpart"))
    fChain->SetBranchAddress("nPFpart", &fPFs.nPFpart, &fPFs.b_nPFpart);
  if (fChain->GetBranch("pfId"))
    fChain->SetBranchAddress("pfId", fPFs.pfId, &fPFs.b_pfId);
  if (fChain->GetBranch("pfPt"))
    fChain->SetBranchAddress("pfPt", fPFs.pfPt, &fPFs.b_pfPt);
  if (fChain->GetBranch("pfVsPtInitial"))
    fChain->SetBranchAddress("pfVsPtInitial", fPFs.pfVsPt, &fPFs.b_pfVsPt);
  if (fChain->GetBranch("pfEta"))
    fChain->SetBranchAddress("pfEta", fPFs.pfEta, &fPFs.b_pfEta);
  if (fChain->GetBranch("pfPhi"))
    fChain->SetBranchAddress("pfPhi", fPFs.pfPhi, &fPFs.b_pfPhi);

  //Printf("nentries: %d",(Int_t)fChain->GetEntries());

  TList *fOutput =  new TList();
  TH1::SetDefaultSumw2();
  
  TH3D *h2CentPtRCEta = new TH3D("h2CentPtRCEta","h2CentPtRCEta;centrality;p_{T,RC};#eta",100,0,100,250,-50.,200.,60,-6,6);
  fOutput->Add(h2CentPtRCEta);
  TH3D *h2CentPtRCEtaVS = new TH3D("h2CentPtRCEtaVS","h2CentPtRCEtaVS;centrality;p_{T,RC}^{VS};#eta",100,0,100,250,-50.,200.,60,-6,6);
  fOutput->Add(h2CentPtRCEtaVS);

  TH3D *h2MultPtRCEta = new TH3D("h2MultPtRCEta","h2MultPtRCEta;centrality;p_{T,RC};#eta",3000,0,6000,150,-50.,100.,60,-6,6);
  fOutput->Add(h2MultPtRCEta);
  TH3D *h2MultPtRCEtaVS = new TH3D("h2MultPtRCEtaVS","h2MultPtRCEtaVS;centrality;p_{T,RC}^{VS};#eta",3000,0,6000,150,-50.,100.,60,-6,6);
  fOutput->Add(h2MultPtRCEtaVS);
  Printf("histos defined");

  Int_t startEntry = 0;
  Int_t lastEntry = fChain->GetEntries();//100;
  Printf("events in chain: %d",lastEntry);
  if(maxEvents<lastEntry)
    lastEntry = maxEvents;
  Printf("lastEntry: %d",lastEntry);

  TRandom3 *rnd = new TRandom3();
 
  for (int j=startEntry; j<lastEntry; j++) {
    fChain->GetEntry(j);
    if(j%100==0) std::cout << "entry: "<< j << std::endl;
    //if(!MinBiasTriggerBit) continue;
    if(!phfCoincFilter) continue;
    
    double etaRC = rnd->Rndm() * (maxEta - minEta) + minEta;
    double phiRC = rnd->Rndm() * (maxPhi - minPhi) + minPhi;

    double ptRC = 0.;
    double ptRCVS = 0.;
    Int_t pfCount = 0;
    for(Int_t i = 0; i<fPFs.nPFpart; i++) {
      double pt = fPFs.pfPt[i];
      double ptVS = fPFs.pfVsPt[i];
      double phi = fPFs.pfPhi[i];
      double eta = fPFs.pfEta[i];

      double dr = deltaR(phi,phiRC,eta,etaRC);
      if(dr<radiusRC) {
        ptRC+=pt;
        ptRCVS+=ptVS;
      }
//.........这里部分代码省略.........

//.........这里部分代码省略.........
       treeHiEvt->SetBranchStatus("run",1);
       treeHiEvt->SetBranchStatus("evt",1);
       treeHiEvt->SetBranchStatus("lumi",1);
       treeHiEvt->SetBranchStatus("vz",1);
       treeHiEvt->SetBranchStatus("hiBin",1);
       treeHiEvt->SetBranchStatus("hiHF",1);
       treeHiEvt->SetBranchStatus("hiHFplus",1);
       treeHiEvt->SetBranchStatus("hiHFminus",1);
       treeHiEvt->SetBranchStatus("hiHFplusEta4",1);
       treeHiEvt->SetBranchStatus("hiHFminusEta4",1);
      //  treeHiEvt->SetBranchStatus("hiNevtPlane",1);

       float vz;
       Int_t hiBin;
       Float_t         hiEvtPlanes[29];   //[hiNevtPlane]
       treeHiEvt->SetBranchAddress("hiEvtPlanes",&hiEvtPlanes);

       treeHiEvt->SetBranchAddress("vz",&vz);
       treeHiEvt->SetBranchAddress("hiBin",&hiBin);

       // specify explicitly which branches to store, do not use wildcard
       treeSkim->SetBranchStatus("*",0);
       treeSkim->SetBranchStatus("ana_step",1);
       treeSkim->SetBranchStatus("pcollisionEventSelection",1);
       treeSkim->SetBranchStatus("pHBHENoiseFilterResultProducer",1);
       treeSkim->SetBranchStatus("HBHENoiseFilterResultRun1",1);
       treeSkim->SetBranchStatus("HBHENoiseFilterResultRun2Loose",1);
       treeSkim->SetBranchStatus("HBHENoiseFilterResultRun2Tight",1);
       treeSkim->SetBranchStatus("HBHENoiseFilterResult",1);
       treeSkim->SetBranchStatus("HBHEIsoNoiseFilterResult",1);

       Int_t pcollisionEventSelection;

       if (treeSkim->GetBranch("pcollisionEventSelection")) {
           treeSkim->SetBranchAddress("pcollisionEventSelection",&pcollisionEventSelection);
       }
       else {   // overwrite to default
           pcollisionEventSelection = 1;
           std::cout<<"could not get branch : pcollisionEventSelection"<<std::endl;
           std::cout<<"set to default value : pcollisionEventSelection = "<<pcollisionEventSelection<<std::endl;
       }

       // event information
       UInt_t run, lumis;
       ULong64_t event;
       treeEvent->SetBranchAddress("run", &run);
       treeEvent->SetBranchAddress("event", &event);
       treeEvent->SetBranchAddress("lumis", &lumis);

       TFile* output = new TFile(outputFile,"RECREATE");

       TTree *configTree = setupConfigurationTreeForWriting(configCuts);

       int centBinWidth = 200/nCentralityBins;  // number of "hiBin"s that a centrality bin covers
       int vertexBinWidth = 30/nVertexBins;     // number of "vz"s    that a vertex     bin covers
       float eventPlaneBinWidth = TMath::Pi()/nEventPlaneBins;     // number of "vz"s    that a vertex     bin covers
                                                // accepted vz range is -15 to 15.

       TTree *outputTreesHLT[nCentralityBins][nVertexBins][nEventPlaneBins];
       TTree *outputTreesEvent[nCentralityBins][nVertexBins][nEventPlaneBins];
       TTree *outputTreesJet[nCentralityBins][nVertexBins][nEventPlaneBins][nJetCollections];
       TTree *outputTreesHiEvt[nCentralityBins][nVertexBins][nEventPlaneBins];
       TTree *outputTreesSkim[nCentralityBins][nVertexBins][nEventPlaneBins];

       Long64_t nEntriesFilled[nCentralityBins][nVertexBins][nEventPlaneBins];    // number of events read for a centrality bin
       for(int i=0; i<nCentralityBins; ++i)